Newton's Method & Square Roots — Quadratic Convergence and the Fast Inverse Sqrt
Newton's method (Newton-Raphson) is one of the most elegant and powerful algorithms in numerical mathematics. Starting from a tangent-line approximation, it refines a root estimate by squaring the number of correct digits with each iteration — quadratic convergence. Applied to the square root function it becomes the ancient Babylonian method. Applied to the inverse square root with a bit-level hack, it becomes the legendary algorithm that powered Quake III Arena's lighting engine.
1. Newton-Raphson Derivation
To find a root of f(x) = 0: start at a guess x₀, draw the tangent line to the curve at (x₀, f(x₀)), and find where it crosses zero.
2. Quadratic Convergence Proof
3. Square Root via Newton's Method
4. The Babylonian Method
The iteration x_{n+1} = (x + S/x)/2 for square roots was known to the Babylonians around 1700 BCE and is recorded on tablet YBC 7289. The tablet shows √2 ≈ 1.41421296... computed to 6 sexagesimal places (60-base).
Heron of Alexandria (c. 50 CE) described the same algorithm and is sometimes credited with it in Western tradition. The method is also called "Heron's method." It's the oldest known example of an iterative numerical algorithm.
5. The Quake III Fast Inverse Square Root
In 1999, id Software's Quake III Arena contained a function
Q_rsqrt that computed 1/√x ≈ 1.6% accurate in just a few
operations — no division, no sqrt call. It was needed for lighting
normal calculations (normalizing vectors requires dividing by their
length).
// Original Quake III source code (C), circa 1999
// Attributed to John Carmack; magic constant likely from Cleve Moler / Greg Walsh
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long *)&y; // bit-reinterpret float as integer
i = 0x5f3759df - (i >> 1); // magic constant initial guess
y = *(float *)&i; // bit-reinterpret integer back to float
y = y * (threehalfs - (x2 * y * y)); // ONE Newton-Raphson iteration
return y;
}
The key insight: IEEE 754 floating point represents a number as (1+m)×2^e in memory. The bit representation is (e + 127) in the high bits and m in the low bits. So log₂(x) ≈ i/L − 127 where i is the integer bit pattern and L = 2²³. The right-shift (i >> 1) approximates dividing log₂(x) by 2 — which in normal space divides the exponent by 2, i.e., takes square root. Subtracting from the magic constant flips the result to 1/√x. One Newton-Raphson iteration (for f(y) = 1/y² − x = 0) gives y = y(3/2 − x/2 · y²), refining to ~1.6% error.
0x5f3759df remains one of the most
famous numbers in game development history.
6. When Newton's Method Fails
- Multiple roots: Near a double root f(r) = f'(r) = 0, convergence degrades to linear (not quadratic). Use modified Newton's: x_{n+1} = xₙ − m·f(xₙ)/f'(xₙ) with multiplicity m to restore quadratic.
- Oscillation: Some functions produce cycles x₀ → x₁ → x₀ → ... Newton maps periodic cycles are a gateway to chaos. Example: f(x) = x³ − 2x + 2 has a 2-cycle for x₀ = 0.
- Divergence: If f'(xₙ) is near zero, the next iterate can be very far away. Always requires a good initial guess; global convergence is not guaranteed.
- Newton fractals: For complex roots of polynomials, the basins of attraction for different roots form beautiful fractal boundaries — small changes in x₀ can converge to completely different roots.