Fast Inverse Square Root and Magic Number 0x5f3759df

This article explores the magic constant 0x5f3759df, made famous by the source code of the video game Quake III Arena, which enables the rapid computation of the inverse square root (\(1/\sqrt{x}\)) of a 32-bit floating-point number. Below, we break down how the algorithm leverages the structure of IEEE 754 floating-point numbers in binary to approximate a logarithmic transformation using integer arithmetic, and how it refines the result using Newton’s method.

The Purpose of Fast Inverse Square Root

In 3D computer graphics, calculating vector lengths and normalizing them is a fundamental operation required for lighting, physics, and camera transformations. Normalizing a vector \((x, y, z)\) requires computing:

\[\frac{1}{\sqrt{x^2 + y^2 + z^2}}\]

In the 1990s, division and square root operations were computationally expensive at the hardware level. The fast inverse square root algorithm bypasses traditional division and root extraction by using bit manipulation and integer subtraction to produce an initial guess accurate to within a few percent, requiring only a fraction of the clock cycles.

IEEE 754 Floating-Point Binary Structure

To understand why the constant works, one must look at how 32-bit single-precision floating-point numbers are represented in the IEEE 754 standard:

• 1 sign bit (\(S\))
• 8 exponent bits (\(E\))
• 23 mantissa/fraction bits (\(M\))

The mathematical value represented by these bits is:

\[x = (1 + m) \times 2^{e}\]

Where \(m = \frac{M}{2^{23}}\) and \(e = E - 127\) (using the 127 exponent bias).

Turning Floating-Point Numbers into Logarithms

If you interpret the exact same 32 bits of a float as a 32-bit integer (\(I_x\)), the integer’s value is:

\[I_x = E \times 2^{23} + M = 2^{23}(E + m) = 2^{23}(e + 127 + m)\]

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

\[\log_2(x) = \log_2((1 + m) \times 2^e) = e + \log_2(1 + m)\]

Because \(m\) lies between \(0\) and \(1\), the function \(\log_2(1 + m)\) can be approximated as a linear function \(m + \mu\), where \(\mu\) is a chosen tuning constant (\(\approx 0.0450465\)) that minimizes error. Substituting this back gives:

\[\log_2(x) \approx e + m + \mu\]

Comparing this to the integer representation \(I_x\):

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

This reveals that reinterpreting a float’s raw bits as an integer automatically computes a scaled and shifted approximation of its base-2 logarithm.

Deriving the Magic Constant

The goal is to compute \(y = \frac{1}{\sqrt{x}} = x^{-1/2}\). Taking the base-2 logarithm of both sides:

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

Using the linear relationship between the integer interpretation and the logarithm:

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

Multiplying through and solving for \(I_y\):

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

The term \(\frac{3}{2} 2^{23}(127 - \mu)\) is a fixed constant. When evaluated with the optimal error-minimizing value of \(\mu\), this decimal value converts directly in hexadecimal to:

\[\text{Constant} \approx \text{0x5f3759df}\]

The division \(-\frac{1}{2}I_x\) is computed in binary using a simple right shift (i >> 1). Thus, the core step of the algorithm:

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

performs an approximate base-2 logarithm, multiplies it by \(-0.5\) (inverting and halving the exponent), and converts it back from logarithmic space to linear floating-point space in a single integer instruction.

Error Correction with Newton-Raphson

The integer subtraction produces an initial approximation with roughly 3% to 4% error. To achieve production accuracy, the algorithm finishes by running one iteration of the Newton-Raphson method for finding roots of \(f(y) = \frac{1}{y^2} - x = 0\):

\[y_{new} = y\left(\frac{3}{2} - \frac{x}{2}y^2\right)\]

In C, this is implemented as:

y = y * ( 1.5F - ( xhalf * y * y ) );

A single iteration brings the maximum relative error down to approximately \(0.175\%\), making it sufficiently precise for computer graphics while running significantly faster than traditional hardware division routines of its era.