Worked example

As an example, the number ${\displaystyle x=0.15625}$ can be used to calculate ${\textstyle {\frac {1}{\sqrt {x}}}\approx 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 ${\displaystyle y=2.61486}$, which has an error of about 3.4%. After one single iteration of Newton's method, the final result is ${\displaystyle y=2.52549}$, an error of only 0.17%.

Avoiding undefined behavior

According to the C standard, reinterpreting a floating point value as an integer by casting then dereferencing the pointer to it is not valid (undefined behavior). Another way would be to place the floating point value in an anonymous union containing an additional 32-bit unsigned integer member, and accesses to that integer provides a bit level view of the contents of the floating point value. However, type punning through a union is also undefined behavior in C++.

# include <stdint.h> // uint32_t

float q_rsqrt(float number)
{
  union {
    float    f;
    uint32_t i;
  } conv = { .f = number };
  conv.i  = 0x5f3759df - (conv.i >> 1);
  conv.f *= 1.5F - (number * 0.5F * conv.f * conv.f);
  return conv.f;
}

In modern C++, the usual method for implementing this function's casts is through C++20's std::bit_cast. This also allows the function to work in constexpr context:

# include <bit>
# include <limits>
# include <cstdint>

constexpr float Q_rsqrt(float number) noexcept
{
  static_assert(std::numeric_limits<float>::is_iec559); // (enable only on IEEE 754)

  float const y = std::bit_cast<float>(
    0x5f3759df - (std::bit_cast<std::uint32_t>(number) >> 1));
  return y * (1.5f - (number * 0.5f * y * y));
}

Floating-point representation

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. To encode a non-zero real number ${\displaystyle x}$ as a single precision float, the first step is to write ${\displaystyle x}$ as a normalized binary number:^[23]

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

where the exponent ${\textstyle e_{x}}$ is an integer, and ${\textstyle 1.b_{1}b_{2}b_{3}\ldots }$ is the binary representation of the significand. Since the single bit before the point in the significand is always 1, it does not need be stored. The equation can be rewritten as:

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

where ${\textstyle m_{x}}$ means ${\textstyle 0.b_{1}b_{2}b_{3}\ldots }$, so ${\textstyle m_{x}\in [0,1)}$. From this form, three unsigned integers are computed:^[24]

• ${\textstyle S_{x}}$, the "sign bit", is ${\displaystyle 0}$ if ${\displaystyle x}$ is positive and ${\displaystyle 1}$ negative or zero (1 bit)
• ${\textstyle E_{x}=e_{x}+B}$ is the "biased exponent", where ${\displaystyle B=127}$ is the "exponent bias"^{[note 2]} (8 bits)
• ${\textstyle M_{x}=m_{x}\times L}$, where ${\textstyle L=2^{23}}$^{[note 3]} (23 bits)

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

As an example, consider again the number ${\displaystyle x=0.15625=0.00101_{2}}$. Normalizing ${\displaystyle x}$ yields:

${\textstyle x=+2^{-3}(1+0.25)}$

and thus, the three unsigned integer fields are:

• ${\displaystyle S=0}$
• ${\displaystyle E=-3+127=124=0111\ 1100_{2}}$
• ${\displaystyle M=0.25\times 2^{23}=2\ 097\ 152=0010\ 0000\ 0000\ 0000\ 0000\ 0000_{2}}$

these fields are packed as shown in the figure below:

Aliasing to an integer as an approximate logarithm

If ${\textstyle {\frac {1}{\sqrt {x}}}}$ were to be calculated without a computer or a calculator, a table of logarithms would be useful, together with the identity ${\textstyle \log _{b}\left({\frac {1}{\sqrt {x}}}\right)=\log _{b}\left(x^{-{\frac {1}{2}}}\right)=-{\frac {1}{2}}\log _{b}(x)}$, which is valid for every base ${\displaystyle 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 ${\displaystyle 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 ${\displaystyle m_{x}\in [0,1)}$, the logarithm on the right-hand side can be approximated by^[26]

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

where ${\displaystyle \sigma }$ is a free parameter used to tune the approximation. For example, ${\displaystyle \sigma =0}$ yields exact results at both ends of the interval, while ${\textstyle \sigma ={\frac {1}{2}}-{\frac {1+\ln(\ln(2))}{2\ln(2)}}\approx 0.0430357}$ yields the optimal approximation (the best in the sense of the uniform norm of the error). However, this value is not used by the algorithm as it does not take subsequent steps into account.

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 ${\displaystyle x}$ as an integer ${\displaystyle I_{x}}$ yields^{[note 4]}

${\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 ${\displaystyle I_{x}}$ is a scaled and shifted piecewise-linear approximation of ${\displaystyle \log _{2}(x)}$, as illustrated in the figure on the right. In other words, ${\displaystyle \log _{2}(x)}$ is approximated by

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

First approximation of the result

The calculation of ${\textstyle y={\frac {1}{\sqrt {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 ${\displaystyle x}$ and ${\displaystyle 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 ${\displaystyle 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 )={\mathtt {0x5F3759DF}}}$

from which it can be inferred that ${\displaystyle \sigma \approx 0.0450466}$. The second term, ${\displaystyle {\frac {1}{2}}I_{x}}$, is calculated by shifting the bits of ${\displaystyle I_{x}}$ one position to the right.^[27]

Newton's method

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 ${\displaystyle y}$ as the inverse square root, ${\displaystyle f(y)={\frac {1}{y^{2}}}-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, ${\displaystyle y_{n}}$ for ${\displaystyle y}$, then a better approximation ${\displaystyle y_{n+1}}$ can be calculated by taking ${\displaystyle y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}}$, where ${\displaystyle f'(y_{n})}$ is the derivative of ${\displaystyle f(y)}$ at ${\displaystyle y_{n}}$.^[28] For the above ${\displaystyle f(y)}$,

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

where ${\displaystyle f(y)={\frac {1}{y^{2}}}-x}$ and ${\displaystyle f'(y)=-{\frac {2}{y^{3}}}}$.

Treating ${\displaystyle y}$ as a floating-point number, y = y*(threehalfs - x/2*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 (${\displaystyle y_{n+1}}$) as the input of the next iteration, the algorithm causes ${\displaystyle y}$ to converge to the inverse square root.^[29] For the purposes of the Quake III engine, only one iteration was used. A second iteration remained in the code but was commented out.^[21]

Accuracy

As noted above, the approximation is very 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, and InvSqrt() gives 9.982522, making the relative difference 0.0017478, or 0.175% of the true value, 10. The absolute error only drops from then on, and the relative error stays within the same bounds across all orders of magnitude.

Magic number

It is not known precisely how the exact value for the magic number was determined. Chris Lomont developed a function to minimize approximation error by choosing the magic number ${\displaystyle R}$ over a range. He first computed the optimal constant for the linear approximation step as 0x5F37642F, close to 0x5F3759DF, but this new constant gave slightly less accuracy after one iteration of Newton's method.^[30] Lomont then searched for a constant optimal even after one and two Newton iterations and found 0x5F375A86, which is more accurate than the original at every iteration stage.^[30] He concluded by asking whether the exact value of the original constant was chosen through derivation or trial and error.^[31] Lomont said that the magic number for 64-bit IEEE754 size type double is 0x5FE6EC85E7DE30DA, but it was later shown by Matthew Robertson to be exactly 0x5FE6EB50C7B537A9.^[32]

Jan Kadlec reduced the relative error by a further factor of 2.7 by adjusting the constants in the single Newton's method iteration as well,^[33] arriving after an exhaustive search at

conv.i = 0x5F1FFFF9 - ( conv.i >> 1 );
conv.f *= 0.703952253f * ( 2.38924456f - x * conv.f * conv.f );
return conv.f;

A complete mathematical analysis for determining the magic number is now available for single-precision floating-point numbers.^[34][35]

Zero finding

Intermediate to the use of one vs. two iterations of Newton's method in terms of speed and accuracy is a single iteration of Halley's method. In this case, Halley's method is equivalent to applying Newton's method with the starting formula ${\displaystyle f(y)={\frac {1}{y^{1/2}}}-xy^{3/2}=0}$. The update step is then

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

where the implementation should calculate ${\displaystyle xy_{n}^{2}}$ only once, via a temporary variable.

Obsolescence

Subsequent additions by hardware manufacturers have made this algorithm redundant for the most part. For example, on x86, Intel introduced the SSE instruction rsqrtss in 1999. In a 2009 benchmark on the Intel Core 2, this instruction took 0.85ns per float compared to 3.54ns for the fast inverse square root algorithm, and had less error.^[1]

Some low-cost embedded systems do not have specialized square root instructions. However, manufacturers of these systems usually provide trigonometric and other math libraries, based on algorithms such as CORDIC. Most applications in this area are not performance-sensitive and a specialized fast inverse square root algorithm is not needed.

Bibliography