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?

我们假设读者已经熟悉数值分析、IEEE 754 浮点数表示，以及基本的 C 语言知识。

Newton 迭代法

首先看到代码注释里面出现了两次 iteration，我们容易想到 Newton 迭代法，其具体内容如下：

假如我们要近似求解 $f (x) = 0$ ，并且有一个近似解 $x_{n}$ ，在 $(x_{n}, f (x_{n}))$ 处进行泰勒展开，解析式为

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

解 $y = 0$ 得到 $x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$ ，依次迭代下去即可接近 $f (x) = 0$ 的真正解。

当然这里有两个重要的问题：

对于初始选取的近似解 $x_{0}$ 有讲究吗？
收敛速度有多快？

对于第一个问题，答案是 yes。简单的例子形如在 $f (x) = x^{2} - 2$ 的 $x = 0$ 处做切线，你不会得到与 x 轴的任何交点，因此方法直接失效。另一个更加不平凡的例子是 $f (x) = x^{3} - 2 x + 2$ ，可以求得迭代式为 $x_{n + 1} = x_{n} - \frac{x_{n}^{3} - 2 x_{n} + 2}{3 x_{n}^{2} - 2}$ ，如果取 $x_{0} = 0$ 开始迭代，则结果会在 0、1 之间横跳，仍然失败。

所以我们有一个直觉：如果我们能取到一个离真正的解比较接近的初始值 $x_{0}$ ，那么通常情况下是比较“讲究”的。

对于第二个问题，答案是二次收敛。为了建立相关方法论的直觉，这里我们不考虑一般情况的证明，而是以 $f (x) = x^{2} - 2$ 作为具体 case 证明为什么这个函数的牛顿迭代时二次收敛的：

我们知道这个方程的正根是 $α = \sqrt{2}$ ，迭代式为：

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

首先我们证明在 $x_{n} > 0$ 的情况下一定收敛到 $\sqrt{2}$ ：首先第一次迭代，均值不等式告诉我们结果一定大于等于 $\sqrt{2}$ ，然后 $x_{n + 1} - x_{n} = \frac{2 - x_{n}^{2}}{2 x_{n}} \leq 0$ ，因此序列单调递减有下界，必定收敛。最后令极限为 $L$ ，代入 $L = \frac{L^{2} + 2}{2 L}$ 得 $L = \sqrt{2}$ ，证毕。

令误差 $e_{n} = x_{n} - α$ ，尝试表示

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

因此我们求得 $\frac{| e_{n + 1} |}{| e_{n} |^{2}} = \frac{1}{2 x_{n}}$ ，前文我们已经得到序列的极限为 $\sqrt{2}$ ，因此右边趋于一个常数，我们就证明了二次收敛。

初始值近似

然后我们回到 1st iteration 的前五行，根据上一节提到的 Newton 的局限性，我们可以合理猜想 i = 0x5f3759df - ( i >> 1 ); 这一行其实也对结果做了一个粗略但足够“讲究”的初始值近似——我们可以先写一段脚本验证几个值。

我们定义相对误差 $ε$ 为：

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

则对于一些常用的输入 $x$ ，结果如下：

输入 $x$	运算后整数 $i$	得到的 $y$	$1 / \sqrt{x}$	相对误差
$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%

可以发现结果虽然不够准确，但对于 $1 / \sqrt{x}$ 这些性质比较单调的函数应该是足够“讲究”了。所以猜想成立：这一段确实做了一个合理的初始值近似功能。

那么紧接着下一个问题就是，为什么这段位运算代码能够用来近似 $1 / \sqrt{x}$ 的结果呢？其实答案并不难想，我们考虑一个 32-bit 浮点数的 IEEE 754 表现形式：最高的一个 bit 是符号位 $S$ ，接下来 8-bit 是指数位 $E$ ，最后 23-bit 是尾数位 $M$ 。IEEE 754 形式转换成浮点形式的映射公式为：

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

我们先考虑一个简单的 case： $\sqrt{x}$ ，那么在不考虑 bias 的情况（也就是 $b i a s = 0$ ）下，直接将IEEE 754 表现形式的结果右移一位便是一个不错的估计——因为符号位不变，指数位差不多变成原来的 $1 / 2$ ，尾数位是一团 mess。但由于 $\log \sqrt{x} = \frac{1}{2} \log x$ ，因此这个估计是合理的。

而实际上， $b i a s = 127$ ，因此简单的右移一位显然不能 work（尝试代入 1，会直接数量级错误）。实际上，对 $E$ 的处理可以形式化为：

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

进一步化简后，我们发现式子可以统一为 E'=(E+127)>>1，而 127>>1 等于 63，而 0x5f3759df 对应的浮点数是 $(1 + \frac{3627487}{8388608}) 2^{63}$ 。因此我们或许可以隐约感觉到，按照这个思路，是有可能弄清楚源代码那一行的含义的。

现在我们考虑要求的式子本身。同样， $\log 1 / \sqrt{x} = - \frac{1}{2} \log x$ ，所以我们有:

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

也就是：

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

而 $⌊ \frac{3 \times 127}{2} ⌋ = 190$ ，刚好就是 0x5f3759df 的指数位！而 $\frac{1}{2} E$ 也刚好对应了后面 - ( i >> 1 ) 那一步。因为整个 IEEE 754 表示右移一位，指数位的确会变成原来的一半（如果是奇数就再去掉 0.5，但实际上也跑到 mantissa 的最高位了，所以影响不大）。

所以我们可以得出结论，0x5f3759df 这个 magic number 真正重要的是指数位的值，mantissa 最多会影响初始值的精度，进而通过 Newton 迭代影响最终取值的精度，对算法本身的正确性影响不大。

我们会在下下节有更好的尾数吗，讨论尾数如何影响精确程度的问题。

再看 Newton 迭代

接下来到了一个整篇文章最头痛的时刻——我们知道 Newton 迭代的公式是 $x_{n + 1} = x_{n} - \frac{f (x_{n})}{f^{'} (x_{n})}$ ，但实际上我们的代码是 y = y * ( threehalfs - ( x2 * y * y ) );，长得有点像，但一眼也没法直接套上去。

问题的瓶颈就出在，我们知道根是 $1 / \sqrt{x}$ ，但满足这个根的方程却有很多，那么哪一个可以对上这段代码呢？或许有人可以直接“注意到”，但至少我这次没有成功。于是我开始直接列出方程：

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

整理后得到

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

这里 $f$ 是关于 $y$ 的函数。于是我们得到了一个 ODE，但是我已经早已忘记了 ODE 是如何求解的，于是我求助了 GPT 并得知了以下方法：

对于形如
$f^{'} (y) = g (y) f (y)$
的微分方程，只要在 $f \neq 0$ 的区间上，两边除以 $f (y)$ ：
$\frac{f^{'} (y)}{f (y)} = g (y)$
注意：
$\frac{d}{d y} \ln | f (y) | = \frac{f^{'} (y)}{f (y)}$
因此直接积分：
$\ln | f (y) | = \int g (y) d y + C$
这就是一阶可分离变量微分方程最基础的求解方法。

你的方程正好属于这一类，并不需要学习高阶 ODE、常微分方程组或数值 ODE。

那么学习完之后，代入 $g (y) = \frac{2}{x y^{3} - y}$ ，即可得到：

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

对于分式积分，我还依稀记得需要拆开，先换成比较熟悉的符号：

\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

通分后也就是 $(a u + v) x^{2} + w x - u \equiv 2$ ，解方程后可得 $u = - 2, v = 2 a, w = 0$ ，代入后继续积分：

\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

把符号换回来，就是：

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

验证一下 $1 / \sqrt{x}$ 的确是 $f$ 的零点，这就证明了 y = y * ( threehalfs - ( x2 * y * y ) ); 这一行代码的确就是标准的 Newton 迭代步骤。

有更好的尾数吗？

假如我们的浮点数都是正规的，现在回到之前忽略的尾数问题。假如我们不考虑尾数的 mess，把代码改成 i = 0x5f000000 - (i >> 1);，对于一些较为极端的例子，误差结果如下：

输入 $x$	近似值 $y$	真值 $1 / \sqrt{x}$	相对误差
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%

并且可以证明，在保持 magic number 为 0x5f000000 的情况，对于任何输入 $x$ ，最大的相对误差就是 $\frac{1 - \sqrt{2}}{\sqrt{2}} \approx - 29.289 %$ 。（证明略）

如果连尾数都可以自由选取（从 0x5f000000 到 0x5f7fffff），那么全局最大误差出现在 $4 / 3$ 处，但实际上受浮点精度限制，对应的 round 值约为 1.33333337，算出来的最大误差 $ε_{m a x}$ 为（相对于 1.33333337）：

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

与前文约 $3 %$ 左右的误差还有数十倍的差距，这也不是靠一次 Newton 迭代就能完全追上的——这就来到了算法准确性的部分。

为了选取最好的 magic number，我们可以设计这样的筛选策略：

直觉告诉我，尾数 $M$ 的轻微变化应该不会让误差 turbulence，因此我们可以先分段选取样本，步长比如说为 0x1000，从 0x5f000000 到 0x5f7ff000 。
对于输入 $x$ ，随机抽 10000 个 binary32 范围内的正数并计算误差；之后的 0x1000 步长内部的遍历环节，再随机抽 10000 个 binary32 范围内的正数并固定，让这 4096 个 magic number 计算各自误差，然后选出误差最低的那几个。

然后这里就引来了一个重要问题——如何定义“最好”，或者说“误差最小”？作为写过论文的人，我当即想出了几种 metric：

geomean 最小：用于衡量整体的误差程度，只不过误差为 0 的时候需要使用极小 $ε$ 替代；
average 最小：考虑 outlier 给数据带来的影响；
minimax（最极端的误差）最小：极端 outlier 的 weight 最重。

我这里选择的第三种 metric，因为这段代码是出自 FPS 游戏 Quake III Arena。作为游戏玩家的哲学，就拿 FPS 来说，肯定是希望“平均略坏”大于“突然卡住”，所以我们要尽可能优化最坏的情况。

选择适合的优化目标后，我们开始第一步的粗筛，并尝试分析 magic number 与 minimax 值的函数关系（前文我猜测应该没有 turbulence）。为了方便复现，这里的随机数种子是 0xC0FFEE，绘图如下（完整实验代码见附录 A：Minimax 粗筛绘图程序）：

可以发现，我的直觉是正确的，图像确实形成了一个 V 形结构。进行 0 次 Newton 迭代（蓝色）和 1 次迭代（黄色），均在 0x5f376000 处分别取到了 3.427736% 和 0.176516% 的极值。

接下来是确定精筛的范围，我们可以有对 V 形函数有以下直觉：如果 $f (a) > f (b) < f (c), a < b < c$ ，那么最小值一定出现在 $(a, c)$ 之间。对于这道题，我们遍历 [0x5f375000, 0x5f377fff] 这个区间就应该足够（代码见附录 B：Minimax 精筛绘图程序）。

然后，有趣的事情发生了：我们发现 0 次迭代的最优值和 1 次迭代的最优值并不一致！前者是 0x5f37642e，后者是 0x5f375a71 和 0x5f375a73 并列第一。两者的函数图象如下：

这就意味着一个非常神奇的现象——牛顿迭代法并不保证对于两个初始误差 $e_{1} < e_{2}$ ，迭代过后的误差 $e_{1}^{'}$ 仍然小于 $e_{2}^{'}$ 。具体原因是什么，可能就涉及到诸如函数本身的性质，以及之前提到过的初始值的“讲究”程度了，这里讲起来比较复杂，文章就跳过了。

当然，另外一个问题是，为什么这两个结果都和原始 magic number 0x5f3759df 虽然很接近，但又不完全一致呢？原因其实是比较 mundane 的：样本量不够，导致选取的样本集合有某种统计意义不可避免的 bias（感性理解就是你投 1000 次币，不太可能刚好正面 500 次反面 500 次，总得偏向某一边几次十几次一样），导致精筛的时候出现针对该样本的“局部最优”。

解决办法也很暴力：覆盖所有不同的相对误差模式。我这里选取的是 $[1, 4)$ （为什么？）的全部 $2^{24}$ 个正规输入并比较 minimax（代码见附录 C：区间 [1,4) 最优求解）。实验组是我们前两个脚本算出来的“最优值”，对照组是 Chris Lomont 的那篇 Fast Inverse Square Root 论文：

Magic number	0 次 Newton	1 次 Newton
Quake III `0x5f3759df`	3.437577282%	0.175233867%
`0x5f37642e`（为 0 次优化）	3.421289021%	0.177588172%
Lomont `0x5f37642f`（同上）	3.421283763%	0.177588946%
`0x5f375a87`（为 1 次优化）	3.436540281%	0.175128778%
Lomont `0x5f375a86`（同上）	3.436546454%	0.175130156%

我们可以发现以下几点（pp = percentage point，百分点）：

我和 Lomont 选择的 0 次 Newton 迭代的 magic number，比 Quake III 0x5f3759df 这个 baseline 小 0.016 pp，但在 1 次迭代后比它误差大 0.002 pp。
我和 Lomont 选择的 1 次 Newton 迭代的 magic number，对比 baseline，实际上无论是 0 次迭代还是 1 次，误差都小 0.001 pp 和 0.0001 pp。
至于我和 Lomont 的比较，则是 0 次 Newton 迭代常数 0x5f37642f 略优于我的 0x5f37642e；而 1 次 Newton 迭代常数，我的 0x5f375a87 又略优于 0x5f375a86——两边都只是 1 bit 的差距！考虑到浮点数舍入本身也有 $10^{- 7}$ 量级的误差，基本上都可以算是不分伯仲的“最优”了。

所以回到这一节的标题，答案是：至少从我们假想的 Quake III 性能优化哲学来看，yes。

参考代码

运行需要依赖 numpy 和 matplotlib：

python3 -m pip install numpy matplotlib

建议使用 Linux 环境，精度高达 $1.084 \times 10^{- 19}$ 。某些 Windows 环境（包括本机 Windows 11）的 long double 可能与 double 无异，从而影响最后几位精度。

附录 A：Minimax 粗筛绘图

下载 appendix_a_coarse.py

运行命令：

python appendix_a_coarse.py

无命令行参数；脚本固定使用 10,000 个随机正规 binary32 输入、种子 0xC0FFEE 和步长 0x1000。输出图片为脚本同目录下的 appendix_a_coarse.png。

附录 B：Minimax 精筛绘图

下载 appendix_b_fine.py

运行命令：

python appendix_b_fine.py --output-dir .

--output-dir . 指定在当前目录保存 CSV 和图片；包括 fisr_minimax_fine_0_newton.png、fisr_minimax_fine_1_newton.png 和 fisr_minimax_fine_comparison.png。

附录 C：区间 [1,4) 最优求解

下载 appendix_c_verify.py

运行命令：

python appendix_c_verify.py --chunk-size 262144

--chunk-size 262144 指定每批验证的输入数量；完整遍历 $[1, 4)$ 的 $2^{24}$ 个 binary32 输入，默认将 CSV 保存为当前目录下的 fisr_full_range_results_python.csv。