Fast inverse square root

As an example, the number x = 0.15625 can be used to calculate ​¹⁄_√x ≈ 2.52982. The first steps of the algorithm are illustrated below:

0011_1110_0010_0000_0000_0000_0000_0000  Bit pattern of both x and i
0001_1111_0001_0000_0000_0000_0000_0000  Shift right one position: (i >> 1)
0101_1111_0011_0111_0101_1001_1101_1111  The magic number 0x5F3759DF
0100_0000_0010_0111_0101_1001_1101_1111  The result of 0x5F3759DF - (i >> 1)

Interpreting as IEEE 32-bit representation:

0_01111100_01000000000000000000000  1.25 × 2−3
0_00111110_00100000000000000000000  1.125 × 2−65
0_10111110_01101110101100111011111  1.432430... × 263
0_10000000_01001110101100111011111  1.307430... × 21

Reinterpreting this last bit pattern as a floating point number gives the approximation y = 2.61486, which has an error of about 3.4%. After one single iteration of Newton's method, the final result is y = 2.52549, an error of only 0.17%.

Since this algorithm relies heavily on the bit-level representation of single-precision floating-point numbers, a short overview of this representation is provided here. In order to encode a non-zero real number x as a single precision float, the first step is to write x as a normalized binary number:[15]

${\displaystyle {\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\\&=\pm 2^{e_{x}}(1+m_{x})\end{aligned}}}$

where the exponent e_x is an integer, m_x ∈ [0, 1), and 1.b₁b₂b₃... is the binary representation of the "significand" (1 + m_x). Since the single bit before the point in the significand is always 1, it need not be stored. From this form, three unsigned integers are computed:[16]

• S_x, the "sign bit", is 0 if x > 0, and 1 if x < 0 (1 bit)
• E_x = e_x + B is the "biased exponent", where B = 127 is the "exponent bias"[note 3] (8 bits)
• M_x = m_x × L, where L = 2²³[note 4] (23 bits)

These fields are then packed, left to right, into a 32-bit container.[17]

As an example, consider again the number x = 0.15625 = 0.00101₂. Normalizing x yields:

x = +2⁻³(1 + 0.25)

and thus, the three unsigned integer fields are:

• S = 0
• E = −3 + 127 = 124 = 01111100₂
• M = 0.25 × 2²³ = 2097152 = 01000000000000000000000₂

these fields are packed as shown in the figure below:

If ​¹⁄_√x was to be calculated without a computer or a calculator, a table of logarithms would be useful, together with the identity log_b(​¹⁄_√x) = −1/2 log_b(x), which is valid for every base b. The fast inverse square root is based on this identity, and on the fact that aliasing a float32 to an integer gives a rough approximation of its logarithm. Here is how:

If x is a positive normal number:

${\displaystyle x=2^{e_{x}}(1+m_{x})}$

then

${\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x})}$

and since m_x ∈ [0, 1), the logarithm on the right hand side can be approximated by[18]

${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }$

where σ is a free parameter used to tune the approximation. For example, σ = 0 yields exact results at both ends of the interval, while σ ≈ 0.0430357 yields the optimal approximation (the best in the sense of the uniform norm of the error).

The integer aliased to a floating point number (in blue), compared to a scaled and shifted logarithm (in gray).

Thus there is the approximation

${\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .}$

Interpreting the floating-point bit-pattern of x as an integer I_x yields[note 5]

${\displaystyle {\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&=L(e_{x}+m_{x}+\sigma +B-\sigma )\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}}$

It then appears that I_x is a scaled and shifted piecewise-linear approximation of log₂(x), as illustrated in the figure on the right. In other words, log₂(x) is approximated by

${\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).}$

The calculation of y = ​¹⁄_√x is based on the identity

${\displaystyle \log _{2}(y)=-{\tfrac {1}{2}}\log _{2}(x)}$

Using the approximation of the logarithm above, applied to both x and y, the above equation gives:

${\displaystyle {\frac {I_{y}}{L}}-(B-\sigma )\approx -{\frac {1}{2}}\left({\frac {I_{x}}{L}}-(B-\sigma )\right)}$

Thus, an approximation of I_y is:

${\displaystyle I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}}$

which is written in the code as

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

The first term above is the magic number

${\displaystyle {\tfrac {3}{2}}L(B-\sigma )={\text{0x5F3759DF}}}$

from which it can be inferred that σ ≈ 0.0450466. The second term, 1/2I_x, is calculated by shifting the bits of I_x one position to the right.[19]

Relative error between direct calculation and fast inverse square root carrying out 0, 1, 2, 3, and 4 iterations of Newton's root-finding method. Note that double precision is adopted and the smallest representable difference between two double precision numbers is reached after carrying out 4 iterations.

With y as the inverse square root, f (y) = 1/y² − x = 0. The approximation yielded by the earlier steps can be refined by using a root-finding method, a method that finds the zero of a function. The algorithm uses Newton's method: if there is an approximation, y_n for y, then a better approximation y_n+1 can be calculated by taking y_n − f (y_n)/f ′(y_n), where f ′(y_n) is the derivative of f (y) at y_n.[20] For the above f (y),

${\displaystyle y_{n+1}={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}}$

where f (y) = 1/y² − x and f ′(y) = −2/y³.

Treating y as a floating-point number, y = y*(threehalfs - x2*y*y); is equivalent to

${\displaystyle y_{n+1}=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right)={\frac {y_{n}\left(3-xy_{n}^{2}\right)}{2}}.}$

By repeating this step, using the output of the function (y_n+1) as the input of the next iteration, the algorithm causes y to converge to the inverse square root.[21] For the purposes of the Quake III engine, only one iteration was used. A second iteration remained in the code but was commented out.[13]

A graph showing the difference between the heuristic fast inverse square root and the direct inversion of square root supplied by libstdc. (Note the log scale on both axes.)

As noted above, the approximation is surprisingly accurate. The single graph on the right plots the error of the function (that is, the error of the approximation after it has been improved by running one iteration of Newton's method), for inputs starting at 0.01, where the standard library gives 10.0 as a result, while InvSqrt() gives 9.982522, making the difference 0.017479, or 0.175% of the true value, 10. The absolute error only drops from then on, while the relative error stays within the same bounds across all orders of magnitude.