Fast Inversed Squareroot Function - Quake 3

This technique gained popularity in the early days of computer graphics when processing power was limited. It's an approximation method used to calculate the reciprocal of a square root \(\displaystyle \frac 1 {\sqrt{x}}\) faster than the standard division and square root functions. This was implemented in Quake 3 game.

Bản tiếng Việt ở đây

The source code

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

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;                       // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

  return y;
}

Let's dig into the source code

The IEEE 754 Standard

According to this standard, a float number is represented as 3 parts of 32 bits: Sign bit (x1), Exponential bits (x8), Mantissa bits (x23).
Then, a positive float number $x$ is calculated by: \(x=(1+\frac M {2^{23}})*2^{E-127} \quad (*)\)

And, bit representation of $x$ would be: \(M + 2^{23}*E \quad (**)\)

Extract $x$ binary information from the \(\log_2x​\)

For $x$ small \(x \in [0,1)\), we could approximate $x$ by: \(x \approx \log_2(1+x)\), then we can claim that:

$$\log_2(1+x) = x + \mu; \quad (\mu \text{ is empirically chosen at 0.043})$$

Looking at an arbitrary value $x$, applying (*), we have:

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

Since \(\displaystyle \frac M {2^{23}} \in [0,1)\), we now apply the approximation rule:

$$\log_2(1 + \frac M {2^{23}}) = \frac M {2^{23}} + \mu$$

Then:

$$\begin{align*} \log_2 x &= \frac M {2^{23}} + \mu + E - 127\\ &= \frac M {2^{23}} + \frac {2^{23} * E} {2^{23}} + \mu - 127 \\ &= \frac 1 {2^{23}}(\underbrace{M + 2^{23}E}_{\text{bit representation (**)}}) + \mu - 127 \end{align*}$$

Calculate \(\frac 1 {\sqrt{x}}\) using \(\log_2​\)

The meaning of the line of code: -(i >>1)

$$\log_2(\frac 1 {\sqrt{x}}) = \log_2(x^{-\frac 1 2}) = -\frac 1 2 \log_2(x)$$

Let \(\Gamma=\frac 1 {\sqrt{x}} \implies \log_2(\Gamma) = \log_2 (\frac 1 {\sqrt{x}}) = -\frac 1 2 \log_2(x)​\)

We replace the logarithm with the bit representation for both sides:

$$\begin{align*} \frac 1 {2^{23}}(M_{\Gamma} + 2^{23}*E_{\Gamma}) + \mu - 127 = -\frac1 2 \left( \frac 1 {2^{23}}(M_x + 2^{23}*E_x) + \mu -127 \right) \\ \iff M_{\Gamma} + 2^{23}E_{\Gamma} = \underbrace{\frac {3} {2} 2^{23}(127-\mu)}_{\text{0x5f3759df}} - \underbrace{\frac 1 2(\underbrace{M_x + 2^{23}*E_x)}{x}}_{i>>1} \end{align*}$$

Here reveals the secret of this line of code:

i  = 0x5f3759df - ( i >> 1 );  // what the fuck?

Newton Iteration Method for Finding Approximation of Root

$$\boxed{x_{new} = x_{old} - \frac {f(x)} {f' (x)}}$$

Let \(\displaystyle f(y)=\frac 1 {y^2} - x\), then: \(\displaystyle f(y)=0 \iff y=\frac 1 {\sqrt{x}}\). We try to solve \(f(y)=0\) using Newton-Ralphson method. First we compute the derivative of $f(y)$:

$$f'(y) = \frac {y-xy^3} {-2}$$

Then

$$y_{new} = y - \frac {f(y)}{f' (y)} = y - \frac {y-xy^3} {-2} = y(\frac 3 2 - \frac x 2 y^2)$$

This is the meaning of the last line of code:

 y  = y * ( threehalfs - ( x2 * y * y ) );

Demo & JS Code

JS Code

function fastInverseSqrt(x) {
            const threeHalfs = 1.5;

            let i = new Float32Array(1);
            let y = new Float32Array(1);

            y[0] = x;

            i = new Int32Array(y.buffer); // Reinterprets the bits as an integer
            i[0] = 0x5f3759df - (i[0] >> 1); // Magic number and bit shift

            y = new Float32Array(i.buffer); // Reinterpret the bits back to float
            y[0] = y[0] * (threeHalfs - (x * 0.5 * y[0] * y[0])); // First iteration of Newton's method
            y[0] = y[0] * (threeHalfs - (x * 0.5 * y[0] * y[0])); 

            return y[0];
        }