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

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

news 2025/7/10 21:03:44

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

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

相应的输出为：