C++ implementation of Java's floatToIntBits() and intBitsToFloat()

Question

I am trying to implement Java's floatToIntBits() and intBitsToFloat() methods in C++. The latter method is the inverse of the former method, and the purpose of the former one is to pack a 32-bit floating point number into a portable 32-bit integer format. More specifically, my replacement function float2ul() for floatToIntBits() aims at the followings:

• It should be portable across platforms that have 8-bit char types and support IEEE floating point arithmetic.
• It returns a representation of the specified floating-point value according to the IEEE 754 floating-point "single format" bit layout.
• Bit 31 (the bit that is selected by the mask 0x80000000UL) represents the sign of the floating-point number (0 means positive and 1 means negative). Bits 30-23 (the bits that are selected by the mask 0x7f800000UL) represent the exponent. Bits 22-0 (the bits that are selected by the mask 0x007fffffUL) represent the significand of the floating-point number.
• If the argument is positive infinity, the result is 0x7f800000UL.
• If the argument is negative infinity, the result is 0xff800000UL.
• If the argument is NaN, the result is 0x7fc00000UL. C++ does not seem to offer a standard way to distinguish a quiet NaN from a signaling NaN, so all NaNs here are encoded with the same value.
• In all cases, if the float type has exactly 32 bits, the result is a unsigned long that, when given to my replacement function ul2float() for Java's intBitsToFloat() method, will produce a floating-point value the same as the argument to ul2float (except all NaN values are collapsed to the single "canonical" NaN value 0x7fc00000UL).
• If the size of the float type on the target platform is greater than 32 bits, the caller is responsible for checking whether the number is within the range of 32-bit floating point finite numbers. If x falls inside the range of 32-bit floating point numbers but its precision is high to fit in 32 bits, the statement ul2float(float2ul(x)) - x == 0.0f is not necessarily true but the difference should be small.

Here is my implementation:

bigpack.h

#ifndef BIGPACK_H
#define BIGPACK_H

#include <limits>
#include <climits>

static_assert(CHAR_BIT==8,
"The char type on this platform has more than 8 bits. \
packnumbers.h requires 8-bit char, i.e. CHAR_BIT==8.");

static_assert(
std::numeric_limits<double>::is_iec559 &&
std::numeric_limits<long double>::is_iec559,
"This platform does not comply to IEC 559 (IEEE 754-1985) \
floating-point arithmetic.");

// The following replacement functions are needed because
// a float type can have more than 32 bits. The function
// isnan() in <cmath> is still usable and no replacement is
// needed.
float FLT32_max();  // replacement of numeric_limits<float>::max()
bool FLT32_isfinite(float x);  // replacement of isfinite()
bool FLT32_isinf(float x);  // replace of isinf()

unsigned long float2ul(float x);  // equivalent of Java's floatToIntBits()
float ul2float(unsigned long i);  // equivalent of Java's intBitsToFloat()

#endif

bigpack.cpp

/*
References:
1. W. Kahan (1997), Lecture Notes on the Status of
   IEEE Standard 754 for Binary Floating-Point Arithmetic
   (http://www.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF),

2. Documentations of the methods floatToIntBits(), intBitsToFloat()
   doubleToLongBits() and longBitsToDouble() in Java
   (http://docs.oracle.com/javase/6/docs/api/java/lang/Float.html)

3. Brian "Beej Jorgensen" Hall (2012),
   Beej's Guide to Network Programming
   (http://beej.us/guide/bgnet/output/html/singlepage/bgnet.html),
   sec. 7.4.
*/

#include <cmath>
#include "bigpack.h"

// Some helper functions for computing powers of 2
template<int N> double _p() { return 2.0 * _p<N-1>(); }
template<> double _p<0>() { return 1.0; }
static const double t2t2[10]  // t2t2[j] = 2^{2^j}
    = {_p<1>(), _p<2>(), _p<4>(), _p<8>(), _p<16>(),
    _p<32>(), _p<64>(), _p<128>(), _p<256>(), _p<512>()};

// Given -126 <= j <= 127, return 2^j
double two_to_power(int j)
{
    if (j<0) return 1.0 / two_to_power(-j);
    double result = 1.0;
    for (int b=0; b<10 && j>0; ++b, j>>=1)
        if ((j&1) > 0)
            result *= t2t2[b];
    return result;
}

static const int FLT32_WIDTH = 32;  // N+K+1 in Kahan
static const int FLT32_K = 7;    // the K in Kahan = no. of exponent bits - 1
static const int FLT32_N1 = 23;  // N-1 in Kahan = no. of significand bits
static const unsigned long FLT32_POS_ZERO = 0UL;
static const unsigned long FLT32_NEG_ZERO = 0x80000000UL;
static const unsigned long FLT32_POS_INF  = 0x7f800000UL;
static const unsigned long FLT32_NEG_INF  = 0xff800000UL;
static const unsigned long FLT32_NAN      = 0x7fc00000UL;
static const unsigned long UINT32_2_TO_N1 = 1UL<<FLT32_N1;  // 2^{N-1}

static const float FLT32_EPSILON = 1.0f/UINT32_2_TO_N1;  // 2^{1-N}
static const float FLT32_MIN_NORMAL = two_to_power(-126);  // 2^{2 - 2^K}
static const float FLT32_MIN_SUBNORMAL = FLT32_MIN_NORMAL * FLT32_EPSILON;  // 2^{3 - 2^K - N}
static const float FLT32_MAX_FINITE  // (1 - 2^{-N}) 2^{2^K}
    = 2.0f * (2.0f - FLT32_EPSILON) / FLT32_MIN_NORMAL;

float FLT32_max() { return FLT32_MAX_FINITE; }
bool FLT32_isfinite(float x) { return fabs(x) <= FLT32_MAX_FINITE; }
bool FLT32_isinf(float x) { return !isnan(x) && !FLT32_isfinite(x); }
bool FLT32_isnormal(float x) { return FLT32_isfinite(x) && fabs(x) >= FLT32_MIN_NORMAL; }

unsigned long float2ul(float x)
{
    using namespace std;
    if (isnan(x)) return FLT32_NAN;
    if (FLT32_isinf(x)) return x<0.0f ? FLT32_NEG_INF : FLT32_POS_INF;
    if (x==0.0f)  return signbit(x) ? FLT32_NEG_ZERO : FLT32_POS_ZERO;

    if (!FLT32_isnormal(x)) {
        // (cf. Kahan 1997) The subnormal number x is encoded as
        // [sign bit] [000...000] [n]
        // where the positive integer n < 2^{N-1} is given by
        // n = 2^{2^K - 2 + N-1} |x|
        unsigned long n = fabs(x) / FLT32_MIN_SUBNORMAL;
        return (x<0.0f ? FLT32_NEG_ZERO : 0UL) | n;
    }

    int e;
    float s = fabs(frexp(x, &e));
    // (cf. Kahan 1997) With the above e and s,
    // the normal number x will be encoded as
    // [sign bit] [ 2^K+k-1 ] [ n-2^{N-1} ]
    // where
    // k = e-1,
    // n = 2^N * s
    return
    (x<0.0f ? FLT32_NEG_ZERO : 0UL)  // sign bit
    | ( ((unsigned long)(e + (1<<FLT32_K) - 2)) << FLT32_N1 )  // 2^K+k-1 shifted
    | (unsigned long)( UINT32_2_TO_N1 * (2.0f * s - 1.0f) );  // n-2^{N-1}
}

float ul2float(unsigned long i)
{
    using namespace std;
    if (i==FLT32_NAN) return numeric_limits<float>::quiet_NaN();
    if (i==FLT32_NEG_INF) return -numeric_limits<float>::infinity();
    if (i==FLT32_POS_INF) return numeric_limits<float>::infinity();
    if (i==FLT32_NEG_ZERO) return -0.0f;
    if (i==FLT32_POS_ZERO) return 0.0f;

    if ((i & FLT32_POS_INF) == 0UL) {
        // (cf. Kahan 1997) The subnormal number x is encoded as
        // [sign bit] [000...000] [n]
        // where the positive integer n < 2^{N-1} is given by
        // n = 2^{2^K - 2 + N-1} |x|
        float x = (i & ~FLT32_NEG_INF) * FLT32_MIN_SUBNORMAL;
        return (i & FLT32_NEG_ZERO) > 0UL ? -x : x;
    }

    // (cf. Kahan 1997) The normal number x encoded as
    // [sign bit] [ 2^K+k-1 ] [ n-2^{N-1} ]
    // is given by
    // x = 2^{k + 1-N} n = 2^{k-1} * (2n) * 2^{1-N}
    unsigned long n = (i & ~FLT32_NEG_INF) + UINT32_2_TO_N1;
    long k1 = (long)((i & FLT32_POS_INF)>>FLT32_N1) - (1L<<(FLT32_K)); // k-1
    float x = two_to_power(k1) * (2.0f * FLT32_EPSILON * n);
    return (i & FLT32_NEG_ZERO) > 0UL ? -x : x;
}

My concerns are:

1. I really don't like the present way of calculating powers of 2, but some other existing ways, I believe, are even worse. For instance, in Beej's Guide to Network Programming, the powers of 2 are essentially obtained by repeated multiplications/divisions of 2, which I think is too inefficient. The standard library function pow() provides another option, but the standard does not guarantee that pow() will give an accurate result when the base is 2. One viable alternative that I can think of is table lookup. Since we only need 127 powers of 2 here, precomputing all of them will not take up much space. However, if I also want to write a double2ull(double) function, I will need 1023 powers of 2 and the table would seem too large.
2. Since I use a template function to generate powers of 2, I have trouble moving some useful static constants (FLT32_MIN_NORMAL, FLT32_MIN_SUBNORMAL and FLT32_MAX_FINITE) to the header file without contaminating it with helper functions.
3. The current implement also makes use of the standard library function frexp(). Unlike the case with pow(), I believe that the use of frexp() should not be a problem, but I am not 100% sure about that.
4. The biggest concern, of course, is how portable is the above pieces of code.

Answer 1

Since most C++ platforms do use IEEE 754 this is much more easily done using a reinterpret_cast. But we do have to make sure it is portable.

Your code basically takes a floating point number and breaks it apart into its representational parts than builds an IEEE 754 format object that it pretends is an int. So we can still use that on systems that don't natively support IEEE 754.

template<bool> std::uint32_t float2ul(float x);
template std::uint32_t float2ul<true>(float x)  {return *reinterpret_cast<std::uint32_t*>(&x);}
template std::uint32_t float2ul<false>(float x) {/* Fancy rebuild here */}


std::uint32_t float2ul(float x)
{
    return float2ul<std::numeric_limits<float>::is_iec559>(x);  // Compile time choice.
}

Note rather than use unisgned long (which may technically be larger than 32 bits) I would use the specific size of integer you want (32 bit). Which is now standardized as std::uint32_t

Your use of variables names. They are a little on the short side. Make them longer and more descriptive.

To answer your questions:

I really don't like the present way of calculating powers of 2, but some other existing ways, I believe, are even worse. For instance, in Beej's Guide to Network Programming, the powers of 2 are essentially obtained by repeated multiplications/divisions of 2, which I think is too inefficient.

Why do you think it is efficient? On integers this is a left shift. On Floats this is incrementing the exponent by 1 (assuming IEEE754 representation). Both relatively trivial operations.

The standard library function pow() provides another option, but the standard does not guarantee that pow() will give an accurate result when the base is 2.

That's because the double can only provide so much accuracy. Your solution suffers from the same problem. Once you get past the number of significant bits that represent the mantissa you will start to loose accuracy.

One viable alternative that I can think of is table lookup. Since we only need 127 powers of 2 here, precomputing all of them will not take up much space. However, if I also want to write a double2ull(double) function, I will need 1023 powers of 2 and the table would seem too large.

1K of table space is not that big on modern machines. I would not worry two much about it.

Since I use a template function to generate powers of 2, I have trouble moving some useful static constants (FLT32_MIN_NORMAL, FLT32_MIN_SUBNORMAL and FLT32_MAX_FINITE) to the header file without contaminating it with helper functions.

You do have to worry about recursion limits in template functions. Different compilers provide different limits. The minimum they must support is 6 (the last time I looked). But your compiler can blow up. You way want to pre-compute the value and save it to a file that you read in.

The current implement also makes use of the standard library function frexp(). Unlike the case with pow(), I believe that the use of frexp() should not be a problem, but I am not 100% sure about that.

If your library is standards compliant I would not worry about (just like I would not worry about pow()).

The biggest concern, of course, is how portable is the above pieces of code.

I am assuming your manually implementing the IEEE754 representation. If that is the case I don't see any problems (as long as your way of transporting integers works and takes account of endianess).

Answer 2

I can't really speak to the portability of such code, as I honestly have no experience with anything other than IEEE754. Since you are targeting C++11, I can suggest that you can replace your compile-time calculation of powers of two in an easier way, using constexpr:

constexpr double pow_two(int i, double d=1)
{
    return (i == 0) ? d : pow_two(i-1, d * 2.0);
}

Which when called like so:

constexpr auto x = pow_two(7);

will be evaluated at compile time.

Also, avoid names like _p (starting with _), as the rules around what is and is not reserved for compiler usage are arcane, and its best to just avoid prefixing anything with _.

Answer 3

Why not use memcpy? gcc and clang are both smart enough to turn this into a virtual no-op, i.e. just transferring the data from one register kind to another. This is a standards-compliant way to do it.

int ieee = /*...*/;
float result;
memcpy(&result, &ieee, sizeof(float));
return result;