The mathematical mystery inside the legendary ’90s shooter Quake 3

Offtopic
1

I enjoyed this and I think some folks around here will too:

I used Gemini to throw together some Julia code so I could see it in action:

function fast_inv_sqrt(x::Float32)
       x2 = x * 0.5f0
       y  = x

       # 1. The "Bit Hack"
       # Reinterpret float bits as a 32-bit integer
       i = reinterpret(Int32, y)

       # The Magic Number: 0x5f3759df
       # This performs the logarithmic approximation
       i = Int32(0x5f3759df) - (i >> 1)

       # Reinterpret back to a float
       y = reinterpret(Float32, i)

       # 2. Newton-Raphson Iteration
       # One iteration is enough for Quake 3's precision requirements
       # Formula: y = y * (1.5 - (x/2 * y * y))
       y = y * (1.5f0 - (x2 * y * y))

       return y
end

# Usage

julia> x = 64.0f0
64.0f0

julia> result = fast_inv_sqrt(x)
0.124788396f0

julia> actual = 1.0f0 / sqrt(x)
0.125f0
2

Relevant reading with regards to strict aliasing, which is shared by Julia and C/C++ standards:
c++ - How to implement fast inverse sqrt without undefined behavior? - Stack Overflow
C++20 gets std::bit_cast to copy bits like Julia’s reinterpret. This copy could be optimized to register moves.

3

No, that’s not relevant to Julia as far as I’m aware. Julia’s reinterpret neither guarantees copied bits nor is it susceptible to UB. Its semantics are quite optimizable and distinct from C++'s behaviors.

4

It’s irrelevant to Julia by design, it’s the original C implementation that violates strict aliasing.

5

That article makes the number far more magical than it is (are they seriously reaching for the decimal 1\,597\,463\,007 and \frac{3}{2} × 2^{23} × (127 - 0.0450465)??). The actual paper by Lomont is better, but I think it can be even simpler.

When you do maths on a floating point number as an integer, you end up doing two computations at once: you’re messing with both the exponent and the significand. The most significant part of that is of course not the significand.

Looking at just the exponent maths, 0x5f3759df has a floating point exponent value of 190. Shifting the integer representation of the input float to the right by a bit divides the exponent value by two. But the value in that exponent slot is offset by 127. And 127 + 127/2 is… 190.5! So a slightly simplified real maths approximation of what that bit twiddling is doing is simply:

\exp\big(\frac{-\log(x)}{2}\big)

Which… if you look at it closely, you’ll realize is a convoluted way to spell x^{\frac{-1}{2}}… or \frac{1}{\sqrt{x}}!

But of course, that only gets the exponent right (or close). It’s guaranteed to be within a factor of 2, though, and then the newton step iteratively refines it from there.

So that explains the 0x5f... leading bits. The mantissa part is indeed quite interesting and I wonder if you could do better by bitmasking out the odd exponents.

4 Likes
6

number that puzzles and inspires experts

Anyone who finds that puzzling today is hardly an expert. As @mbauman was faster to show there’s obviously a number around that size which gets the result about right, then it’s just a refinement search to find the optimum when combined with a Newton-Raphson step. Maybe that was a bit trickier in the 90s but today you can just brute force it.