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

🕗 发布于 2024-09-23 10:42 牛顿法牛顿拉夫逊平方根倒数 John Carmack NR-iteration

平方根倒数快速算法 --- 向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，因此： $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 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

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

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

上一篇：Java接口详解
下一篇：Springboot指定扫描路径

【已为网站上传证书,却显示不安全】
【已为网站上传证书,却显示不安全】
阅读更多2024-11-16
Jtti：服务器总是自动重启怎么办?
1.查看系统日志：通过查看系统日志文件(如/var/log/syslog、/var/log/messages或Windows Event Viewer)获取错误信息，找出重启原因。解决方案：查看系统日
阅读更多2024-11-16
Scala的Array(1)
Scala的Array表示长度不可变的数组，若需要定义可变数组需要倒包。
阅读更多2024-11-16
Scala中的Array
Array是一种在Scala中用于存储固定大小，相同类型元素的数据结构，Array是可变。类型有两种：Array表示长度不可变的数组，和ArrayBuffer表示长度可变的数组。(一)Array的定义
阅读更多2024-11-16
llama factory lora 微调 qwen2.5 7B Instruct模型
项目背景甲方提供一台三卡4080显卡需要进行qwen2.5 7b Instruct模型进行微调。以下为整体设计。要使用对进行，流程与之前提到的模型类似。LoRA 微调是一种高效的微调方法，通过低秩
阅读更多2024-11-16
大数据中的时序数据，以及时序数据架构设计解决方案
时序数据架构设计需根据具体的应用场景、数据量和实时性要求来选择合适的方案。如果数据量较小且不需要实时处理，关系型数据库可能足够；如果数据量大且需要高效实时查询，可以选择时序数据库；而对于大规模分布式处
阅读更多2024-11-16
Linux网络——网络初识
在谈及网络之前，我们要先对学的知识有一个了解网络的本质就是获取数据，而系统的本质就是加工数据。
阅读更多2024-11-16
革新预测领域：频域融合时间序列预测，深度学习新篇章，科研涨点利器
2024发论文&模型涨点之——频域+时间序列频域+时间序列不仅能提供更丰富的信息，还能提高模型性能和预测准确性。对于论文er来说，是个可发挥空间大、可挖掘创新点多的研究方向。当前，学术界的目光
阅读更多2024-11-16
本地启动浏览器，并禁用web安全性，解决本地启动时，服务端强制要求https协议导致请求不通的问题
本地启动浏览器，并禁用web安全性，解决本地启动时，服务端强制要求https协议导致请求不通的问题
阅读更多2024-11-16
推荐一款来自韩国的影音全能播放器：KMPlayer播放器
来自韩国的影音全能播放器，与Mplayer一样从linux平台移植而来的Kmplayer(简称KMP)几乎可以播放您系统上所有的影音文件。通过各种插件扩展KMP可以支持层出不穷的新格式。强大的插件功能
阅读更多2024-11-16

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

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

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

1. 十进制转二进制

2. IEEE 754 浮点数表示

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

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

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

1. 十进制转二进制

2. IEEE 754 浮点数表示

相关文章