THE Fast inverse sqrt

(Or: How To Get Profanity Into The Comments Of Famous, Well-Publicised Code)

Why?

3D graphics programming uses loads of inverse square root calculation
Distance calculations!
Normalisation of vectors!
Which is used for lighting and shadows and basically everything ever

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

\frac{​ 1 ​ ​}{​ \sqrt ​ x ​ 2 ​ ​ + y ​ 2 ​ ​ + z ​ 2 ​ ​ ​ ​ ​ ​}

enter Fastinvsqrt...

Finding the square root is computationally expensive
Division is also computationally expensive
Approximate?

float FastInvSqrt(float x) {
  float xhalf = 0.5f * x;
  int i = *(int*)&x;         // evil floating point bit level hacking
  i = 0x5f3759df - (i >> 1);  // what the f***?
  x = *(float*)&i;
  x = x*(1.5f-(xhalf*x*x));
  return x;
}

(actual, verbatim code from the Quake III Arena source, edited only to align with the Synacor HR policy...)

wtf?

int i = *(int*)&x;

Reinterpret the floating point encoding of x as an integer!

wtf?

Floating point numbers have three parts: sign, exponent, and mantissa:

s e e e e e e e e m m m m m m m m m m m m m m m m m m m m m m m

Mantissa encodes a value between 0 and 1, exponent is a signed int between -127 and 128.

M = integer rep. of mantissa, E = integer rep. of exponent

m = \frac{M}{L}

m = \frac{​ M ​ ​}{​ L ​}

e = E- 127

e = E - 127

(for 32-bit floats, L = 2^23)

(1+m)\cdot2^e

(1 + m) \cdot 2 ​ e ​ ​

M + L\cdot E

M + L \cdot E

....

y = \frac{1}{\sqrt{x}} = x^{-\frac{1}{2}}

y = \frac{​ 1 ​ ​}{​ \sqrt ​ x ​ ​ ​ ​} = x ​ - \frac{​ 1 ​ ​}{​ 2 ​} ​ ​

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

l o g ​ 2 ​ ​ (y) = - \frac{​ 1 ​ ​}{​ 2 ​} l o g ​ 2 ​ ​ (x)

We can replace x and y with their floating point components...

log_2(1+m_y) + e_y = -\frac{1}{2}(log_2(1+m_x)+e_x)

l o g ​ 2 ​ ​ (1 + m ​ y ​ ​) + e ​ y ​ ​ = - \frac{​ 1 ​ ​}{​ 2 ​} (l o g ​ 2 ​ ​ (1 + m ​ x ​ ​) + e ​ x ​ ​)

Helpfully, for v between 0 and 1, log2(1+v) looks a lot like a straight line...

log_2(1+v) \approx v+\sigma

l o g ​ 2 ​ ​ (1 + v) \approx v + σ

m_y + \sigma + e_y \approx-\frac{1}{2}(m_x + \sigma + e_x)

m ​ y ​ ​ + σ + e ​ y ​ ​ \approx - \frac{​ 1 ​ ​}{​ 2 ​} (m ​ x ​ ​ + σ + e ​ x ​ ​)

M_y + LE_y \approx \frac{3}{2}L(B-\sigma) - \frac{1}{2}(M_x + LE_x)

M ​ y ​ ​ + L E ​ y ​ ​ \approx \frac{​ 3 ​ ​}{​ 2 ​} L (B - σ) - \frac{​ 1 ​ ​}{​ 2 ​} (M ​ x ​ ​ + L E ​ x ​ ​)

Suddenly, integer representations!

I_y \approx \frac{3}{2}L(B-\sigma) - \frac{1}{2}I_x

I ​ y ​ ​ \approx \frac{​ 3 ​ ​}{​ 2 ​} L (B - σ) - \frac{​ 1 ​ ​}{​ 2 ​} I ​ x ​ ​

Or, back in C code...

i = K - (i >> 1); // for some magic K

Looks familiar~?

For sigma = 0.0450465...

\frac{3}{2}2^{23}(127 - 0.0450465) = 1597463007

\frac{​ 3 ​ ​}{​ 2 ​} 2 ​ 23 ​ ​ (127 - 0.0450465) = 1597463007

...which, in hex, is 0x5f3759df!

newton's method

Given a decent approximation, we can do a single iteration of Newton's method to get an /unreasonably/ accurate approximation of the inverse square root.

Used to use a lookup table for the approximation, before they figured out the mad science magic number.

More information at:

Hummus and Magnets

Further explanation and pretty graphs

(seriously, check it out, excellent explanation)