当前位置：首页 > news >正文

数值计算 --- 平方根倒数快速算法(0x5f3759df,这是什么鬼!!!)

news 2024/11/15 15:11:22

平方根倒数快速算法 --- 向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可能永远不会看到原文。

正文

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

为了达到这一目的，Greg Walsh在计算/选择初值的时候用了代码中的那个神秘的magic number“5F3759DF”，这使得牛顿拉夫逊法在这段代码中只用迭代一次。而带有magic number的那句计算初值的代码，又巧妙地利用了浮点数x在计算中默认的存储方式/形式。

现在，我把上面的所说的两段重新梳理一下，并换一种方式来陈述也许更有助于理解。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$

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

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

如此一来，我们最开始的问题就转化为要找到能够令函数f(a)等于0的a。如果把函数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，用牛顿拉夫逊法求x=1时的 $1/\sqrt{1}$

函数f(a):

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

函数f(a)的导数:

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

初次迭代

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

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

$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代入式3计算出该点处的斜率，并利用点斜式求出点(0.5,3.0)处切线line0的方程

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

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

step3: 令y=0，找到切线line0与a轴的交点得到新的a，并称其为a1

$0-3.0=-16*(a-0.5)\Rightarrow a=0.6875$

python code:

#guess一个初值x0,并求出相应的f(a0)
a0=0.5
fa0=1/(a0**2)-x
print(f"a={a0}, Err:f(a)={fa0}")#利用点斜式求出点(a0,f(a0))处的切线line0
#y-f(a0)=s*(a-a0)
aa=np.linspace(-1, 8, 1000)
s=-2*(a0**(-3))
print(f"slope={s}")
y0=s*(aa-a0)+fa0#令y=0,找到切线line0与a轴的交点得到新的a1
a1=a0-fa0/s
print(f"a_new={a1}")# plot
fig, ax = plt.subplots()
ax.plot(a, fa, linewidth=2.0,label='f(a)')
ax.plot(aa, y0,'--',label='line0')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)
plt.show()

相应的log:

第二次迭代

基于新的a，即上面的a1，在函数图像上找到新的点(a1，f(a1))。画该点处的切线line1，得到line1与a轴的交点又会得到一个新的a，最终完成第二次迭代。具体的计算过程我这里不再赘述，我的python代码注释中有详细的说明。

python code:

#把a1代入函数f(a)求出相应的f(a1)
fa1=1/(a1**2)-x
print(f"a={a1}, Err:f(a)={fa1}")#利用点斜式求出点(a1,f(a1))处的切线line1
#y-f(a1)=s*(a-a1)
s=-2*(a1**(-3))
y1=s*(aa-a1)+fa1
print(f"slope={s}")#令y=0,找到切线line1与a轴的交点得到新的a2
a2=a1-fa1/s
print(f"a_new={a2}")# plot
fig, ax = plt.subplots()
ax.plot(a, fa, linewidth=2.0,label='f(a)')
ax.plot(aa, y0,'--',label='line0')
ax.plot(aa, y1,'--',label='line1')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)
plt.show()

相应的log:

第三次迭代

后续的迭代无非就是个重复上面的过程，求点，画切线和求交点更新a，直到误差err的值足够小，到那个时候，我们就声称已经求得最终的a了，即， $1/\sqrt{1}$ 的值(准备的说是达到一定精度的近似值)。

python code:

#把a2代入函数f(a)求出相应的f(a2)
fa2=1/(a2**2)-x
print(f"a={a2}, Err:f(a)={fa2}")#利用点斜式求出点(a2,f(a2))处的切线line2
#y-f(a2)=s*(a-a2)
s=-2*(a2**(-3))
y2=s*(aa-a2)+fa2
print(f"slope={s}")#令y=0,找到切线line2与a轴的交点得到新的a3
a3=a2-fa2/s
print(f"a_new={a3}")# plot
fig, ax = plt.subplots()
ax.plot(a, fa, linewidth=2.0,label='f(a)')
ax.plot(aa, y0,'--',label='line0')
ax.plot(aa, y1,'--',label='line1')
ax.plot(aa, y2,'--',label='line2')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)
plt.show()

相应的log:

第四次迭代

python code:

#把a3代入函数f(a)求出相应的f(a3)
fa3=1/(a3**2)-x
print(f"a={a3}, Err:f(a)={fa3}")#利用点斜式求出点(a3,f(a3))处的切线line3
#y-f(a3)=s*(a-a3)
s=-2*(a3**(-3))
y3=s*(aa-a3)+fa3
print(f"slope={s}")#令y=0,找到切线line3与a轴的交点得到新的a4
a4=a3-fa3/s
print(f"a_new={a4}")# plot
fig, ax = plt.subplots()
ax.plot(a, fa, linewidth=2.0,label='f(a)')
ax.plot(aa, y0,'--',label='line0')
ax.plot(aa, y1,'--',label='line1')
ax.plot(aa, y2,'--',label='line2')
ax.plot(aa, y3,'--',label='line3')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)
plt.show()

第五次迭代

python code:

#把a4代入函数f(a)求出相应的f(a4)
fa4=1/(a4**2)-x
print(f"a={a4}, Err:f(a)={fa4}")#利用点斜式求出点(a4,f(a4))处的切线line4
#y-f(a4)=s*(a-a4)
s=-2*(a4**(-3))
y4=s*(aa-a4)+fa4
print(f"slope={s}")#令y=0,找到切线line4与a轴的交点得到新的a5
a5=a4-fa4/s
print(f"a_new={a5}")# plot
fig, ax = plt.subplots()
ax.plot(a, fa, linewidth=2.0,label='f(a)')
ax.plot(aa, y0,'--',label='line0')
ax.plot(aa, y1,'--',label='line1')
ax.plot(aa, y2,'--',label='line2')
ax.plot(aa, y3,'--',label='line3')
ax.plot(aa, y4,'--',label='line4')
ax.set(xlim=(-1, 8),ylim=(-2, 6))
ax.legend()
ax.grid(True)
plt.show()

下面这张表格记录了上面每一步的结果，可以看到牛顿拉夫逊法的收敛速度还是比较快的。经过5次迭代基本上已经非常接近真实值1了。

牛顿拉夫逊迭代法的更新过程
iteration num	a	a_new	Err
1	0.5	0.687	3
2	0.687	0.868	1.12
3	0.868	0.975	0.32
4	0.975	0.9990923	0.0513
5	0.9990923	0.9999988	0.0018

在完全熟悉了牛顿拉夫逊法以后，上面的5次迭代过程都可以通过下面的计算通式来表示，从而免去了前面繁琐的计算过程。其中， $x_{n}$ 表示前一次的计算结果， $x_{n+1}$ 表示本次迭代更新后的x。

就本例而言，则上式变成:

$a_{n+1}=a_{n}-f(a_{n})/f'(a_{n})$ ，(迭代方程)

其中：

$f(a_{n})= 1/{a_{n}}^{2}-x$

$f'(a_{n})= -2*1/{a_{n}}^{3}$

代入迭代方程，得到：

$a_{n+1}=a_{n}-f(a_{n})/f'(a_{n})$

$\Rightarrow a_{n+1}=a_{n}(1.5-x{a_{n}}^{2}/2)$ ，(式4)

现在我们来试一下这个公式，实际上他得到结果将和我们上面的计算结果是一样的。令a的初值a=0.5代入式3，然后不断的把新得到的a代入式4迭代，得到如下结果：

上面的迭代过程透露了一个信息，如果我们初值选择的足够好，我们就能很快完成迭代。这里我们以x=2为例，即求 $1/\sqrt{2}$ 。因为我们都知道根号2的值大约是1.414，所以根号2的倒数也就是约等于1/1.414约等于0.707。为了用牛顿法求得更加精确的值，我们令初值a=0.707并代入迭代公式，得到：

第一次迭代后的估计值：

用python计算器算出来的精确值：

相当于只迭代了一次就已经精确到了小数点后第7位！！！而如果初值仍然是选择0.5的话，需迭代5次才能超过上面迭代一次后的精度。

而这就是平方根倒数快速算法的两大核心：1，采用牛顿拉夫逊法计算近似值。2，尽可能精确的选择牛顿法的初值。现在我们再回到前面推导的结果式4：

$a_{n+1}=a_{n}(1.5-x{a_{n}}^{2}/2)$ ，(式4)

如果我们令 $a_{n+1},a_{n}$ 为y，再令x/2为x2，则式4变为：

$y=y*(1.5-x2*y^{2})=y*(1.5-x2*y*y)$

这正是平方根倒数快速算法的code中的最后一行：

这说明了以下几点：

1，这段代码使用了牛顿拉夫逊法(毕竟code中的代码和推导出的公式一样)

2，这段代码中所使用的初值y，必然是一个非常精确的初始值

3，这个非常精确的初始值就是批注为“what the fuck？”那一行求出来的。

4，也许是因为初始值选的好，这段代码中的第二次牛顿迭代被注释掉了。也就是说整个函数只迭代了一次就达到了一个非常高的精度。

如何选择正确的初值?以及WTF中的神秘数字究竟是怎么来的

花开两朵，各表一枝。选择正确而相对精确的初值的关键在于如何求出 $1/\sqrt{x}$ 的近似值，而求 $1/\sqrt{x}$ 的快速算法在于充分的利用浮点数x在计算机中的表示/编码方式。

因为x为浮点数(注意：我这里的x就是代码中的number)。则根据标准IEEE 754，x的二进制浮点数表示如下(准确的说应该叫normal number的表示)：

$(-1)^{S}\times 2^{E-b}\times (1+T\cdot 2^{1-p})$

又因为x不能为负(负数没法进行开根号运算)，符号位S默认为0，则浮点数x在计算机中的二进制可表示如下：

$x= (1+T\cdot 2^{1-p})\times 2^{E-b}$

对于单精度float而言，p=24，b=127，则：

$x= (1+T\cdot 2^{-23})\times 2^{E-127}$

我们对x取以2为底的对数，得到：

${log_{2}}^{x}={log_{2}}^{(1+T\cdot 2^{-23})\times 2^{E-127}}={log_{2}}^{(1+T\cdot 2^{-23})}+{E-127}$

再令：

$M=T\cdot 2^{-23}$

则上式变为：

${log_{2}}^{x}={log_{2}}^{(1+T\cdot 2^{-23})\times 2^{E-127}}$

$={log_{2}}^{(1+T\cdot 2^{-23})}+{E-127}={log_{2}}^{(1+M)}+{E-127}$

即：

${log_{2}}^{x}={log_{2}}^{(1+M)}+{E-127}$ ，(式5)

注意，T字段所保存的是trailing significand，即，放大一定精度后的有效数字的尾数/有效数字的小数部分(默认隐含了首位1)。计算机在保存T时把小数点右移了23位，即，乘以 $2^{23}$ 。因此，在读取T时才有了上面的 $T\cdot 2^{-23}$ 。这就是说上面的M实际上是“1.xxxxx...”中的“0.xxxxx...”部分，是一个介于0~1之间的数。

为了更好的理解M，这里插播一个例子，1/3是如何被保存成二进制浮点数的？

        计算机使用二进制浮点数表示小数时，采用的是 IEEE 754 浮点数标准。由于1/3是一个无限循环小数，在二进制中它也不能被精确表示，所以计算机只能以有限的精度近似存储它。

1. 十进制转二进制

在十进制下，1/3=0.33333...是一个循环小数。转换到二进制后为：

$\frac{1}{3}_{10}=0.333..._{10}=0.01010101..._{2}$

        也是一个二进制的循环小数，但由于计算机只能只能保存有限的位数，这个循环小数在保存时会被截断，得到一个近似值。

2. IEEE 754 浮点数表示

在 IEEE 754 单精度浮点数标准中，32位浮点数的表示结构如下：

1 位符号位：表示正数或负数
8 位指数：存储实际指数的偏移量（偏移 127）
23 位尾数（有效数字）：存储归一化的尾数，隐含首位为 1 的小数部分

对于1/3计算机会将其转换为二进制表示，然后使用以下步骤：

标准化二进制小数：将二进制小数表示成规范形式。规范形式要求小数点左侧只能有一位，且必须是1，因此：
$0.01010101..._{2}=1.01010101..._{2}\times 2^{-2}$

计算指数E：指数部分需要加上偏移量(127)。所以，计算机所保存的指数E等于上面的实际指数−2加上127。−2+127=125，再转换为二进制后为 01111101。

有效数字的尾数T：有效数字尾数的精度共 23 位，因此我们在保存小数部分时，去掉整数部分的1不保存： $1.01010101..._{2}\Rightarrow 0.01010101..._{2}$

        然后再把小数点右移23位，得到：

$0.01010101..._{2}\Rightarrow 01010101010101010101010_{2}$

4.最终存储形式：将符号位（0，正数）、指数（125 的二进制表示 011111010111110101111101）和尾数组合起来，得到：

$00111110101010101010101010101010$

这就是1/3的IEEE 754单精度浮点数表示。

由于M是一个0~1之间的小数，人们发现当M=0~1时，函数y=log2(1+M)与y=M的函数值差异很小。

Matlab code:

close all
clear allx=0.01:0.01:pi/2;
f1=log2(1+x);
f2=x;
plot(x,f1,x,f2);
grid on;
legend("y=log2(1+M)","y=x")diff=abs(f1-f2);
figure
plot(x,diff)
legend("diff")

因此，我们认为在x=0~1之间:

${log_{2}}^{(1+M)}\approx M$

基于这一近似，式5变为：

${log_{2}}^{x}={log_{2}}^{(1+M)}+E-127=M+E-127$

$=T\cdot 2^{-23}+E-127=1/2^{-23}\cdot (E\times 2^{23}+T)-127$

又因为括号中的 $E\times 2^{23}+T$ ，正好是浮点数x在计算机中的存储形式（我们这里用 $x_{B}$ 来表示），即：

$x_{B}=E\times 2^{23}+T$

这里，我们再插播一下。如果还是以上面插播信息中的1/3为例的话。我文章中的x就是1/3(十进制)，而 $x_{B}$ 就是上面那个例子中最终保存的 $00111110101010101010101010101010$ 。他们是一个数，只不过一个是实际数，一个是在计算机中存的数。

如此一来，我们利用浮点数x在计算机中默认的二进制存储方式，得到了log2(x)的表示方式：

${log_{2}}^{x}=1/2^{-23}\cdot (E\times 2^{23}+T)-127=1/2^{-23}\cdot x_{B}-127$

${log_{2}}^{x}=x_{B}/2^{23}-127$ ，（式6）

现在我们再回到计算 $1/\sqrt{x}$ 的近似值问题。根据（式1）我们知道：

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

对上式两边同时取以2为底的对数，得到：

${log_{2}}^{a}={log_{2}}^{x^{-1/2}}$

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

根据前面推导出的log2(x)的表示方式（式6）：

${log_{2}}^{a}=a_{B}/2^{23}-127$ ， $-1/2\cdot {log_{2}}^{x}=-1/2\cdot (x_{B}/2^{23}-127)$

$a_{B}/2^{23}-127=-1/2\cdot (x_{B}/2^{23}-127)$

$a_{B}=381\times 2^{22}-x_{B}/2$

其中 $381\times 2^{22}=1598029824$ 这个数，如果用十六进制来表示的话就是：

$381\times 2^{22}=5f400000$

则上式变为：

$a_{B}=5f400000-x_{B}/2$ ，（式7）

这个十六进制的数code中的那个神秘数字“5f3759df”已经比较接近了，而这个数表示成十进制是1597463007。

这里我们暂时先不讨论这两个十六进制常数的差异，先看看(式7)究竟表示什么意思：

$a_{B}=5f400000-x_{B}/2$ ，（式7）

我们知道a就是我们要求的十进制数x的平方根的倒数，而我们又知道不论十进制数a或x是多少，他在计算机中都要以二进制浮点数的方式被保存为 $a_{B}$ 和 $x_{B}$ 的形式。因此，(式7)的意思是说，对于一个已经按照IEEE 754标准被保存好的十进制浮点数x，他在计算机中换了个样子，变成了 $x_{B}$ ，但他仍然等于x。而要想求得 $x_{B}$ 的平方根的倒数，只需按照(式7)就能快速求出近似值 $a_{B}$ ，这个 $a_{B}$ 是与之对应的十进制浮点数a，保存在计算机中的样子。而要想把 $a_{B}$ 再变成a，只需按照浮点数的编码方式解析出来即可。

现在让我们再回到原代码，我们注意到评论为WTF的上下两句所做的正是我在上文中所描述的过程。所不同的是代码中的y是我文中的x，代码中的i是我文中的 $x_{B}$ ，代码中的经过神秘数字“5f3759df”计算后的新i是我文中的 $a_{B}$ ,而把新i重新解码后的浮点数y是我文中的a：

现在，我们有了能够快速求解出较为精确的 $1/\sqrt{x}$ 的公式(式7)，再加上之前根据牛顿拉夫逊法求得的(式4) $a_{n+1}=a_{n}(1.5-x{a_{n}}^{2}/2)$ 。至此，我们基本上复现了平方根倒数快速算法的全部过程，且和原始code一致(除了magic number之外)。

我们来试试我们现有的快速算法，看看他的效果究竟怎么样，还是以x=1为例，求 $1/\sqrt{1}$ 。

C code:

# include <stdio.h>
# include <math.h>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 hackingi = 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 removedreturn y;
}float myQ_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 hackingi = 0x5f400000 - (i >> 1);y = *(float*)&i;y = y * (threehalfs - (x2 * y * y));   // 1st iteration// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removedreturn y;
}int main() {float x = 4.0f;float y = 0,yy=0;y=Q_rsqrt(x);yy = myQ_rsqrt(x);printf("input x=%f\n", x);printf("ideal result=%f\n", 1/sqrt(x));printf("calc with 5f3759df=%f\n", y);printf("calc with 5f400000=%f\n", yy);return 0;
}

相应的输出为：