自学内容网 自学内容网

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

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

1,平方根倒数快速算法是如何选择初值的?WTF中的神秘数字究竟是怎么来的?

        花开两朵,各表一枝。在前面的介绍中,我们已经知道了这段代码的作者在函数的最后使用了NR-iteration,且只用了一次NR-iteration。这样一来,选择正确的初值就显得尤为重要了。在源代码中,求解NR-iteration所需初始值的关键在于充分的利用浮点数x在计算机中的表示/编码方式

        由于code中的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. 标准化二进制小数:将二进制小数表示成规范形式。规范形式要求小数点左侧只能有一位,且必须是1,因此:0.01010101..._{2}=1.01010101..._{2}\times 2^{-2}

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

  3. 有效数字的尾数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 all

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

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 hacking
i = 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 removed

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

 相应的输出为:

        就本例而言,二者的计算结果都非常接近准确值1,但5f400000的精度要更高,5f3759df的误差约为0.002。 如果以x=4为例,准确值为0.5,再看看二者的表现:

        结果还是5f400000的准确性更高,5f3759df的误差约为0.0009。但上面的两个例子都是平方根的结果正好是整数的情况,例如\sqrt{1}=1\sqrt{4}=2。但如果碰到平方根为无理数的情况呢,我们分别试试x=2和x=3的情况。

        有趣的是,在这两个例子中基于magic number的计算结果要比5f400000的精度高。对于x=2而言,5f400000的误差约为0.004,5f3759df的误差约为0.0002。 对于x=3而言,5f400000的误差约为0.006,5f3759df的误差约为0.0005。 

        这样看来,不论是采取哪种常数去估算初值,基本上都已经能够得到较为准确的结果。毕竟,这一初值还要用牛顿拉夫逊法再迭代一次才是最终的结果。


(全文完) 

--- 作者,松下J27

参考文献(鸣谢):

1,https://en.wikipedia.org/wiki/Newton%27s_method#Examples

2,什么代码让程序员之神感叹“卧槽”?改变游戏行业的平方根倒数算法_哔哩哔哩_bilibili

3,[算法] 平方根倒数速算法中的魔数0x5f3759df的来源 | GeT Left

4,https://i.hsfzxjy.site/uncover-the-secret-of-fast-inverse-square-root-algorithm/

5,https://www.youtube.com/watch?v=p8u_k2LIZyo

6,计算机中的浮点数(一)_浮点表示法-CSDN博客

7,计算机中的浮点数(二)-CSDN博客 

8,Beyond3D - Origin of Quake3's Fast InvSqrt()

9,Beyond3D - Origin of Quake3's Fast InvSqrt() - Part Two

版权声明:所有的笔记,可能来自很多不同的网站和说明,在此没法一一列出,如有侵权,请告知,立即删除。欢迎大家转载,但是,如果有人引用或者COPY我的文章,必须在你的文章中注明你所使用的图片或者文字来自于我的文章,否则,侵权必究。 ----松下J27  


原文地址:https://blog.csdn.net/daduzimama/article/details/142444260

免责声明:本站文章内容转载自网络资源,如本站内容侵犯了原著者的合法权益,可联系本站删除。更多内容请关注自学内容网(zxcms.com)!