A fast method to round a double to a 32-bit int explained
Solution 1
A value of the double
floating-point type is represented like so:
and it can be seen as two 32-bit integers; now, the int
taken in all the versions of your code (supposing it’s a 32-bit int
) is the one on the right in the figure, so what you are doing in the end is just taking the lowest 32 bits of mantissa.
Now, to the magic number; as you correctly stated, 6755399441055744 is 251 + 252; adding such a number forces the double
to go into the “sweet range” between 252 and 253, which, as explained by Wikipedia, has an interesting property:
Between 252 = 4,503,599,627,370,496 and 253 = 9,007,199,254,740,992, the representable numbers are exactly the integers.
This follows from the fact that the mantissa is 52 bits wide.
The other interesting fact about adding 251 + 252 is that it affects the mantissa only in the two highest bits—which are discarded anyway, since we are taking only its lowest 32 bits.
Last but not least: the sign.
IEEE 754 floating point uses a magnitude and sign representation, while integers on “normal” machines use 2’s complement arithmetic; how is this handled here?
We talked only about positive integers; now suppose we are dealing with a negative number in the range representable by a 32-bit int
, so less (in absolute value) than (−231 + 1); call it −a. Such a number is obviously made positive by adding the magic number, and the resulting value is 252 + 251 + (−a).
Now, what do we get if we interpret the mantissa in 2’s complement representation? It must be the result of 2’s complement sum of (252 + 251) and (−a). Again, the first term affects only the upper two bits, what remains in the bits 0–50 is the 2’s complement representation of (−a) (again, minus the upper two bits).
Since reduction of a 2’s complement number to a smaller width is done just by cutting away the extra bits on the left, taking the lower 32 bits gives us correctly (−a) in 32-bit, 2’s complement arithmetic.
Solution 2
This kind of "trick" comes from older x86 processors, using the 8087 intructions/interface for floating point. On these machines, there's an instruction for converting floating point to integer "fist", but it uses the current fp rounding mode. Unfortunately, the C spec requires that fp->int conversions truncate towards zero, while all other fp operations round to nearest, so doing an
fp->int conversion requires first changing the fp rounding mode, then doing a fist, then restoring the fp rounding mode.
Now on the original 8086/8087, this wasn't too bad, but on later processors that started to get super-scalar and out-of-order execution, altering the fp rounding mode generally serializes the CPU core and is quite expensive. So on a CPU like a Pentium-III or Pentium-IV, this overall cost is quite high -- a normal fp->int conversion is 10x or more expensive than this add+store+load trick.
On x86-64, however, floating point is done with the xmm instructions, and the cost of converting
fp->int is pretty small, so this "optimization" is likely slower than a normal conversion.
Related videos on Youtube
Yu Hao
GFW (Great Firewall of China) is one of the most notorious inventions in the history of Internet. Anyone working for it should be ashamed.
Updated on February 18, 2022Comments
-
Yu Hao about 2 years
When reading Lua’s source code, I noticed that Lua uses a macro to round
double
values to 32-bitint
values. The macro is defined in theLlimits.h
header file and reads as follows:union i_cast {double d; int i[2]}; #define double2int(i, d, t) \ {volatile union i_cast u; u.d = (d) + 6755399441055744.0; \ (i) = (t)u.i[ENDIANLOC];}
Here
ENDIANLOC
is defined according to endianness: 0 for little endian, 1 for big endian architectures; Lua carefully handles endianness. Thet
argument is substituted with an integer type likeint
orunsigned int
.I did a little research and found that there is a simpler format of that macro which uses the same technique:
#define double2int(i, d) \ {double t = ((d) + 6755399441055744.0); i = *((int *)(&t));}
Or, in a C++-style:
inline int double2int(double d) { d += 6755399441055744.0; return reinterpret_cast<int&>(d); }
This trick can work on any machine using IEEE 754 (which means pretty much every machine today). It works for both positive and negative numbers, and the rounding follows Banker’s Rule. (This is not surprising, since it follows IEEE 754.)
I wrote a little program to test it:
int main() { double d = -12345678.9; int i; double2int(i, d) printf("%d\n", i); return 0; }
And it outputs
-12345679
, as expected.I would like to understand how this tricky macro works in detail. The magic number
6755399441055744.0
is actually 251 + 252, or 1.5 × 252, and 1.5 in binary can be represented as 1.1. When any 32-bit integer is added to this magic number—Well, I’m lost from here. How does this trick work?
Update
-
As @Mysticial points out, this method does not limit itself to a 32-bit
int
, it can also be expanded to a 64-bitint
as long as the number is in the range of 252. (Although the macro needs some modification.) -
Some materials say this method cannot be used in Direct3D.
-
When working with Microsoft assembler for x86, there is an even faster macro written in assembly code (the following is also extracted from Lua source):
#define double2int(i,n) __asm {__asm fld n __asm fistp i}
-
There is a similar magic number for single precision numbers: 1.5 × 223.
-
Cory Nelson almost 11 years"fast" compared to what?
-
Mysticial almost 11 years@CoryNelson Fast compared to a simple cast. This method, when implemented properly (with SSE intrinsics) is quite literally a hundred times faster than a cast. (which invokes a nasty function call to a rather expensive conversion code)
-
Cory Nelson almost 11 yearsRight -- I can see it being faster than
ftoi
. But if you're talking SSE, why not just use the single instructionCVTTSD2SI
? -
Mysticial almost 11 years@CoryNelson For
double -> int32
conversions, that instruction is fine. But fordouble -> int64
, this trick becomes king. -
tmyklebu almost 11 years@Mysticial: It doesn't quite seem to work for
double -> int64
; the addition maps 2<sup>52</sup>+ and 2<sup>52</sup>+1 to the same number. (I guess it's still good for rounding things that aren't in that troublesome range, though. 50-bit numbers or something.) -
Mysticial almost 11 years@tmyklebu Many of the use cases that go
double -> int64
are indeed within the2^52
range. These are particularly common when performing integer convolutions using floating-point FFTs. -
Oak almost 11 yearsRelated: stackoverflow.com/questions/429632/…. You can also find explanation + benchmarks here.
-
MSalters almost 11 years@Mysticial: This method, when implemented properly should be just as fast as a cast, when implemented properly. Obviously your implementation poorly implements this particular cast.
-
Mysticial almost 11 years@MSalters Not necessarily true. A cast must live up to the language's specification - including proper handling of overflow and NAN cases. (or whatever the compiler specifies in the case IB or UB) These checks tend to be very expensive. The trick mentioned in this question completely ignores such corner cases. So if you want the speed and your application doesn't care (or never encounters) such corner cases, then this hack is perfectly appropriate.
-
Pascal Cuoq almost 11 yearsThis same technique is explained in this blog post, except that it is for single-precision instead of double-precision, and with the additional simplification that the number to convert is assumed positive: blog.frama-c.com/index.php?post/2013/05/04/nearbyintf3
-
brian beuning over 9 yearsIf the double is 1e-100 then it will return 0 which is right. If the double is 1e100 I think you will get the wrong answer.
-
Mysticial over 7 yearsFWIW, a similar question about efficient
double <-> int64
conversions recently came up. -
Suraj Jain about 7 yearsSir, i need to talk with you can we chat ?
-
-
YvesgereY almost 11 years"""The other interesting fact about adding 2^51+2^52 is that it affects the mantissa only in the two highest bits - which are discarded anyway, since we are taking only its lowest 32 bits""" What is that ? Adding this may shift all the mantissa !
-
Matteo Italia almost 11 years@John: of course, the whole point of adding them is to force the value to be in that range, which obviously can result in shift the mantissa (between the other things) in respect to the original value. What I was saying here is that, once you are in that range, the only bits that differ from the corresponding 53 bits integer are bit 51 and 52, which are discarded anyway.
-
Wojciech Migda almost 8 yearsFor those who'd like to convert to
int64_t
you can do that by shifting the mantissa left and then right by 13 bits. This will clear the exponent and the two bits from the 'magic' number, but will keep and propagate the sign to the entire 64-bit signed integer.union { double d; int64_t l; } magic; magic.d = input + 6755399441055744.0; magic.l <<= 13; magic.l >>= 13;
-
Kentzo almost 4 yearsDo I understand correctly that 2^51 is only needed to handle negative values?
-
rici over 3 yearsThe point of the trick described in the OP is that it is fast. Or, at least, that 20 years ago it was fast. When Pentium IV's were still common, using type punning to get an integer out of a double was a huge win (see the benchmarks here, copied from a comment to the OP.) Your version is, of course, simpler. But it doesn't solve the problem of getting that result into an integer variable without the slow cast. (There's also nothing C++-specific about your code.)
-
Amr Ali over 3 yearsI find the trick useful to round floating point numbers, and not to cast to integers.
-
Peter Cordes over 2 yearsSpecifically, x86-64 makes it cheap because SSE / SSE2 have convert-with-truncation instructions, like
cvtTsd2si
as well as convert with current rounding modecvtsd2si
. SSE3 even added x87fisttp
for truncating conversion in legacy x87 registers, useful for 32-bit code that's using 32-bit calling conventions which pass float/double in x87. And/or for functions using x87 80-bit for extra precision.