Is Carmack / Wales Inverse Square Root Algorithm Biased

Question

When implementing the "Carmack Inverse Square Root" algorithm, I noticed that the results seemed biased. The following code seems to give the best results:

float InvSqrtF(float x) { // Initial approximation by Greg Walsh. int i = * ( int* ) &x; i = 0x5f3759df - ( i >> 1 ); float y = * ( float * ) &i; // Two iterations of Newton-Raphson method to refine the initial estimate. x *= 0.5f; float f = 1.5F; y = y * ( f - ( x * y * y ) ); y = y * ( f - ( x * y * y ) ); * ( int * )(&y) += 0x13; // More magic. return y; } 

The key difference is the penultimate “more magical” line. Since the original results were too low with a fairly constant coefficient, this adds 19 * 2 ^ (exponent (y) -bias) to the result with only one instruction. It seems like they give me about 3 extra bits, but I don’t notice something?

Answer 1

Newton's method causes an offset. The function whose zero you need to find

 f(y) = x - 1/y² 

is concave, therefore - unless you start with y ≥ √(3/x) - the Newton method only produces approximations ≤ 1/√x (and strictly less, unless you start with an exact result) with exact arithmetic.

Floating-point arithmetic sometimes causes too large approximations, but usually not in the first two iterations (since the original assumption is usually not close enough).

So yes, there is bias, and adding a small amount usually improves the result. But not always. For example, in the region of about 1.25 or 0.85, the results without adjustment are better than with. In other regions, tuning gives one bit of extra precision, while in others it gives more.

In any case, the added magic constant should be adjusted to the area from which x most often taken for the best results.

Answer 2

Since this method is an approximation, the result will be somewhat overstated and underestimated by some others. You can find some good numbers on McEniry paper about how this error propagates for different configurations and the math behind them.

So, if you have no convincing evidence that the result is clearly biased in your application area, I would prefer to adjust the magic constant as suggested in the Lomont document :-)