The Quake-3 fast inverse square root algorithm

How do you compute \(\frac{1}{\sqrt{x}}\) quickly and efficiently? There’s a cool piece of code found in Quake III source code attributed to the legendary game programmer John Carmack. It actually goes back to SGI, 3dfx and was first written in mid 1980s by Greg Walsh.

This code works roughly \(4\) times faster than the naive 1.0/sqrt(x), and the maximum relative error over all floating point numbers is \(0.175228\) percent. How does it work? What bit master designed such an incredible hack?

// Fast inverse square root algorithm from Quake III Arena
float Q_rsqrt(float x){
    float xhalf = 0.5f * x;
    int i = *(int*)(&x);        // get bits for floating value
    i = 0x5f3759df - (i >> 1);  // gives an initial guess y_0
    x = *(float*)&i;            // convert bits back to float
    x = x*(1.5f - xhalf * x * x);   // Newton step, repeating increases accuracy.
}

Background

This blog post assumes 32-bit architecture. In particular, the floating point representation is IEEE 754. The above code is reported to be endian-neutral.

\(32\)-bit floating-point numbers are stored in the following format:

where \(s\) is a \(1\)-bit sign(\(1\) denotes negative), \(E\) is an \(8\)-bit exponent, and \(M\) is a \(23\)-bit mantissa. The exponent is biased by \(127\) to accomodate positive and negative exponents, and the mantissa does not store the leading \(1\), so think of \(M\) as a binary number with the decimal point to the left, thus \(M\) is a value in \(I=[0,1)\). The represented value is:

\[ x = (-1)^s(1+M)2^{E-127} \]

These bits can be viewed as the floating point representation of a real number, or thinking only of bits, as an integer. Thus, \(M\) will be considered a real number in \(I\) or as an integer, depending on context. \(M\) as a real number is \(M\) divided by \(2^{23}\).

The algorithm.

The central idea is to use an iterative root solver such as Newton-Raphson and the magic constant is used to compute a good initial guess. Given the floating point value \(x > 0\), we want to compute \(\frac{1}{\sqrt{x}}\). We let:

\[ \begin{align*} y &= \frac{1}{\sqrt{x}}\\ y^2 &= \frac{1}{x} \end{align*} \]