Can you access the bits that represent an IEEE754 single precision float in a C++14 constant expression (constexpr)?<\/p>\n
(Why would you want to do that? Maybe you want to run a fast inverse square root<\/a> at compile time. Or maybe you want to do something that is actually useful. I wanted to know if it could be done.)<\/p>\n For context: this article is based on experiences using gcc-5.3.0 and clang-3.7.1 with -std=c++14 -march=native on a Sandy Bridge Intel i7. Where I reference sections from the C++ standard, I’m referring to the November 2014 draft<\/a>.<\/p>\n Before going\u00a0further, I’ll quote 5.20.6 from the standard:<\/p>\n Since this International Standard imposes no restrictions on the accuracy of floating-point operations, it is unspecified whether the evaluation of a floating-point expression during translation yields the same result as the evaluation of the same expression (or the same operations on the same values) during program execution.88 <\/sup><\/p>\n 88) Nonetheless, implementations are encouraged to provide consistent results, irrespective of whether the evaluation was performed during translation and\/or during program execution.<\/p><\/blockquote>\n In this post, I document things that worked (and didn’t work) for me. You may have a different experience.<\/p>\n (Error text from g++-5.3.0)<\/p>\n You can’t access the bits of a float via a typecast pointer [which is undefined behavior, and covered by 5.20.2.5]:<\/p>\n You can’t convert it via a reinterpret cast [5.20.2.13]<\/p>\n (gcc reports an error with the memory access, but does not object to the You can’t convert it through a union [gcc, for example,\u00a0permits this for non-constant expressions<\/a>, but the standard forbids it in 5.20.2.8]:<\/p>\n You can’t use And you can’t define a constexpr So what can you do?<\/p>\n For\u00a0 So I wrote a function ( When And we can have it:<\/p>\n Nothing difficult about that.<\/p>\n In IEEE754 land, negative zero is a thing. Ideally, we’d like\u00a0this behavior:<\/p>\n But the check for zero also matches negative zero.\u00a0Negative zero is not something that the C++ standard has anything to say about, given that IEEE754 is\u00a0an implementation choice [3.9.1.8: “The value representation of floating-point types is implementation defined”]. My compilers treat negative zero the same as zero for all comparisons and arithmetic operations. As such,\u00a0 I did look into other methods\u00a0for detecting negative zero in C++, without finding something that would work in a constant expression. I have seen divide by zero used as a way to detect negative zero (resulting in\u00a0\u00b1infinity, depending on the sign of the zero), but that doesn’t compile in a constant expression:<\/p>\n and divide by zero is explicitly named as undefined behavior in 5.6.4, and so by 5.20.2.5 is unusable in a constant expression.<\/p>\n Result: negative zero becomes positive zero.<\/p>\n We want this:<\/p>\n And this:<\/p>\n works as expected.<\/p>\n Same idea, different sign:<\/p>\n Also works.<\/p>\n There’s no way to generate arbitrary NaN constants in a constant expression that I can see (not least because casting bits to floats isn’t possible in a constant expression, either), so it seems impossible to get this right in general.<\/p>\n In practice, maybe this is good enough:<\/p>\n NaN values can be anywhere in the range of It is possible to detect a NaN value in a constant expression, but not its payload. (At least that I’ve been able to find). So there’s this:<\/p>\n Which means that of the 2*(223<\/sup>-1) possible NaNs, one will be handled correctly (in this case, So… partial success?\u00a0NaNs are correctly detected, but the bits for only one NaN value will be returned correctly.<\/p>\n (On the other hand, the probability that it will ever matter could be stored\u00a0as a denormalized float)<\/p>\n [That It is possible to detect and correctly construct bits for every normalized value. It does requires a little care to avoid truncation and undefined behavior. I wrote a few different implementations — the one that I describe here requires relatively little code, and doesn’t perform terribly [0].<\/p>\n The first step is to find and clear the sign bit. This simplifies subsequent steps.<\/p>\n Now we have What happens when a float\u00a0is cast\u00a0to an integral type?<\/p>\n The value of If truncation takes place, data is lost.\u00a0If the number is too large, the result is (probably) meaningless.<\/p>\n For our conversion function, if we can scale Single precision IEEE 754 floating point numbers have, at most, (23+1) bits of precision (23 in the significand, 1 implicit).\u00a0This means that we can scale down large numbers and scale up small numbers into the required range.<\/p>\n Multiplying by powers of two change only the exponent of the float, and leave the significand unmodified. As such, we can arbitrarily scale a float by a power of two and — so long as we don’t over- or under-flow the float — we will not lose any of the bits in the significand.<\/p>\n For the sake of simplicity (believe it or not [1]), my approach is to scale If Next, All of the bits of adjust the exponent accordingly:<\/p>\n and construct the result:<\/p>\n With this, we have correct results for every normalized float.<\/p>\n The final detail is denormalized values. Handling of normalized values as presented so far fails because denormals will have additional leading zeroes. They are fairly easy to account for:<\/p>\n To attempt to demystify that And that’s it. All the possible values of a float are accounted for.<\/p>\n (Side note: If you’re compiling with -ffast-math, passing denormalized numbers to You can indeed obtain the bit representation of a floating point number at compile time. Mostly. Negative zero is wrong, NaNs are detected but otherwise not accurately converted.<\/p>\n Enjoy your compile-time bit-twiddling!<\/p>\n The whole deal:<\/p>\n [0] Why does runtime performance matter? Because that’s how I tested the conversion function while implementing it. I was applying\u00a0Bruce Dawson’s advice for testing floats<\/a>\u00a0and the quicker I found out that I’d broken the conversion the better. For the implementation described in this post, it takes about 97 seconds to test all four billion float values on my laptop — half that time if I wasn’t testing negative numbers (which are unlikely to cause problems due to the way I handle the sign bit). The implementation I’ve described in this post is not the fastest solution to the problem, but it is relatively compact, and well behaved in the face of Admission buried in a footnote: I have not validated correct behavior of this code for every floating point number in actual compile-time constant expressions. Compile-time evaluation of four billion invocations of [1] It is conceptually simpler to multiply [2] Update 2016-03-28: Thanks to satbyy<\/a>\u00a0for pointing out the missing ampersand — it was lost sometime after copying the code into the article.<\/p>\n [3] Update 2016-03-28: Thanks to louiswins<\/a> for pointing out additional code errors.<\/p>\n","protected":false},"excerpt":{"rendered":" Can you access the bits that represent an IEEE754 single precision float in a C++14 constant expression (constexpr)? (Why would you want to do that? Maybe you want to run a fast inverse square root at compile time. Or maybe you want to do something that is actually useful. I wanted to know if it … Continue reading “floats, bits, and constant expressions”<\/span><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[5,26],"tags":[],"_links":{"self":[{"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/posts\/1518"}],"collection":[{"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1518"}],"version-history":[{"count":172,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/posts\/1518\/revisions"}],"predecessor-version":[{"id":1691,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=\/wp\/v2\/posts\/1518\/revisions\/1691"}],"wp:attachment":[{"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1518"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1518"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/brnz.org\/hbr\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1518"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}Methods of conversion that won’t work<\/h3>\n
constexpr uint32_t bits_cast(float f)\r\n {\r\n return *(uint32_t*)&f; \/\/ [2]\r\n }\r\n<\/pre>\nerror:\r\n accessing value of 'f' through a 'uint32_t {aka unsigned int}' glvalue\r\n in a constant expression<\/pre>\nconstexpr uint32_t bits_reinterpret_cast(float f) \r\n{ \r\n \u00a0const unsigned char* cf = reinterpret_cast<const unsigned char*>(&f); \r\n \/\/ endianness notwithstanding\r\n \u00a0return (cf[3] << 24) | (cf[2] << 16) | (cf[1] << 8) | cf[0]; \r\n} \r\n<\/pre>\nerror:\r\n '*(cf + 3u)' is not a constant expression\r\n<\/pre>\n
reinterpret_cast<\/code>. clang produces a specific error for the cast.)<\/p>\nconstexpr uint32_t bits_union(float f) \r\n{ \r\n union Convert { \r\n uint32_t u;\r\n float f;\r\n constexpr Convert(float f_) : f(f_) {}\r\n };\r\n return Convert(f).u;\r\n}\r\n<\/pre>\nerror:\r\n accessing 'bits_union(float)::Convert::u' member instead of \r\n initialized 'bits_union(float)::Convert::f' member \r\n in constant expression\r\n<\/pre>\n
memcpy()<\/code> [5.20.2.2]:<\/p>\nconstexpr uint32_t bits_memcpy(float f) \r\n{\r\n uint32_t u = 0;\r\n memcpy(&u, &f, sizeof f);\r\n return u;\r\n}\r\n<\/pre>\nerror:\r\n 'memcpy(((void*)(&u)), ((const void*)(&f)), 4ul)' \r\n is not a constant expression\r\n<\/pre>\n
memcpy()<\/code>-like function that is capable of the task [5.20.2.11]:<\/p>\nconstexpr void* memcpy(void* dest, const void* src, size_t n)\r\n{\r\n char* d = (char*)dest;\r\n const char* s = (const char*)src;\r\n while(n-- > 0)\r\n *d++ = *s++;\r\n return dest;\r\n}\r\n\r\nconstexpr uint32_t bits_memcpy(float f)\r\n{\r\n uint32_t u = 0;\r\n memcpy(&u, &f, sizeof f);\r\n return u;\r\n}\r\n<\/pre>\nerror:\r\n accessing value of 'u' through a 'char' glvalue\r\n in a constant expression\r\n<\/pre>\n
Floating point operations in constant expressions<\/h3>\n
constexpr float f = 2.0f, g = 2.0f<\/code>\u00a0the following operations are available [as they are not ruled out\u00a0by anything I can see in 5.20]:<\/p>\n\n
\nstatic_assert(f == g, \"not equal\");<\/code><\/li>\n
\nstatic_assert(f * 2.0f == 4.0f, \"arithmetic failed\");<\/code><\/li>\n
\nconstexpr int i = (int)2.0f; static_assert(i == 2, \"conversion failed\");<\/code><\/li>\n<\/ul>\nuint32_t bits(float)<\/code>) that will return the binary representation of an IEEE754 single precision float. The full function is at the end of this post. I’ll go through the various steps required to produce (my best approximation of) the desired result.<\/p>\n1. Zero<\/h3>\n
bits()<\/code> is passed the value zero, we want this behavior:<\/p>\nstatic_assert(bits(0.0f) == 0x00000000);<\/pre>\n
if (f == 0.0f)\r\n return 0;\r\n<\/pre>\n
2. Negative zero<\/h3>\n
static_assert(bits(-0.0f) == 0x80000000)<\/pre>\n
bits()<\/code>\u00a0returns the wrong value when considering negative zero, returning 0x00000000<\/code> rather than the desired\u00a00x80000000<\/code>.<\/p>\nconstexpr float r = 1.0f \/ -0.0f;<\/pre>\n
error: '(1.0e+0f \/ -0.0f)' is not a constant expression\r\n<\/pre>\n
3. Infinity<\/h3>\n
static_assert(bits(INFINITY) == 0x7f800000);<\/pre>\n
else if (f == INFINITY)\r\n return 0x7f800000;\r\n<\/pre>\n
4. Negative Infinity<\/h3>\n
static_assert(bits(-INFINITY) == 0xff800000);<\/pre>\n
else if (f == -INFINITY)\r\n return 0xff800000;\r\n<\/pre>\n
5. NaNs<\/h3>\n
static_assert(bits(NAN) == 0x7fc00000);<\/pre>\n
0x7f800001 -- 0x7fffffff<\/code> and 0xff800001 -- 0xffffffff<\/code>. I have no idea as to the specific values that are seen in practice, nor what they mean. 0x7fc00000<\/code> shows up in \/usr\/include\/bits\/nan.h<\/code> on the system I’m using to write this, so — right or wrong — I’ve chosen that as the reference value.<\/p>\n else if (f != f) \/\/ NaN\r\n return 0x7fc00000; \/\/ This is my NaN...\r\n<\/pre>\n
0x7fc00000<\/code>). For the other 16,777,213 values, the wrong value will be returned (in this case, 0x7fc00000<\/code>).<\/p>\n6. Normalized Values<\/h3>\n
\/\/ pseudo-code\r\nstatic_assert(bits({ 0x1p-126f, ..., 0x1.ffff7p127})\r\n == { 0x00800000, ..., 0x7f7fffff});\r\nstatic_assert(bits({ -0x1p-126f, ..., -0x1.ffff7p127})\r\n == { 0x80800000, ..., 0xff7fffff});\r\n<\/pre>\n0x1pnnn<\/em>f<\/code> format<\/a>\u00a0happens to be a convenient way to represent exact values that can be stored as binary floating point numbers]<\/p>\n bool sign = f < 0.0f; \r\n float abs_f = sign ? -f : f;\r\n<\/pre>\n
abs_f<\/code> — it’s positive, non-zero, non-infinite, and not a NaN.<\/p>\nuint64_t i = (uint64_t)f;\r\n<\/pre>\n
f<\/code> will be stored in i<\/code>, according to the following rules:<\/p>\n\n
f<\/code> is too large to be represented as a uint64_t<\/code>\u00a0(i.e. f<\/code> > 264<\/sup>-1) the result is undefined.<\/li>\n<\/ul>\nabs_f<\/code>\u00a0into a range where it is not larger than (264<\/sup>-1), and it has no fractional part, we have access to an exact representation of the bits that make up the float. We just need to keep track of the amount of scaling being done.<\/p>\nabs_f<\/code>\u00a0in steps of 241<\/sup> so that (abs_f<\/code>\u00a0\u2265 287<\/sup>) like so:<\/p>\n int exponent = 254; \r\n\r\n while(abs_f < 0x1p87f) \r\n { \r\n abs_f *= 0x1p41f; \r\n exponent -= 41; \r\n }\r\n<\/pre>\nabs_f<\/code> \u2265\u00a0287<\/sup>, the least significant bit of abs_f<\/code>, if set, is 2(87-23)<\/sup>==264.<\/sup><\/p>\nabs_f<\/code> is scaled back down by 264<\/sup>\u00a0(which adds no fractional part as the least significant bit is 264<\/sup>) and converted to an unsigned 64 bit integer.<\/p>\n uint64_t a = (uint64_t)(abs_f * 0x1p-64f);\r\n<\/pre>\n
abs_f<\/code> are now present in a<\/code>, without overflow or truncation. All that is needed now is to determine where they are:<\/p>\n int lz = count_leading_zeroes(a);\r\n<\/pre>\n
exponent -= lz;\r\n<\/pre>\n
uint32_t significand = (a << (lz + 1)) >> (64 - 23); \/\/ [3]\r\n return (sign << 31) | (exponent << 23) | significand;\r\n<\/pre>\n
7. Denormalized Values<\/h3>\n
\/\/ pseudo-code\r\nstatic_assert(bits({ 0x1.0p-149f, ..., 0x1.ffff7p-127f})\r\n == { 0x00000001, ..., 0x007fffff});\r\nstatic_assert(bits({ -0x1.0p-149f, ..., -0x1.ffff7p-127f})\r\n == { 0x80000001, ..., 0x807fffff});\r\n<\/pre>\n if (exponent <= 0) \r\n { \r\n exponent = 0; \r\n lz = 8 - 1; \r\n }\r\n<\/pre>\nlz = 8 - 1<\/code>\u00a0a little: there are 8 leading bits that aren’t part of the significand of a denormalized single precision float after the repeated 2-41<\/sup> scaling that has taken place. There is also no leading 1 bit that is present in all normalized numbers (which is accounted for in the calculation of significand<\/code> above as (lz + 1)<\/code>). So the leading zero count (lz<\/code>) is set to account for the 8 bits of offset to the start of the denormalized significand, minus the one that the subsequent calculation assumes it needs to skip over.<\/p>\nbits()<\/code> will return invalid results. That’s -ffast-math for you. With gcc or clang, you could add an #ifdef __FAST_MATH__<\/code> around the test for negative exponent.)<\/p>\nConclusion<\/h3>\n
\n\/\/ Based on code from \r\n\/\/ https:\/\/graphics.stanford.edu\/~seander\/bithacks.html\r\nconstexpr int count_leading_zeroes(uint64_t v) \r\n{ \r\n constexpr char bit_position[64] = { \r\n 0, 1, 2, 7, 3, 13, 8, 19, 4, 25, 14, 28, 9, 34, 20, 40, \r\n 5, 17, 26, 38, 15, 46, 29, 48, 10, 31, 35, 54, 21, 50, 41, 57, \r\n 63, 6, 12, 18, 24, 27, 33, 39, 16, 37, 45, 47, 30, 53, 49, 56, \r\n 62, 11, 23, 32, 36, 44, 52, 55, 61, 22, 43, 51, 60, 42, 59, 58 }; \r\n \r\n v |= v >> 1; \/\/ first round down to one less than a power of 2 \r\n v |= v >> 2; \r\n v |= v >> 4; \r\n v |= v >> 8; \r\n v |= v >> 16; \r\n v |= v >> 32; \r\n v = (v >> 1) + 1; \r\n \r\n return 63 - bit_position[(v * 0x0218a392cd3d5dbf)>>58]; \/\/ [3]\r\n}\r\n \r\nconstexpr uint32_t bits(float f) \r\n{ \r\n if (f == 0.0f) \r\n return 0; \/\/ also matches -0.0f and gives wrong result \r\n else if (f == INFINITY) \r\n return 0x7f800000; \r\n else if (f == -INFINITY) \r\n return 0xff800000; \r\n else if (f != f) \/\/ NaN \r\n return 0x7fc00000; \/\/ This is my NaN...\r\n \r\n bool sign = f < 0.0f; \r\n float abs_f = sign ? -f : f; \r\n \r\n int exponent = 254; \r\n \r\n while(abs_f < 0x1p87f) \r\n { \r\n abs_f *= 0x1p41f; \r\n exponent -= 41; \r\n } \r\n \r\n uint64_t a = (uint64_t)(abs_f * 0x1p-64f); \r\n int lz = count_leading_zeroes(a);\r\n exponent -= lz;\r\n \r\n if (exponent <= 0) \r\n { \r\n exponent = 0; \r\n lz = 8 - 1;\r\n } \r\n \r\n uint32_t significand = (a << (lz + 1)) >> (64 - 23); \/\/ [3]\r\n return (sign << 31) | (exponent << 23) | significand; \r\n}\r\n<\/pre>\n-ffast-math<\/code>.<\/p>\nbits()<\/code>\u00a0takes more time than I’ve been willing to invest so far.<\/p>\nabs_f<\/code> by two (or one half) until the result is exactly positioned so that no leading zero count is required after the cast — at least, that was what I did in my first attempt. The approach described here was found to be significantly faster. I have no doubt that better-performing constant-expression-friendly approaches\u00a0exist.<\/p>\n