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。

64 位的情况又如何？

我们概括一下之前 32 位的分析，我们主要探索了以下几个问题：

什么是 Newton 迭代？
这段代码前面的 shift 部分的大概直觉，确认指数部分的取值。
理解最后那一行确实在做 Newton 迭代。
筛选尾数部分，找到在 minimax 意义下的最优解。

而 Newton 迭代公式本身并不会因为从 binary32 变成 binary64 就有所改变，因此问题 1 和 3 无需重新探索，只用看 2 和 4 即可。

对于问题 2，solution 是 immediate 的，因为公式

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

仍然成立，只不过 bias 从 127 变成了 1023，经过计算后得知指数位为 1534。

对于问题 4，此时 magic number 的搜索空间从原来的 $2^{23}$ 变成了 $2^{52}$ ，即使使用之前的粗筛+精筛法，成本也基本不可行。但注意到我们之前在 binary32 中发现 magic number 和误差的函数关系是一个 V 字形，是一种单峰函数。除了最 naive 的遍历之外，我也立刻想出了三分的做法——那是因为五年前我校测机试的 OI 部分也考到了这个，那真是一段令人怀念的岁月啊。

当时那段回忆录我写：

“回到2题，我抛弃了数学方法并尝试三分，仍然因为精度问题爆零，我气愤地把所有的变量全部改为double，事情便有了进展，获得了8分的好成绩。”

这里那道题的满分是 10 分。当年我用 double 来优化当时的三分算法，现在我又用三分来优化 FISR 在 double 情况下的 magic number，挺美妙的。

回到正题，三分的思路本身很简单，对于有极小值的单峰函数 $f$ ，初始区间 $[a, b]$ ，令两个三等分点 $m_{1} = \frac{2}{3} a + \frac{1}{3} b$ ， $m_{2} = \frac{1}{3} a + \frac{2}{3} b$ ：

若 $f (m_{1}) > f (m_{2})$ ，则更新 $a = m_{1}$ ；
若 $f (m_{1}) \leq f (m_{2})$ ，则更新 $b = m_{2}$ ；
为了 robustness，当 $b - a$ 小于某个设定的距离 $Δ$ 后，就可以恢复穷举了。

对于本文的 case， $f$ 的 evaluate 方式仍然是固定随机种子 0xC0FFEE，计算 10000 个随机正规 binary64 值的 minimax。选取的 $Δ = 512$ ，并分开评估 0 次和一次 Newton 迭代的情形。实现见附录 D：三分搜索。初步实验结果如下：

优化目标	论文结果	我们的随机搜索结果	$dist / 2^{52}$	$ε_{m a x}$
0 次 Newton	`0x5FE6EC85E7DE30DA`	`0x5FE6EC853A893173`	$6.457 \times 10^{- 7}$	3.421%
1 次 Newton	`0x5FE6EB50C7B537A9`	`0x5FE6EB50BC36EFCB`	$4.282 \times 10^{- 8}$	0.175%

这里的对照组是论文 Improving the Accuracy of the Fast Inverse Square Root by Modifying Newton–Raphson Corrections。对照组跟 Lomont 的论文一样，采用了详细的不等式推导进行了极小值的求解。观察如下：

虽然我们选出的 magic number 与论文并非完全相同，但在 binary32 尺度上也就是 1 bit 以内的差距。
两种 magic number 算出来的结果，与真实值的 $ε_{m a x}$ 的差距也在 $10^{- 5}$ 量级以内。

当然，为了尽量规避先前提到的 selection bias，10000 个随机正规 binary64 值肯定不够，样本量仍然需要扩大。但遍历 $[1, 4)$ 之间的所有正规浮点数仍然极不现实，于是我们同时采用了三种方案：

$[1, 4)$ 内的等距 $2^{24}$ 输入。
$[1, 4)$ 内的随机 $2^{24}$ 输入。
直接将先前 binary32 测的 $[1, 4)$ 转换到 binary64 作为输入。

具体实现见附录 E，最终结果如下：

Newton	方法	Magic number	$dist / 2^{52}$	$ε_{max}$	$Δ ε_{max}$ (pp)
0 次	论文	`0x5fe6ec85e7de30da`	—	3.421281331784%	—
	等距	`0x5fe6ec85e7c18a4c`	$4.169 \times 10^{- 10}$	3.421281298320%	$- 3.346 \times 10^{- 8}$
	随机	`0x5fe6ec85e7b577a0`	$5.926 \times 10^{- 10}$	3.421281284218%	$- 4.757 \times 10^{- 8}$
	float32 网格	`0x5fe6ec85e7a5ce2a`	$8.205 \times 10^{- 10}$	3.421281265925%	$- 6.586 \times 10^{- 8}$
1 次	论文	`0x5fe6eb50c7b537a9`	—	0.175118367122%	—
	等距	`0x5fe6eb50c7b33619`	$2.919 \times 10^{- 11}$	0.175118366879%	$- 2.429 \times 10^{- 10}$
	随机	`0x5fe6eb50c736bce6`	$1.841 \times 10^{- 9}$	0.175118351810%	$- 1.531 \times 10^{- 8}$
	float32 网格	`0x5fe6eb50c77c4635`	$8.286 \times 10^{- 10}$	0.175118360228%	$- 6.894 \times 10^{- 9}$

其中：

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

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

$Δ ε_{max}$ 越小越好。

另一方面，为了检验方案是否能够避免过拟合，我同时也做了 cross-evaluation：把每种采样方案选出的冠军，分别放到另外两种输入集合上测试，并同时对比论文的 baseline，结果如下图所示：

因此我们得到的最终结论是：

对于做 0 次 Newton 迭代，我们的筛选结果稳定优于论文的 magic number。
对于做 1 次 Newton 迭代，至少在等距、随机与 binary32 网格这几种 case 中，等距样本筛选出来的 magic number 稳定优于论文的 magic number。

做 1 次 Newton 迭代出现的 discrepancies，解释也比较容易：

等距的话就相当于完全公平的抛硬币正 500 反 500 情形，没有 bias；
random 就是最不均匀的，容易出现过拟合，并在等距 case 中出现明显退化；
binary32 网格由于 IEEE 754 表示本身就不 even，但也远好于随机数选择的 bias，所以效果介于两者之间。

参考代码

运行需要依赖 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。

附录 D：三分搜索

下载 appendix_d_ternary.py

运行命令：

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

--samples 10000 指定随机样本数，--seed 0xC0FFEE 固定随机种子，--bracket-width 512 指定三分搜索停止后穷举的区间宽度。结果默认保存为 fisr_binary64_ternary_results.csv。

附录 E：binary64 最终实验

下载 appendix_e_final.py

运行命令：

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

--stage all 执行三种采样方案的优化、交叉验证和绘图，--out-dir . 将 CSV 和图片保存到当前目录。每种数据集包含 $2^{24}$ 个输入；输出图片为 appendix_e_cross_0_newton.png 和 appendix_e_cross_1_newton.png。