April 20, 2023

In 2005, id Software released the source code for their 1999 game Quake III
Arena under the GPL-2 license. In the file code/game/q_math.c,
there is a function for calculating the reciprocal square root of a number
which at first glance seems to use a very peculiar algorithm:

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

Many articles have been written about this particular algorithm and it has its
own well written Wikipedia page where it is referred to as the
fast inverse square root. The algorithm actually appeared on various forums
before the Q3 source code was released. Ryszard of Beyond3D did some
investigating in 2004-2005 and eventually tracked down
the original author of the algorithm to be Greg Walsh at Ardent Computer who
created it more than a decade earlier.

How does it work?

So how does the method work, anyway? It is performed in two steps:

1. obtain a rough approximation y for the reciprocal square root of our
   number:

   y  = number;
   i  = * ( long * ) &y;
   i  = 0x5f3759df - ( i >> 1 );
   y  = * ( float * ) &i;
   

2. improve the approximation using a single step of the Newton-Raphson (NR) method:

   const float threehalfs = 1.5F;
   x2 = number * 0.5F;
   y  = y * ( threehalfs - ( x2 * y * y ) );
   

First approximation

The most interesting part is the first one. It uses a seemingly magic number
0x5f3759df and some bit shifting and somehow ends up with the reciprocal
square root. The first line stores the 32-bit floating-point number y as a
32-bit integer i by taking a pointer to y, converting it to a long
pointer and dereferencing it. So y and i hold two identical 32-bit vectors,
but one is interpreted as a floating-point number and the other is interpreted
as an integer number. Then, the integer number is shifted one step to the
right, negated, and the constant 0x5f3759df is added. Finally, the resulting
value is interpreted as a floating number again by dereferencing a float
pointer that points to the integer i value.

Here, shifting, negation and addition is performed in the integer domain, how
do these operations affect the number in the floating-point domain? In order to
understand how this can yield an approximation of the reciprocal square root we
must be familiar with how floating point numbers are represented in memory. A
floating-point number consists of a sign $$s in {0,1}$$, exponent $$ein mathbb{Z}$$ and a fractional part $$0leq{f}<1$$. The value of the
floating-point numbe