Fast Inverse Square Root Algorithm Explained

The fast inverse square root algorithm, popularized by the source code of Quake III Arena, computes \(y = 1/\sqrt{x}\) up to four times faster than traditional floating-point division and square root operations of its era. It achieves this by treating the bit representation of a standard 32-bit floating-point number as a regular integer, performing a right bit-shift and subtracting it from a specific “magic number,” and then refining the estimate using a single step of the Newton-Raphson method.

IEEE 754 Floating-Point Format in Binary

To understand how the integer trick works, one must look at how a 32-bit float (float in C) is laid out according to the IEEE 754 standard:

• 1 Sign bit (\(S\)): Represents positive (\(0\)) or negative (\(1\)). Since square roots require positive numbers, \(S = 0\).
• 8 Exponent bits (\(E\)): Biased by \(127\).
• 23 Mantissa/Significand bits (\(M\)): Represents the fractional part with an implicit leading 1.

Mathematically, the real value \(x\) represented by the float is:

\[x = \left(1 + \frac{M}{2^{23}}\right) \times 2^{E - 127}\]

When the same 32 bits are read directly as a standard 32-bit integer \(I_x\), the binary pattern is interpreted as:

\[I_x = E \times 2^{23} + M\]

Connecting Floating-Point Bits to Logarithms

Taking the base-2 logarithm of the floating-point value \(x\):

\[\log_2(x) = \log_2\left(1 + \frac{M}{2^{23}}\right) + E - 127\]

Because \(\frac{M}{2^{23}}\) falls between \(0\) and \(1\), we can use the linear approximation \(\log_2(1 + t) \approx t + \mu\), where \(\mu\) is a small tuning constant (\(\approx 0.0450465\)):

\[\log_2(x) \approx \frac{M}{2^{23}} + \mu + E - 127\]

Factoring out \(2^{23}\):

\[\log_2(x) \approx \frac{1}{2^{23}} (E \times 2^{23} + M) - 127 + \mu\]

Notice that the term \(E \times 2^{23} + M\) is identical to the integer representation \(I_x\):

\[\log_2(x) \approx \frac{I_x}{2^{23}} - 127 + \mu\]

This reveals the core insight: the raw integer value of an IEEE 754 float is approximately proportional to the logarithm of the floating-point number.

Deriving the Integer Operation for \(1/\sqrt{x}\)

We want to calculate \(y = x^{-1/2}\). Taking the logarithm of both sides:

\[\log_2(y) = -\frac{1}{2} \log_2(x)\]

Substitute the integer approximation of the logarithm for both \(y\) and \(x\):

\[\frac{I_y}{2^{23}} - 127 + \mu \approx -\frac{1}{2} \left( \frac{I_x}{2^{23}} - 127 + \mu \right)\]

Multiply through by \(2^{23}\) and solve for the integer representation of the result, \(I_y\):

\[I_y \approx \frac{3}{2} \times 2^{23} (127 - \mu) - \frac{1}{2} I_x\]

The expression \(\frac{3}{2} \times 2^{23} (127 - \mu)\) evaluates to a constant integer. When choosing an optimal \(\mu\) to minimize mathematical error across normalized floats, this constant evaluates in hexadecimal to 0x5f3759df (or 0x5f375a86 in alternate optimizations).

Dividing an integer by two is performed in hardware via a single right bit-shift (I_x >> 1). The resulting formula in C code becomes:

i = 0x5f3759df - ( i >> 1 );

Step-by-Step Breakdown in Code

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

    x2 = number * 0.5F;
    y  = number;
    
    // 1. Bit-level reinterpretation: float to integer
    i  = * ( long * ) &y;                       
    
    // 2. Logarithmic approximation using the magic constant and bit-shift
    i  = 0x5f3759df - ( i >> 1 );               
    
    // 3. Bit-level reinterpretation: integer back to float
    y  = * ( float * ) &i;                      
    
    // 4. 1st iteration of Newton-Raphson approximation
    y  = y * ( threehalfs - ( x2 * y * y ) );   

    return y;
}

1. Type Punning: The pointer cast * ( long * ) &y copies the exact 32 binary bits of the floating-point variable into a standard 32-bit integer without running floating-point conversion routines.
2. Integer Arithmetic: The right shift (i >> 1) halves the logarithm, and subtracting it from 0x5f3759df applies the inversion and exponent biasing simultaneously.
3. Restoring the Float: The bits are cast back to a floating-point structure. At this stage, the estimate is accurate to within a few percent.
4. Newton-Raphson Refinement: The formula \(y_{n+1} = y_n \left(1.5 - 0.5 x y_n^2\right)\) applies one iteration of root-finding for \(f(y) = \frac{1}{y^2} - x = 0\), bringing the final error down to less than 0.175% across the entire range of positive inputs.