数值计算 --- 平方根倒数快速算法(上)

平方根倒数快速算法(上) --- 向Greg Walsh致敬！

写在最前面 --- 一场关于平方根倒数快速算法作者的讨论：

上图中的这段代码出自一个早期的3D游戏<雷神之锤>的源代码，它实现的功能就是计算一个数x的平方根的倒数：

$1/\sqrt{x}$

这段代码之所以称之为经典，私以为主要是因为以下几点原因：

1，这段代码是出自于程序员大神John Carmack之手(实际上是另有其人：Greg Walsh)。

2，这段代码之简洁，且计算速度非常快。

3，这段代码中John Carmack遗留下来的那句锐评论“what the fuck”。这里，我把他翻译成“这是什么鬼？！！！”

4，这代码中神秘的magic number“5F3759DF”究竟是干什么的？

为了彻底的弄明白这些问题，比如说那个magic number是什么，怎么来的，有什么用？还有为什么这个算法的计算速度如此之快，精度如此之高，等等。我查阅了一些相关资料和视频。下面是我自己的总结，作为我对这段代码的学习备忘录。

By the way：那有人要说了，这都2024年了，为什么还要研究这个算法。根据我自己的学习体会，我觉得主要就是兴趣。再其次就是这段代码从最初被创作到现已经过去这么多年了，被研究已经比较充分了。因此，在网上能够查到的相应资料很多，更有利于个人学习，再加上现在有了chatGPT的加持：）。当年很多不太清楚的地方，在后面这些年中也因为大家的钻研变得慢慢清晰了，尤其是这段代码的原作者(Greg Walsh)究竟是谁的争论，也是直到2007年随着大家的激烈讨论才逐渐浮出水面的！下面的这个名为Beyond3D的外国论坛记录了当年的讨论，以及这段代码的原作者出现的过程。

是John Carmack或Michael Abrash吗?

在讨论之前提到的《Doom3》引擎中 NV40 的渲染路径时，这段代码被提到并归因于 John Carmack；他是显而易见的选择，因为它出现在他引擎的源代码中。Michael Abrash 也被认为可能是作者。Michael 以 x86 汇编优化专家的身份脱颖而出，他撰写了传奇的《Zen of Assembly Language》和《Zen of Graphics Programming》等书，并且在《Quake》开发期间与 John 一起工作，负责优化当时的 CPU 上的 Quake 软件渲染器。

为了确定这段代码的究竟是不是他写的，有人写信问他，于是就有了下面的这段问答。(这一讨论发生在2004年)

上述回信的翻译：

你好，John，

Beyond3D.com 论坛上正在讨论这个问题，谁是以下代码的作者：

这可以归因于你吗？分析显示它的原理极其巧妙，据说来自 Quake 3 的源代码。大多数人说这是你的作品，也有人说是 Michael Abrash 的。你知道是谁写的吗？可以讲一下它的历史吗？

不是我，我也不认为是 Michael。也许是 Terje Mathison？

John Carmack

是Terje Mathisen吗？

尽管John Carmack在他的回信中否认自己是这段代码的作者，同时他也不认为这段代码出自Michael Abrash之手。但他提到了另一个可能的线索/可能的作者Terje Mathison，但随后针对这一问题的讨论便没了下文。

       直到最近(时间来了2005年的八月)，John 在今年的 QuakeCon 演讲中提到完全开源 Quake 3 v1.32 源代码，包括 3D 渲染器，得到了在场观众的热烈欢呼。随着《Doom3》的发布并赢得了好评，黑客们开始关注 id 公司，问 Quake 3 这个相对较古老的渲染器什么时候可以供人们学习和使用。

        QuakeCon 后的一周，Quake 3 源代码被正式公开，Slashdot 也报道了这个显而易见的新闻，并再次提到了快速近似倒数平方根代码的作者问题。很快我想到了 Terje，但当时却忘了问他！

        Terje Mathisen是 x86 汇编语言优化领域的大师之一。回到 3D 图形的早期软件优化时代，像 Michael、Terje（当然还有 John）这样的人会花费大量时间对关键的性能代码进行手写汇编优化。这个投资在《Doom》和《Quake》的时代取得了显著的回报。你可以在 comp.lang.asm.x86 中看到 Terje 分享的建议、优化、轶事和代码片段，这个人不是一般的厉害。那么Terje是真正的作者吗？

上述回信的翻译：

Ryszard 写道：

嘿 Terje，

自从 id 公司公开了 Quake 3 Arena 源代码后，这个问题再次出现。

你是写出这个快速倒数平方根实现的人吗？如果是，你能讲讲它的来历和你如何想出这个算法吗？一大批黑客和极客都想知道这个问题。由于 John 说不是他，可能也不是 Michael，那是你吗？

        你好，Ryszard，还有你好 John，好几年没见了。:-(

        谢谢你把我列为可能的作者，当我第一次看到这个问题时，我确实以为这是我写的代码。:-)

        五年前，我确实为了帮助一位瑞典朋友解决流体化学计算问题，写过一个非常快速且准确的invsqrt()函数。这使得他的模拟运行时间从原来的一周缩短到了一半，且能保证8-10位有效的数字。

        然而，文中的代码却不是我写的，我猜它的准确度应该只有百分之一以内？瑞典朋友需要至少 48 位有效位的精度，为此我采用了更直接的查表法和牛顿拉夫逊迭代法。由于水分子包含三个原子，我可以同时计算三个invsqrt()值，从而避免了 FP 流水线中几乎所有的气泡。

        不过，我感觉Q3A的代码风格更类似于MIT的HAKMEM文件中的一些代码。:-)

Terje

是Gary Tarolli吗？

这样看来Terje也不是真正的作者，尽管如此，我们也在他的回信中看到了他的实力。他基于查表法和NR迭代法所写的invsqrt()汇编版的精度达到了惊人的48位。而且他回忆说，这些代码的风格类似于 MIT 的 HAKMEM 文档中的代码。

随着John Carmack和Terje Mathisen都否定了他们是平方根倒数快速算法这段代码的作者，且Michael Abrash也间接被否定了，看来要想找到原作者似乎需要另谋出路。通过 Google搜索，我们发现NVIDIA也可能与之有关，并找到了一篇提到了NVIDIA的某位员工---Slashdot的文章。经过一番询问，我们得知 3dfx 的创始人之一Gary Tarolli是最有可能知道这段代码来源的人。

因此，我又向 Gary 发了一封邮件，询问他是否是本代码的作者？

上述回信的翻译：



一段过去的记忆！(时间来了2005年的九月)

我确实认得下面的代码，但这不是我的功劳。

我记得大约十年前碰巧看到过它，我当时还重新推导过它。这是牛顿-拉夫逊迭代法，加上个非常巧妙的近似初值。

        我记得我当时尝试过除了0x5f3759df的其他值。我可能是在IRIS Indigo做了与之相关的工作，或者是为Kubota咨询时而做的，我现在不太确定了。

        鉴于它所执行的数学运算量、准确度以及不需要查表的特性，这真是段非常出色的代码。

        我特别喜欢“i = 0x5f3759df - (i >> 1);”这一整数运算。实际上它是在用整数做浮点运算——我当时花了很长时间才弄明白它究竟是如何计算的，但现在我已经记不清细节了。

啊，那些日子——快的整数运算，慢的浮点运算……

这样看来，我确实在很多年前在我的键盘上敲过这段代码，我还稍微调整过那个十六进制常数，但除了我经常使用它，并且可能因此促成了它的广为流行和大量使用。我对这段代码再没有什么其他的贡献了。

附言：抱歉花了这么久才回复。

从Gary的回复中我们得知，他确实认得这段代码，但他并不认为这是他的原创。他回忆起大约 10 年前曾经遇到过这段代码，并重新推导过它。他推测这段代码是使用牛顿迭代法和一个非常巧妙的初始近似值。虽然他曾经在调试这段代码的十六进制常量，但他并不记得详细的内容了。但他确实是为该代码的广为流行做出贡献的人之一。

虽然，我们到目前为止都还没有找到代码的原作者，但经过 John Carmack、Michael Abrash、Terje Mathison 和 Gary Tarolli 的多方反馈，我们可能已经尽可能接近了真相。这段代码是 3D 图形编程领域的一个经典案例，它展示了软件在 3D 性能分析中的重要性。

最后，我想感谢 John Carmack、Terje Mathison、Gary Tarolli以及其他为这段代码的传播和优化做出贡献的人。

至此，关于代码原作者的讨论也告以段落了，且上面的讨论也被一个名叫Slashdot的报导了。巧的是，正是由于他的报导，这段代码的真正作者看到了报导，并主动联系了相关人员，这才有了下面的故事。

真正的作者找到了 --- 是Greg Walsh！

        我们本以为这就是故事的结局，于是发表了文章并宣传出去。Slashdot报道了此事，给文章带来了相当大的曝光度。正是这种曝光度让真正的作者站了出来并承认了这段代码。原来，这段代码并非出自Gary之手，而是另一位名字以G开头的硅谷老将——Greg Walsh。

        Greg自上世纪70年代以来一直是计算机行业的老兵，他的第一批大项目包括在“互联网”诞生前研究互联网和分布式计算技术，并在斯坦福大学工作期间，协助Xerox PARC开发了第一个WYSIWYG文字处理器。之后，Greg参与创立了Ardent Computer，关于该公司的更多信息，你可以在它的维基百科词条中找到。

        Ardent专门开发使用定制矢量处理器和Intel i860 RISC矢量处理器的并行图形小型计算机（Intel i860曾用于老NeXT机器中的图形子系统）。Ardent后来在Kubota（当时Ardent的资助方）的强迫下与其最大竞争对手Stellar合并，形成了Stardent。



        Ardent的Titan图形小型计算机当时在实现其性能目标时遇到了困难，而Greg最初是为了加速无法利用矢量硬件的软件而设计了快速1/sqrt(x)函数。当时，Greg与MathWorks的创始人兼MatLab作者Cleve Moler一起在Ardent工作。正是与Cleve的合作让Greg萌生了编写该函数的想法。Cleve当时正在研究使用牛顿-拉夫森迭代法进行近似计算，而Greg也在同一时期为Titan编写了类似的1/cuberoot(x)函数，不过方法更加复杂。

        与此同时，Gary Tarolli也在为Kubota提供咨询服务，因此我们终于明白了代码如何进入公开源码。因为Gary Tarolli也在使用并调整该代码以适应自己的需求，而他在3dfx的工作至少可以通过Brian Hook与id Software建立联系。Gary和Greg后来在Accel Graphics公司再次合作（Accel最终被Evans and Sutherland收购，另一个涉及3D领域的公司，关于它的故事你也可以写一本书）。如果有谁有多余的AccelSTAR硬件，欢迎联系我！

        因此，最终我们认定Greg是这些年来参与编写或修改这段代码的人之一。我们可以理直气壮地称这段代码灵感来自于Cleve，Greg为作者，现在大概可以为这段故事画上句号了。

        Greg后来离开3D行业，创办了一家商业软件公司，成功上市并在1999年引起轰动；Kubota最终为DEC工作站生产了插入式3D板（Denali现在大多是罕见的收藏品），然后彻底退出了3D行业；而Cleve仍然在MathWorks担任首席科学家兼创始人。

        非常感谢Greg与我们联系并提供了这些信息，也感谢Slashdot（以及提交该文章的Geo），因为没有那次大曝光，Greg可能永远不会看到原文。

上述故事中的文字，都是我用chatGPT从外国论坛中翻译的，我只是对译文做了一丁点修改和排版。

正文

他(Greg Walsh)的这套代码，我学习下来。其核心思想就是牛顿拉夫逊法，而众所周知，牛顿拉夫逊法的核心是通过多次迭代提升计算精度/缩小误差，最终得到我们想要的答案。而Greg Walsh这段代码的精妙之处在于出色的选择牛顿拉夫逊法的初值(对牛顿拉夫逊法而言，初值选择的越是接近真实答案就越理想，迭代后的准确率就越高)，使得这段代码只需迭代一次就能接近使用其他初值多次迭代后的精度。

为了达到这一目的，Greg Walsh在计算/选择初值的时候使用了代码中的那个神秘的magic number“5F3759DF”。而带有magic number的那句计算初值的代码，又巧妙地利用了浮点数基于IEEE 754标准在计算中的默认存储方式/形式。

现在，我把上面两段我所要表达的重点重新梳理一下，并换一种方式来陈述。Greg Walsh的这套计算平方根倒数的快速算法，实际上还包含了一个“内嵌的一个平方根倒数的快速算法”。“内嵌的平方根倒数快速算法”专门用于计算 $1/\sqrt{x}$ 的初值(一个比较精确的近似值)。由于下面几句代码实现：

有了前面的初值后，后面只用了一次牛顿拉夫逊法就完成了全部的计算，即，下面的这一句代码：

下面是对代码原理的详细剖析：

1，问题的转化

给定一个任意数x，要计算 $1/\sqrt{x}$ ，假设 $1/\sqrt{x}$ 的值为a，则：

$1/\sqrt{x}=a$ ，（式1）

$\Rightarrow 1/a^{2}-x=0$

如果定义一个自变量为a的函数f(a):

$f(a)= 1/a^{2}-x$ ，（式2）

则，令函数f(a)等于0的a就是我们要求的 $1/\sqrt{x}$ ：

$f(a)= 0\Rightarrow 1/a^{2}-x=0$

如此一来，我们最开始的问题就转化为要找到能够令函数f(a)等于0的a，即: $a=?\Rightarrow f(a)=0$ 。如果把函数f(a)画出来，则函数f(a)与a轴的交点(a，0)就是我们要找的根。

函数的根指的是使函数的值为零的自变量，在函数的图像中，一个函数的根即这个函数与横坐标的交点。

以x=1为例：

$1/\sqrt{1}=1$ 得到a=1，画出函数 $f(a)= 1/a^{2}-1$ 的图像可见f(a)与横坐标轴a的交点正好是(a=1,y=0)。

python code:

import matplotlib.pyplot as plt
import numpy as np

# make data
x=1
a = np.linspace(0.01,8, 1000)
y = 1/(a**2)-x

# plot
fig, ax = plt.subplots()
ax.plot(a, y, linewidth=2.0,label='f(a)')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)

2，用牛顿拉夫逊法(NR-iteration)求x=1时的 $1/\sqrt{1}$

在上面的例子中，由于x取得比较特殊，我们直接求出了 $1/\sqrt{1}$ ，即，令函数f(a)为0的根。现在我们依然以x=1为例，依然是求 $1/\sqrt{1}$ ，依然是通过求解f(a)的根。只不过是用NR-iteration去求f(a)=0的近似解。

令x=1，代入(式2)有：

$f(a)= 1/a^{2}-x=1/a^{2}-1$ ,（式3）

对(式3)求导，得到:

$f'(a)= -2/a^{3}$ ,（式4）

初次迭代

step1: 随机给定一个初始值a0，并求出相应的f(a0)

用牛顿拉夫逊法需要逐级逼近，虽然我们已经知道最终的答案是a=1，但我这里假设我并不知道。所以，我随便取了一个真实答案左边的值a0=0.5(注意：这里不能取1左边的值，否则牛顿拉夫逊法会失效)，代入(式3)有：

$a_{0}=0.5\Rightarrow f(a_{0})=1/{0.5}^{2}-1=3.0$

得到点(0.5，3.0)，这里的3.0表示实际答案与真实答案之间的误差err。

step2: 把a0=0.5代入(式4)计算出该点处的斜率，并利用点斜式求出函数f(a)在点(0.5，3.0)处的切线方程(称切线为line0)

$s= -2/a_{0}^{3}= -2/0.5^{3}=-16$

line0: $y-3.0=s*(a-0.5)\Rightarrow y-3.0=-16*(a-0.5)$