可以免费做试卷题目的网站,h5网站模板下载,seo单页面优化,室内设计公司排名十强平方根倒数快速算法 --- 向Greg Walsh致敬#xff01; 1#xff0c;牛顿拉夫逊 已知x#xff0c;要计算#xff0c;假设的值为a#xff0c;则#xff1a;
#xff0c;#xff08;式1#xff09; 如果定义一个自变量为a的函数f(a): 则#xff0c;令函数f(a)等于0的a就…平方根倒数快速算法 --- 向Greg Walsh致敬 1牛顿拉夫逊 已知x要计算假设的值为a则
式1 如果定义一个自变量为a的函数f(a): 则令函数f(a)等于0的a就是我们要求的 如此一来我们最开始的问题就转化为要找到能够令函数f(a)等于0的a。如果把函数f(a)画出来则函数f(a)与a轴的交点(a0)就是我们要找的根(一个函数的根是指使函数的值为零的自变量)。 以x1为例即a1画出函数的图像可见f(a)与横坐标轴a的交点正好是(a1,y0)。 python code: import matplotlib.pyplot as plt
import numpy as np# make data
x1
a np.linspace(0.01,8, 1000)
y 1/(a**2)-x# plot
fig, ax plt.subplots()
ax.plot(a, y, linewidth2.0,labelf(a))
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True) 2用牛顿拉夫逊法求x1时的 函数f(a):
,式2
函数f(a)的导数:
,式3 初次迭代 step1: 随机给定一个初始值a0并求出相应的f(a0) 用牛顿拉夫逊法需要逐级逼近虽然我们已经知道最终的答案是a1但我这里假设我并不知道。所以我随便取了一个真实答案左边的值a00.5(注意这里不能取1左边的值否则牛顿拉夫逊法会失效)代入式2有 得到点(0.5,3.0)这里的3.0表示实际答案与真实答案之间的误差err。 step2: 把a00.5代入式3计算出该点处的斜率并利用点斜式求出点(0.5,3.0)处切线line0的方程 line0: step3: 令y0找到切线line0与a轴的交点得到新的a并称其为a1 python code: #guess一个初值x0,并求出相应的f(a0)
a00.5
fa01/(a0**2)-x
print(fa{a0}, Err:f(a){fa0})#利用点斜式求出点(a0,f(a0))处的切线line0
#y-f(a0)s*(a-a0)
aanp.linspace(-1, 8, 1000)
s-2*(a0**(-3))
print(fslope{s})
y0s*(aa-a0)fa0#令y0,找到切线line0与a轴的交点得到新的a1
a1a0-fa0/s
print(fa_new{a1})# plot
fig, ax plt.subplots()
ax.plot(a, fa, linewidth2.0,labelf(a))
ax.plot(aa, y0,--,labelline0)
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True)
plt.show() 相应的log: 第二次迭代 基于新的a即上面的a1在函数图像上找到新的点(a1f(a1))。画该点处的切线line1得到line1与a轴的交点又会得到一个新的a最终完成第二次迭代。具体的计算过程我这里不再赘述我的python代码注释中有详细的说明。 python code: #把a1代入函数f(a)求出相应的f(a1)
fa11/(a1**2)-x
print(fa{a1}, Err:f(a){fa1})#利用点斜式求出点(a1,f(a1))处的切线line1
#y-f(a1)s*(a-a1)
s-2*(a1**(-3))
y1s*(aa-a1)fa1
print(fslope{s})#令y0,找到切线line1与a轴的交点得到新的a2
a2a1-fa1/s
print(fa_new{a2})# plot
fig, ax plt.subplots()
ax.plot(a, fa, linewidth2.0,labelf(a))
ax.plot(aa, y0,--,labelline0)
ax.plot(aa, y1,--,labelline1)
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True)
plt.show()
相应的log: 第三次迭代 后续的迭代无非就是个重复上面的过程求点画切线和求交点更新a直到误差err的值足够小到那个时候我们就声称已经求得最终的a了即的值(准备的说是达到一定精度的近似值)。 python code: #把a2代入函数f(a)求出相应的f(a2)
fa21/(a2**2)-x
print(fa{a2}, Err:f(a){fa2})#利用点斜式求出点(a2,f(a2))处的切线line2
#y-f(a2)s*(a-a2)
s-2*(a2**(-3))
y2s*(aa-a2)fa2
print(fslope{s})#令y0,找到切线line2与a轴的交点得到新的a3
a3a2-fa2/s
print(fa_new{a3})# plot
fig, ax plt.subplots()
ax.plot(a, fa, linewidth2.0,labelf(a))
ax.plot(aa, y0,--,labelline0)
ax.plot(aa, y1,--,labelline1)
ax.plot(aa, y2,--,labelline2)
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True)
plt.show() 相应的log: 第四次迭代 python code: #把a3代入函数f(a)求出相应的f(a3)
fa31/(a3**2)-x
print(fa{a3}, Err:f(a){fa3})#利用点斜式求出点(a3,f(a3))处的切线line3
#y-f(a3)s*(a-a3)
s-2*(a3**(-3))
y3s*(aa-a3)fa3
print(fslope{s})#令y0,找到切线line3与a轴的交点得到新的a4
a4a3-fa3/s
print(fa_new{a4})# plot
fig, ax plt.subplots()
ax.plot(a, fa, linewidth2.0,labelf(a))
ax.plot(aa, y0,--,labelline0)
ax.plot(aa, y1,--,labelline1)
ax.plot(aa, y2,--,labelline2)
ax.plot(aa, y3,--,labelline3)
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True)
plt.show() 第五次迭代 python code: #把a4代入函数f(a)求出相应的f(a4)
fa41/(a4**2)-x
print(fa{a4}, Err:f(a){fa4})#利用点斜式求出点(a4,f(a4))处的切线line4
#y-f(a4)s*(a-a4)
s-2*(a4**(-3))
y4s*(aa-a4)fa4
print(fslope{s})#令y0,找到切线line4与a轴的交点得到新的a5
a5a4-fa4/s
print(fa_new{a5})# plot
fig, ax plt.subplots()
ax.plot(a, fa, linewidth2.0,labelf(a))
ax.plot(aa, y0,--,labelline0)
ax.plot(aa, y1,--,labelline1)
ax.plot(aa, y2,--,labelline2)
ax.plot(aa, y3,--,labelline3)
ax.plot(aa, y4,--,labelline4)
ax.set(xlim(-1, 8),ylim(-2, 6))
ax.legend()
ax.grid(True)
plt.show() 下面这张表格记录了上面每一步的结果可以看到牛顿拉夫逊法的收敛速度还是比较快的。经过5次迭代基本上已经非常接近真实值1了。
牛顿拉夫逊迭代法的更新过程 iteration numaa_newErr10.50.687320.6870.8681.1230.8680.9750.3240.9750.99909230.051350.99909230.99999880.0018 在完全熟悉了牛顿拉夫逊法以后上面的5次迭代过程都可以通过下面的计算通式来表示从而免去了前面繁琐的计算过程。其中表示前一次的计算结果表示本次迭代更新后的x。 就本例而言则上式变成: (迭代方程) 其中 代入迭代方程得到 (式4) 现在我们来试一下这个公式实际上他得到结果将和我们上面的计算结果是一样的。令a的初值a0.5代入式3然后不断的把新得到的a代入式4迭代得到如下结果 上面的迭代过程透露了一个信息如果我们初值选择的足够好我们就能很快完成迭代。这里我们以x2为例即求。因为我们都知道根号2的值大约是1.414所以根号2的倒数也就是约等于1/1.414约等于0.707。为了用牛顿法求得更加精确的值我们令初值a0.707并代入迭代公式得到 第一次迭代后的估计值 用python计算器算出来的精确值 相当于只迭代了一次就已经精确到了小数点后第7位而如果初值仍然是选择0.5的话需迭代5次才能超过上面迭代一次后的精度。 而这就是平方根倒数快速算法的两大核心1采用牛顿拉夫逊法计算近似值。2尽可能精确的选择牛顿法的初值。现在我们再回到前面推导的结果式4
(式4)
如果我们令为y再令x/2为x2则式4变为 这正是平方根倒数快速算法的code中的最后一行 这说明了以下几点 1这段代码使用了牛顿拉夫逊法(毕竟code中的代码和推导出的公式一样) 2这段代码中所使用的初值y必然是一个非常精确的初始值 3这个非常精确的初始值就是批注为“what the fuck”那一行求出来的。 4也许是因为初始值选的好这段代码中的第二次牛顿迭代被注释掉了。也就是说整个函数只迭代了一次就达到了一个非常高的精度。 如何选择正确的初值?以及WTF中的神秘数字究竟是怎么来的 花开两朵各表一枝。选择正确而相对精确的初值的关键在于如何求出的近似值而求的快速算法在于充分的利用浮点数x在计算机中的表示/编码方式。 因为x为浮点数(注意我这里的x就是代码中的number)。则根据标准IEEE 754x的二进制浮点数表示如下(准确的说应该叫normal number的表示) 又因为x不能为负(负数没法进行开根号运算)符号位S默认为0则浮点数x在计算机中的二进制可表示如下 对于单精度float而言p24b127则 我们对x取以2为底的对数得到 再令 则上式变为 即
(式5) 注意T字段所保存的是trailing significand即放大一定精度后的有效数字的尾数/有效数字的小数部分(默认隐含了首位1)。计算机在保存T时把小数点右移了23位即乘以。因此在读取T时才有了上面的。这就是说上面的M实际上是“1.xxxxx...”中的“0.xxxxx...”部分是一个介于0~1之间的数。 为了更好的理解M这里插播一个例子1/3是如何被保存成二进制浮点数的 计算机使用二进制浮点数表示小数时采用的是 IEEE 754 浮点数标准。由于1/3是一个无限循环小数在二进制中它也不能被精确表示所以计算机只能以有限的精度近似存储它。 1. 十进制转二进制 在十进制下1/30.33333...是一个循环小数。转换到二进制后为 也是一个二进制的循环小数但由于计算机只能只能保存有限的位数这个循环小数在保存时会被截断得到一个近似值。 2. IEEE 754 浮点数表示 在 IEEE 754 单精度浮点数标准中32位浮点数的表示结构如下 1 位符号位表示正数或负数8 位指数存储实际指数的偏移量偏移 12723 位尾数有效数字存储归一化的尾数隐含首位为 1 的小数部分 对于1/3计算机会将其转换为二进制表示然后使用以下步骤 标准化二进制小数将二进制小数表示成规范形式。规范形式要求小数点左侧只能有一位且必须是1因此 计算指数E指数部分需要加上偏移量(127)。所以计算机所保存的指数E等于上面的实际指数−2加上127。−2127125再转换为二进制后为 01111101。 有效数字的尾数T有效数字尾数的精度共 23 位因此我们在保存小数部分时去掉整数部分的1不保存 然后再把小数点右移23位得到 4.最终存储形式将符号位0正数、指数125 的二进制表示 011111010111110101111101和尾数组合起来得到 这就是1/3的IEEE 754单精度浮点数表示。 由于M是一个0~1之间的小数人们发现当M0~1时函数ylog2(1M)与yM的函数值差异很小。 Matlab code: close all
clear allx0.01:0.01:pi/2;
f1log2(1x);
f2x;
plot(x,f1,x,f2);
grid on;
legend(ylog2(1M),yx)diffabs(f1-f2);
figure
plot(x,diff)
legend(diff) 因此我们认为在x0~1之间: 基于这一近似式5变为 又因为括号中的正好是浮点数x在计算机中的存储形式我们这里用来表示即 这里我们再插播一下。如果还是以上面插播信息中的1/3为例的话。我文章中的x就是1/3(十进制)而就是上面那个例子中最终保存的。他们是一个数只不过一个是实际数一个是在计算机中存的数。 如此一来我们利用浮点数x在计算机中默认的二进制存储方式得到了log2(x)的表示方式 式6 现在我们再回到计算的近似值问题。根据式1我们知道 对上式两边同时取以2为底的对数得到 根据前面推导出的log2(x)的表示方式式6 其中这个数如果用十六进制来表示的话就是 则上式变为 式7 这个十六进制的数code中的那个神秘数字“5f3759df”已经比较接近了而这个数表示成十进制是1597463007。 这里我们暂时先不讨论这两个十六进制常数的差异先看看(式7)究竟表示什么意思
式7 我们知道a就是我们要求的十进制数x的平方根的倒数而我们又知道不论十进制数a或x是多少他在计算机中都要以二进制浮点数的方式被保存为和的形式。因此(式7)的意思是说对于一个已经按照IEEE 754标准被保存好的十进制浮点数x他在计算机中换了个样子变成了但他仍然等于x。而要想求得的平方根的倒数只需按照(式7)就能快速求出近似值这个是与之对应的十进制浮点数a保存在计算机中的样子。而要想把再变成a只需按照浮点数的编码方式解析出来即可。 现在让我们再回到原代码我们注意到评论为WTF的上下两句所做的正是我在上文中所描述的过程。所不同的是代码中的y是我文中的x代码中的i是我文中的代码中的经过神秘数字“5f3759df”计算后的新i是我文中的,而把新i重新解码后的浮点数y是我文中的a 现在我们有了能够快速求解出较为精确的的公式(式7)再加上之前根据牛顿拉夫逊法求得的(式4)。至此我们基本上复现了平方根倒数快速算法的全部过程且和原始code一致(除了magic number之外)。 我们来试试我们现有的快速算法看看他的效果究竟怎么样还是以x1为例求。 C code: # include stdio.h
# include math.hfloat 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,yy0;yQ_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;
} 相应的输出为 就本例而言二者的计算结果都非常接近准确值1但5f400000的精度要更高5f3759df的误差约为0.002。 如果以x4为例准确值为0.5再看看二者的表现 结果还是5f400000的准确性更高5f3759df的误差约为0.0009。但上面的两个例子都是平方根的结果正好是整数的情况例如和。但如果碰到平方根为无理数的情况呢我们分别试试x2和x3的情况。 有趣的是在这两个例子中基于magic number的计算结果要比5f400000的精度高。对于x2而言5f400000的误差约为0.0045f3759df的误差约为0.0002。 对于x3而言5f400000的误差约为0.0065f3759df的误差约为0.0005。 这样看来不论是采取哪种常数去估算初值基本上都已经能够得到较为准确的结果。毕竟这一初值还要用牛顿拉夫逊法再迭代一次才是最终的结果。 二者之间在十进制上的差为
1598029824(5f400000) -1597463007(5f3759df)566817 全文完
--- 作者松下J27 参考文献(鸣谢)
1https://en.wikipedia.org/wiki/Newton%27s_method#Examples
2什么代码让程序员之神感叹“卧槽”改变游戏行业的平方根倒数算法_哔哩哔哩_bilibili
3[算法] 平方根倒数速算法中的魔数0x5f3759df的来源 | GeT Left
4https://i.hsfzxjy.site/uncover-the-secret-of-fast-inverse-square-root-algorithm/
5https://www.youtube.com/watch?vp8u_k2LIZyo
6计算机中的浮点数(一)_浮点表示法-CSDN博客
7计算机中的浮点数(二)-CSDN博客 版权声明所有的笔记可能来自很多不同的网站和说明在此没法一一列出如有侵权请告知立即删除。欢迎大家转载但是如果有人引用或者COPY我的文章必须在你的文章中注明你所使用的图片或者文字来自于我的文章否则侵权必究。 ----松下J27
