Why does fast inverse square root work?

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;
}

So what the fuck?

I'll assume you're already familiar with numerical analysis, IEEE 754 floating-point representation, and basic C.

Newton's method

The comments mention "iteration" twice, so Newton's method naturally comes to mind. Here's how it works.

Suppose we want to approximately solve $f (x) = 0$ , and we already have an approximation $x_{n}$ . A Taylor expansion at $(x_{n}, f (x_{n}))$ gives the tangent line

y = f^{'} (x_{n}) (x - x_{n}) + f (x_{n})

Setting $y = 0$ gives $x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$ . Repeat this, and we can approach an actual root of $f (x) = 0$ .

Of course, two important questions remain:

Do we need to be picky about the initial guess $x_{0}$ ?
How fast does it converge?

For the first question, yes. A simple example: draw the tangent to $f (x) = x^{2} - 2$ at $x = 0$ , and it won't intersect the x-axis at all. The method fails immediately. A slightly less trivial example is $f (x) = x^{3} - 2 x + 2$ , whose iteration is $x_{n + 1} = x_{n} - \frac{x_{n}^{3} - 2 x_{n} + 2}{3 x_{n}^{2} - 2}$ . Start at $x_{0} = 0$ , and it just bounces between 0 and 1. Still no luck.

So here's some intuition: if we can pick an initial guess $x_{0}$ that's reasonably close to the actual root, that's usually a pretty decent choice.

The answer to the second question is quadratic convergence. To get a feel for why, I'll skip the general proof and use $f (x) = x^{2} - 2$ as a concrete example.

We know its positive root is $α = \sqrt{2}$ , and the iteration is

x_{n + 1} = x_{n} - \frac{x_{n}^{2} - 2}{2 x_{n}} = \frac{x_{n}^{2} + 2}{2 x_{n}}

First, let's show that it converges to $\sqrt{2}$ when $x_{n} > 0$ . After the first iteration, AM–GM tells us that the result is at least $\sqrt{2}$ . Then $x_{n + 1} - x_{n} = \frac{2 - x_{n}^{2}}{2 x_{n}} \leq 0$ , so the sequence is decreasing and bounded below. It must converge. Call the limit $L$ ; substituting it into $L = \frac{L^{2} + 2}{2 L}$ gives $L = \sqrt{2}$ . Done.

Now let the error be $e_{n} = x_{n} - α$ , and try expressing the next error:

e_{n + 1} = x_{n + 1} - α = \frac{x_{n}^{2} - 2 x_{n} α + 2}{2 x_{n}} = \frac{(x_{n} - α)^{2}}{2 x_{n}} = \frac{e_{n}^{2}}{2 x_{n}}

This gives $\frac{| e_{n + 1} |}{| e_{n} |^{2}} = \frac{1}{2 x_{n}}$ . We already know that the sequence converges to $\sqrt{2}$ , so the right-hand side approaches a constant. That's our quadratic convergence.

Approximating the initial guess

Back to the five lines before the first iteration. Given the limitations of Newton's method we just discussed, it's reasonable to guess that i = 0x5f3759df - ( i >> 1 ); produces a rough but decent initial approximation. We can write a little script to check a few values first.

Define the relative error $ε$ as

ε = \frac{y - 1 / \sqrt{x}}{1 / \sqrt{x}} \times 100 % .

For some familiar inputs $x$ , we get:

Input $x$	Integer $i$ after the operation	Resulting $y$	$1 / \sqrt{x}$	Relative error
$10^{- 12}$	`0x4970fb79`	987063.6	1000000	−1.294%
$10^{- 6}$	`0x44743e01`	976.9688	1000	−2.303%
0.001	`0x41f5d0a8`	30.72688	31.62278	−2.833%
0.01	`0x41256e5a`	10.33944	10	+3.394%
0.1	`0x4050f379`	3.264860	3.162278	+3.244%
0.5	`0x3fb759df`	1.432430	1.414214	+1.288%
1	`0x3f7759df`	0.966215	1	−3.378%
2	`0x3f3759df`	0.716215	0.707107	+1.288%
3	`0x3f1759df`	0.591215	0.577350	+2.401%
4	`0x3ef759df`	0.483108	0.5	−3.378%
10	`0x3ea759df`	0.326858	0.316228	+3.361%
100	`0x3dd359df`	0.103199	0.1	+3.199%
1000	`0x3cfa59df`	0.030560	0.031623	−3.359%
$10^{6}$	`0x3a7d47df`	0.000966189	0.001	−3.381%
$10^{12}$	`0x3582ef8d`	$9.75547 \times 10^{- 7}$	$10^{- 6}$	−2.445%
$10^{30}$	`0x2692607a`	$1.01569 \times 10^{- 15}$	$10^{- 15}$	+1.569%

The results aren't exactly accurate, but for a fairly well-behaved monotonic function like $1 / \sqrt{x}$ , they should be decent enough. So the guess checks out: this code does produce a reasonable initial approximation.

The next question, then, is why this bit-twiddling code can approximate $1 / \sqrt{x}$ at all. The idea isn't too hard to see. Consider the IEEE 754 representation of a 32-bit float: the top bit is the sign $S$ , followed by the 8-bit exponent $E$ and the 23-bit mantissa $M$ . The mapping from this representation to the floating-point value is

(- 1)^{S} \times 2^{E - b i a s} \times (1. M)

Let's start with something simpler: $\sqrt{x}$ . If we ignore the bias for now, taking $b i a s = 0$ , shifting the IEEE 754 representation right by one bit gives a fairly good estimate. The sign stays the same, the exponent is roughly halved, and the mantissa becomes a mess. But since $\log \sqrt{x} = \frac{1}{2} \log x$ , the estimate makes sense.

In reality, though, $b i a s = 127$ , so a plain right shift obviously won't work. Try plugging in 1, and you'll be off by orders of magnitude. What we actually want to do with $E$ is

E^{'} - b i a s = \frac{1}{2} (E - b i a s)

Simplify this, and we get E'=(E+127)>>1. Now, 127>>1 is 63, while the floating-point value represented by 0x5f3759df is $(1 + \frac{3627487}{8388608}) 2^{63}$ . We can start to see that this line of thought might get us somewhere with that mysterious line of code.

Now consider the expression we actually want. Similarly, $\log 1 / \sqrt{x} = - \frac{1}{2} \log x$ , so we have:

E^{'} - b i a s = - \frac{1}{2} (E - b i a s)

In other words:

E^{'} = \frac{3}{2} \times 127 - \frac{1}{2} E

And $⌊ \frac{3 \times 127}{2} ⌋ = 190$ , which is exactly the exponent field of 0x5f3759df! The $\frac{1}{2} E$ part also lines up with the - ( i >> 1 ) that follows. Shifting the entire IEEE 754 representation right by one really does halve the exponent. If it's odd, we lose another 0.5, but that bit actually ends up in the top bit of the mantissa, so it's not a huge problem.

So we can conclude that the really important part of the magic number 0x5f3759df is its exponent. The mantissa affects the accuracy of the initial guess, and therefore the final accuracy after Newton's method, but doesn't make much difference to the correctness of the algorithm itself.

We'll come back to how the mantissa affects accuracy in Can we find a better mantissa?, two sections down.

Back to Newton's method

Now for the most headache-inducing part of this whole article. We know that Newton's iteration is $x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$ , but our code says y = y * ( threehalfs - ( x2 * y * y ) );. It looks vaguely similar, but I can't just match them up at a glance.

The problem is that we know the root is $1 / \sqrt{x}$ , but plenty of equations have that root. Which one gives us this code? Maybe someone can just "notice" the answer. I couldn't, at least not this time. So I started by writing down the equation directly:

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

Rearranging gives

2 f (y) = (x y^{3} - y) f^{'} (y)

Here $f$ is a function of $y$ . Great, we've got an ODE. Except I've long since forgotten how to solve ODEs, so I asked GPT and learned the following:

For a differential equation of the form
$f^{'} (y) = g (y) f (y)$
on any interval where $f \neq 0$ , divide both sides by $f (y)$ :
$\frac{f^{'} (y)}{f (y)} = g (y)$
Note that
$\frac{d}{d y} \ln | f (y) | = \frac{f^{'} (y)}{f (y)}$
So we can integrate directly:
$\ln | f (y) | = \int g (y) d y + C$
This is the basic method for solving a first-order separable differential equation.

Your equation is exactly this kind. You don't need higher-order ODEs, systems of ODEs, or numerical ODE methods.

Lesson learned. Substitute $g (y) = \frac{2}{x y^{3} - y}$ , and we get:

\ln | f (y) | = \int \frac{2}{x y^{3} - y} d y + C

For integrals of rational functions, I vaguely remember that we're supposed to split the fraction up. Let's first switch to some more familiar notation:

\int \frac{2}{a x^{3} - x} d x = \int \frac{2}{a x^{3} - x} d x = \int (\frac{u}{x} + \frac{v x + w}{a x^{2} - 1}) d x

Putting everything over a common denominator gives $(a u + v) x^{2} + w x - u \equiv 2$ . Solving for the coefficients gives $u = - 2, v = 2 a, w = 0$ . Substitute them back and keep integrating:

\int (\frac{- 2}{x} + \frac{2 a x}{a x^{2} - 1}) d x = - 2 \ln | x | + \ln | a x^{2} - 1 | + C = \ln | a - \frac{1}{x^{2}} | + C

Switching the notation back:

f (y) = C (x - \frac{1}{y^{2}})

Check that $1 / \sqrt{x}$ is indeed a root of $f$ , and there we have it: y = y * ( threehalfs - ( x2 * y * y ) ); really is a standard Newton iteration.

Can we find a better mantissa?

Assume all our floating-point inputs are normal. Let's return to the mantissa we ignored earlier. If we throw away that whole mess and change the code to i = 0x5f000000 - (i >> 1);, here are the errors for a few fairly extreme examples:

Input $x$	Approximation $y$	True value $1 / \sqrt{x}$	Relative error
0.5	1	1.414214	−29.289%
1	0.75	1	−25.000%
2	0.5	0.707107	−29.289%
3	0.4375	0.577350	−24.223%
4	0.375	0.5	−25.000%

We can also show that, with the magic number fixed at 0x5f000000, the worst relative error over inputs $x$ is $\frac{1 - \sqrt{2}}{\sqrt{2}} \approx - 29.289 %$ . Proof omitted.

If we let the mantissa vary freely too, from 0x5f000000 to 0x5f7fffff, the global worst case occurs at $4 / 3$ . Floating-point precision rounds this to about 1.33333337, and the resulting maximum error $ε_{m a x}$ , relative to 1.33333337, is:

ε_{m a x} = \frac{1.333333253 - 0.866025391}{0.866025391} \approx + 53.960065 %

That's tens of times the roughly $3 %$ error we saw earlier. One Newton iteration isn't going to make up for all of that. This brings us to the algorithm's accuracy.

To pick the best magic number, we can try the following search strategy:

My intuition says that tiny changes in the mantissa $M$ shouldn't make the error fluctuate wildly. So we can start by sampling magic numbers at a step size of 0x1000, from 0x5f000000 to 0x5f7ff000.
For the input $x$ , randomly sample 10000 positive binary32 values and calculate the errors. Then, when searching within a 0x1000-wide interval, draw another 10000 positive binary32 values and keep them fixed. Evaluate all 4096 magic numbers on that same set and pick the few with the lowest errors.

Which brings up an important question: what exactly do we mean by "best", or "lowest error"? Having written papers before, I immediately came up with a few metrics:

Minimize the geometric mean: a measure of overall error, though zero errors need to be replaced with a tiny $ε$ .
Minimize the average: account for the effect of outliers on the data.
Minimize the minimax error, the worst error: extreme outliers carry the most weight.

I chose the third metric, because this code comes from the FPS game Quake III Arena. Speaking as a gamer, I'd rather things be slightly worse on average than suddenly freeze up. So let's optimize the worst case as much as we can.

With an objective chosen, we can start the coarse search and see how the minimax error varies with the magic number. As I guessed earlier, it shouldn't be too turbulent. For reproducibility, the random seed is 0xC0FFEE. Here's the plot; the full experiment code is in Appendix A: Coarse minimax search.

My intuition checks out: the curves do form a V shape. With 0 Newton iterations in blue and 1 iteration in yellow, both reach their minima at 0x5f376000, with errors of 3.427736% and 0.176516%, respectively.

Next we need a range for the fine search. For a V-shaped function, the intuition is that if $f (a) > f (b) < f (c), a < b < c$ , the minimum must lie somewhere in $(a, c)$ . In this case, searching [0x5f375000, 0x5f377fff] should be enough. See Appendix B: Fine minimax search for the code.

Then something interesting happens: the best constant for 0 iterations isn't the best constant for 1 iteration! The former is 0x5f37642e; for the latter, 0x5f375a71 and 0x5f375a73 tie for first place. Here are the two curves:

This means something rather curious: given two initial errors $e_{1} < e_{2}$ , Newton's method doesn't guarantee that the errors after an iteration still satisfy $e_{1}^{'} < e_{2}^{'}$ . Why, exactly? That probably comes down to the properties of the function and how "decent" our initial guesses really are. It gets a bit complicated, so I'll skip it here.

Of course, there's another question: why are both results close to the original magic number 0x5f3759df, but not exactly the same? The reason is fairly mundane. Our sample isn't large enough, so the selected set inevitably has some statistical bias. Think of tossing a coin 1000 times: you're unlikely to get exactly 500 heads and 500 tails. It'll lean one way by a handful of tosses, or a dozen or so. The fine search then finds an optimum specific to that sample.

The fix is equally brute-force: cover all the distinct relative-error patterns. I chose all $2^{24}$ normal inputs in $[1, 4)$ (why this interval?) and compared their minimax errors. The code is in Appendix C: Searching the interval [1,4). Our experimental candidates are the "best" values from the first two scripts, and the reference is Chris Lomont's Fast Inverse Square Root paper:

Magic number	0 Newton iterations	1 Newton iteration
Quake III `0x5f3759df`	3.437577282%	0.175233867%
`0x5f37642e` (optimized for 0 iterations)	3.421289021%	0.177588172%
Lomont `0x5f37642f` (same objective)	3.421283763%	0.177588946%
`0x5f375a87` (optimized for 1 iteration)	3.436540281%	0.175128778%
Lomont `0x5f375a86` (same objective)	3.436546454%	0.175130156%

A few observations (pp = percentage point):

The constants Lomont and I chose for 0 Newton iterations reduce the error by about 0.016 pp compared with the Quake III 0x5f3759df baseline, but their errors after 1 iteration are about 0.002 pp worse.
The constants Lomont and I chose for 1 Newton iteration actually improve on the baseline both before and after the iteration, by about 0.001 pp and 0.0001 pp, respectively.
Comparing my results with Lomont's, his 0-iteration constant 0x5f37642f is slightly better than my 0x5f37642e; for 1 iteration, my 0x5f375a87 is slightly better than 0x5f375a86. In both cases, it's just a 1-bit difference! Given that floating-point rounding itself introduces errors on the order of $10^{- 7}$ , I'd say these "optimal" constants are pretty much neck and neck.

So, back to this section's title: at least under our imagined Quake III optimization philosophy, yes.

What about 64 bits?

Let's recap what we explored for 32 bits:

What is Newton's method?
What's the rough intuition behind the shift at the beginning, and how do we choose the exponent?
Why is the last line really a Newton iteration?
How do we search the mantissa for the best minimax result?

The Newton iteration itself doesn't change when we move from binary32 to binary64, so there's no need to revisit 1 and 3. We only need to look at 2 and 4.

For question 2, the answer is immediate. The equation

E^{'} - b i a s = - \frac{1}{2} (E - b i a s)

still holds; the bias just changes from 127 to 1023. Work it out, and the exponent field becomes 1534.

For question 4, the magic-number search space has grown from $2^{23}$ to $2^{52}$ . Even our earlier coarse-then-fine approach is pretty much impractical. But remember that V-shaped relationship we saw between the binary32 magic number and the error? That's a unimodal function. Besides the most naive exhaustive search, ternary search immediately came to mind. Why? Because it turned up in the OI part of my university entrance programming test five years ago. Ah, those were the days.

In the memoir I wrote back then:

"Back to problem 2. I gave up on the mathematical approach and tried ternary search. Still got zero because of precision issues. In a fit of rage, I changed every variable to double. Things finally started moving, and I achieved the excellent score of 8 points."

That problem was worth 10 points. Back then, I used double to improve my ternary search. Now I'm using ternary search to improve the FISR magic number for double. Kind of beautiful, really.

Anyway, ternary search itself is simple. For a unimodal function $f$ with a minimum, start with an interval $[a, b]$ and two trisection points $m_{1} = \frac{2}{3} a + \frac{1}{3} b$ and $m_{2} = \frac{1}{3} a + \frac{2}{3} b$ :

If $f (m_{1}) > f (m_{2})$ , update $a = m_{1}$ .
If $f (m_{1}) \leq f (m_{2})$ , update $b = m_{2}$ .
For robustness, once $b - a$ is smaller than some chosen distance $Δ$ , switch back to exhaustive search.

For this article, evaluating $f$ still means calculating the minimax error on 10000 random normal binary64 values with the fixed seed 0xC0FFEE. I chose $Δ = 512$ and evaluated 0 and 1 Newton iterations separately. The implementation is in Appendix D: Ternary search. Here are the preliminary results:

Objective	Paper's result	Our random-sample search result	$dist / 2^{52}$	$ε_{m a x}$
0 Newton iterations	`0x5FE6EC85E7DE30DA`	`0x5FE6EC853A893173`	$6.457 \times 10^{- 7}$	3.421%
1 Newton iteration	`0x5FE6EB50C7B537A9`	`0x5FE6EB50BC36EFCB`	$4.282 \times 10^{- 8}$	0.175%

The reference here is Improving the Accuracy of the Fast Inverse Square Root by Modifying Newton–Raphson Corrections. Like Lomont's paper, it derives the minima through detailed inequalities. A couple of observations:

Our magic numbers aren't exactly the same as the paper's, but the differences amount to less than 1 bit at binary32 scale.
The differences in $ε_{m a x}$ between the two choices are also within the order of $10^{- 5}$ .

Of course, 10000 random normal binary64 values still aren't enough if we want to avoid the selection bias mentioned earlier. We need a larger sample. Exhaustively checking all normal floating-point values in $[1, 4)$ is still wildly impractical, so I tried three schemes:

$2^{24}$ equally spaced inputs in $[1, 4)$ .
$2^{24}$ random inputs in $[1, 4)$ .
The same binary32 grid in $[1, 4)$ we tested earlier, converted directly to binary64.

The implementation is in Appendix E. Here are the final results:

Newton iterations	Method	Magic number	$dist / 2^{52}$	$ε_{max}$	$Δ ε_{max}$ (pp)
0	Paper	`0x5fe6ec85e7de30da`	—	3.421281331784%	—
	Equally spaced	`0x5fe6ec85e7c18a4c`	$4.169 \times 10^{- 10}$	3.421281298320%	$- 3.346 \times 10^{- 8}$
	Random	`0x5fe6ec85e7b577a0`	$5.926 \times 10^{- 10}$	3.421281284218%	$- 4.757 \times 10^{- 8}$
	float32 grid	`0x5fe6ec85e7a5ce2a`	$8.205 \times 10^{- 10}$	3.421281265925%	$- 6.586 \times 10^{- 8}$
1	Paper	`0x5fe6eb50c7b537a9`	—	0.175118367122%	—
	Equally spaced	`0x5fe6eb50c7b33619`	$2.919 \times 10^{- 11}$	0.175118366879%	$- 2.429 \times 10^{- 10}$
	Random	`0x5fe6eb50c736bce6`	$1.841 \times 10^{- 9}$	0.175118351810%	$- 1.531 \times 10^{- 8}$
	float32 grid	`0x5fe6eb50c77c4635`	$8.286 \times 10^{- 10}$	0.175118360228%	$- 6.894 \times 10^{- 9}$

Here:

dist = | M_{ours} - M_{paper} |

Δ ε_{max} = ε_{max}^{ours} - ε_{max}^{paper} .

Smaller $Δ ε_{max}$ is better.

To check whether these schemes avoid overfitting, I also did some cross-evaluation: take the winner from each sampling scheme, test it on the other two input sets, and compare everything with the paper's baseline. The plots look like this:

So our final conclusions are:

With 0 Newton iterations, our selected constants consistently outperform the paper's magic number.
With 1 Newton iteration, at least across these equally spaced, random, and binary32-grid cases, the magic number selected from the equally spaced samples consistently outperforms the paper's magic number.

The discrepancies after 1 Newton iteration aren't too hard to explain:

Equally spaced inputs are like that perfectly balanced coin-toss case with 500 heads and 500 tails. No sampling imbalance.
Random sampling is the least uniform, making it prone to overfitting; it also degrades noticeably on the equally spaced test set.
The binary32 grid isn't evenly spaced because of the IEEE 754 representation itself, but it's still much better than the bias in the random sample. So it lands somewhere in between.

Reference code

You'll need NumPy and Matplotlib:

python3 -m pip install numpy matplotlib

I'd recommend running this on Linux, where the precision can reach $1.084 \times 10^{- 19}$ . In some Windows environments, including this machine's Windows 11 setup, long double may have the same precision as double, which affects the last few digits.

Appendix A: Coarse minimax search

Download appendix_a_coarse.py

Run:

python appendix_a_coarse.py

No command-line arguments. The script uses 10,000 random normal binary32 inputs, seed 0xC0FFEE, and step size 0x1000. It saves appendix_a_coarse.png next to the script.

Appendix B: Fine minimax search

Download appendix_b_fine.py

Run:

python appendix_b_fine.py --output-dir .

--output-dir . saves the CSV and plots in the current directory, including fisr_minimax_fine_0_newton.png, fisr_minimax_fine_1_newton.png, and fisr_minimax_fine_comparison.png.

Appendix C: Searching the interval [1,4)

Download appendix_c_verify.py

Run:

python appendix_c_verify.py --chunk-size 262144

--chunk-size 262144 sets the number of inputs processed per batch. The script checks all $2^{24}$ binary32 inputs in $[1, 4)$ and, by default, saves the CSV as fisr_full_range_results_python.csv in the current directory.

Appendix D: Ternary search

Download appendix_d_ternary.py

Run:

python appendix_d_ternary.py --samples 10000 --seed 0xC0FFEE --bracket-width 512

--samples 10000 sets the random sample size, --seed 0xC0FFEE fixes the random seed, and --bracket-width 512 sets the width of the interval searched exhaustively after ternary search stops. Results are saved to fisr_binary64_ternary_results.csv by default.

Appendix E: Final binary64 experiment

Download appendix_e_final.py

Run:

python appendix_e_final.py --stage all --out-dir .

--stage all runs the optimization, cross-evaluation, and plotting for all three sampling schemes. --out-dir . saves the CSV files and plots in the current directory. Each dataset contains $2^{24}$ inputs; the output plots are appendix_e_cross_0_newton.png and appendix_e_cross_1_newton.png.