Fast Inverse Square Root

The article “fast inverse sqrt” came to my attention. It shows a small function written in C which is amazingly fast and approximates sqrt(1/x) pretty well. Appearently it was used in the Quake source code to speed up vector normalizations. But how does it work? Also, can it be improved?

Once you realize that the bit pattern of a positive 32-bit IEEE 754 floating point number interpreted as integer (f2i) “satisfies”

f2i : float int (reinterpretation of bit pattern)

f2i

for some constant around it explains why the initial estimate of sqrt(1/x) can be computed the way it’s done in the following code snippet:

float InvSqrt(float x){
   float xhalf = 0.5f * x;
   int i = *(int*)&x; // store floating-point bits in integer
   i = 0x5f3759d5 - (i >> 1); // initial guess for Newton's method
   x = *(float*)&i; // convert new bits into float
   x = x*(1.5f - xhalf*x*x); // One round of Newton's method
   return x;
}

The magic line “i = 0x5f3759d5 – (i >> 1)” is due to

The accuracy of the above C implementation is best described by the possible range of relative errors:

This is a maximum relative error of 0.001735. After some testing and experimenting I derived a modified version by replacing the constants:

float InvSqrt(float x){
   uint32_t i = 0x5F1F1412 - (*(uint32_t*)&x >> 1);
   float tmp = *(float*)&i;
   return tmp * (1.69000231f - 0.714158168f * x * tmp * tmp);
}

This function returns floats with a relative error from the range

So, the maximum relative error is limited to 0.00065 which is nearly a third of the original version’s error. This improvement comes with no extra cost. I just replaced the constants.

What’s the magic behind it? Here’s a quick explanation of how I got to these coefficients: First, I tried to find another integer for the initial estimate. The goal was not to minimize the relative error but to minimize the size of the error range measured via where and are the maximum respective minimum possible relative errors. A simple bisection performed by hand led to the above constant of 0x5F1F1412. It could be a local minimum, though. There may be another slightly better choice. Note that since I didn’t minimize the relative error for the initial estimate it’s off by some factor of around 0.88. But this bias is compensated for in the next step. The other two coefficients are finally chosen to minimize the maximum relative error.

– P

Advertisement