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

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

数值计算 --- 平方根倒数快速算法(上)_开根号倒数快速-CSDN博客文章浏览阅读712次，点赞31次，收藏30次。由于平方根倒数快速算法实在是太过经典，出于对code中magic number"0x5f3759df"的好奇，出于对John Carmack“what the fuck!”这一锐评的热爱，还有出于对于这一算法的原作者Greg Walsh的致敬。我写了这篇文章，作为该算法的自学笔记。_开根号倒数快速https://blog.csdn.net/daduzimama/article/details/142306873

数值计算 --- 平方根倒数快速算法(中)-CSDN博客文章浏览阅读752次，点赞17次，收藏27次。由于平方根倒数快速算法实在是太过经典，出于对code中magic number"0x5f3759df"的好奇，出于对John Carmack“what the fuck!”这一锐评的热爱，还有出于对于这一算法的原作者Greg Walsh的致敬。我写了这篇文章，作为该算法的自学笔记。https://blog.csdn.net/daduzimama/article/details/142444260

在我的上一篇文章中，我停在了神秘数---5f3759df这一关乎到如何才能尽可能精确且快速的计算出用于NR-iteration初值的常数的讨论。在这篇文章中，我将继续讨论“5f3759df”这个数。这个让人着迷的神秘数究竟是怎么来的已经无从得知了(根据我在网上查到的现有资料中还没有发现)。但，既然code中明确用到了这个数，我在这里可以简单的“就数论数”。看看这个数究竟对初值的计算究竟有什么帮助。

注意：为了便于行文顺畅，以及对照前后文看时不容易引起误解和混淆，我这里所有公式的编号都将延续上文公式的index。但会给这篇文章中的新公式给予新的index。本来这一系列的文章我最初都是放在一篇文章中写的，但那样的文章实在是太过于冗长，这才分成了上中下三篇。

1，再论log2(1+M)的近似

为了说明这个问题，我们先回到上篇文章中的“log2(1+M)的近似”。在上一篇文章中我写道：

由于 $M_{x}$ 是一个0~1之间的小数。以M为自变量，当M=0~1时，以M为自变量的函数y=log2(1+M)与函数y=M的函数值差异很小。故而，有如下近似：

${log_{2}}^{(1+M_{x})}\approx M_{x}$ ， $when \; M_{x}=0\sim 1$ ，(式10)

基于这一近似，再加上后续一系列的推导依次得到：

${log_{2}}^{x}\approx 1/2^{-23}\cdot (E_{x}\times 2^{23}+T_{x})-127$ ，(式12)

${log_{2}}^{x}\approx x_{B}/2^{23}-127$ ，(式14)

${log_{2}}^{a}=-1/2\cdot {log_{2}}^{x}$ ，(式16)

$a_{B}=381\times 2^{22}-x_{B}/2$ ，(式17)

$a_{B}=5f400000-x_{B}/2$ ，(式18)

如果详细观察上述近似的函数图像可知，这种近似只有在函数两端的近似效果比较好，越是到函数中部，误差越大。

为了改善这一问题，我们对原来的近似做一个修正。令y=M的函数上移一部分，即，让函数加上一个很小的正数σ，得到如下函数和对应的函数图像(下图中的黄线)：

$y=M+\sigma$

2，σ是多少才合适呢？

基于上面的结论，我们为了尽可能均衡的分散计算误差，选择了用M+σ的近似替换原来的近似得到：

${log_{2}}^{(1+M_{x})}\approx M_{x}+\sigma$ ， $when \; M_{x}=0\sim 1$ ，(式19)

原来的近似是：

${log_{2}}^{(1+M_{x})}\approx M_{x}$ ， $when \; M_{x}=0\sim 1$ ，(式10)

Tips:为了加以区别，文中引用的所有前文中的公式，都沿用前文中所使用的背景颜色(淡红色或黄色)。而，本文中新公式的背景色都是淡蓝色。

基于新的近似(式19)，浮点数x的对数表达式(式9)变为：

${log_{2}}^{x}={log_{2}}^{(1+M_{x})}+E-127\approx M_{x}+\sigma +E_{x}-127$

得到：

${log_{2}}^{x}\approx M_{x}+\sigma +E_{x}-127$ ，(式20)

同样的，把之前的整体替换(式8)代到(式20)中：

$M_{x}=T_{x}\cdot 2^{-23}$ ，(式8)

(式20) 变为：

${log_{2}}\approx M_{x}+\sigma+E_{x}-127=T_{x}\cdot 2^{-23}+\sigma+E_{x}-127=1/2^{23}\cdot (E_{x}\times 2^{23}+T_{x})+\sigma-127$

得到：

${log_{2}}\approx 1/2^{23}\cdot (E_{x}\times 2^{23}+T_{x})+\sigma-127$ ，(式21)

再把浮点数x在计算机中的存储形式(式13)代入上式：

$x_{B}=E_{x}\times 2^{23}+T_{x}$ ，(式13)

(式21) 变为：

${log_{2}}\approx 1/2^{23}\cdot (E_{x}\times 2^{23}+T_{x})+\sigma-127=1/2^{23}\cdot x_{B}+\sigma-127$

得到:

${log_{2}}\approx x_{B}/2^{23}+\sigma-127$ ，(式22)

把(式22)代入(式16)的左右两边：

${log_{2}}^{a}=-1/2\cdot {log_{2}}^{x}$ ，(式16)

左边：

${log_{2}}^{a}=a_{B}/2^{23}+\sigma-127$

右边：

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

最终(式16)变为：

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

得到：

$a_{B}=3\times (127-\sigma )\times 2^{22}-x_{B}/2$ ，(式23)

令(式23)中的“ $3\times (127-\sigma )\times 2^{22}$ ”等于十六进制的5f3759df：

$3\times (127-\sigma )\times 2^{22}=5f3759df_{16}$

最终求出函数的上移幅度σ：

$\sigma =127-5f3759df_{16}/3/(2^{22})=127-1597463007_{10}/3/(2^{22})=0.04504656791687012$

把这一结果用matlab画出来得到新的近似函数图像和均匀分布后的误差：

通过把原来的近似函数上移了一个约等于0.045的σ，近似函数与原函数的误差从原来的中间误差特别大而两端误差很小，变成了较为均衡的分布，并在两个点处达到最小。

3，还有比“5f3759df”更好的选择吗？

由于这套算法被很多人使用和研究过，又因为这套代码最令人着迷的就是那个神秘数。因此，在我自己学习这套代码的时候，除了之前提到的粗略近似所使用的常数“5f400000”和Grag Walsh文中给出的“5f3759df”，还有一个比较著名的数学家Chris Lomont也提出了他认为较为理想的常数。并通过理论计算和穷举法分别找到了0x5f37642f和0x5f375a86。

此外，在他的论文中，Lomont亦指出了64位的IEEE754浮点数（即双精度类型）所对应的魔术数字是0x5fe6ec85e7de30da，但后来的研究表明0x5fe6eb50c7aa19f9的结果精确度更高。McEniry得出的结果则是0x5FE6EB50C7B537AA，精度介于两者之间，Matthew Robertson给出的精度更高的值是0x5FE6EB50C7B537A9。

这里我想通过3个数的实验来结束关于这一问题的讨论。并基于这三个数的实验看看“5f400000”，“5f3759df”，“5f37642f”和“5f375a86”这四个数究竟谁的精度更好。这里稍微提一嘴，由于“5f3759df”，“5f37642f”和“5f375a86”这三个数都小于“5f400000”，这说明前面三者干的都是一件事，都在把“5f400000”所对应的直线向上做平移。

我们分别以x=2,3,4为例，即用上边四个常数分别求 $1/\sqrt{2},1/\sqrt{3},1/\sqrt{4}$ ，并且给出初值和迭代1，2，3次后的结果，看看究竟哪个常数好？

# include <stdio.h>
# include <math.h>

void rsqrt_5f400000(float number,float *y)
{
	long i;
	float x2;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y[0] = number;
	i = *(long*)&y[0];                       // evil floating point bit level hacking
	i = 0x5f400000 - (i >> 1);
	y[0] = *(float*)&i;
	y[1] = y[0] * (threehalfs - (x2 * y[0] * y[0]));   // 1st iteration
	y[2] = y[1] * ( threehalfs - ( x2 * y[1] * y[1]) );   // 2nd iteration, this can be removed
	y[3] = y[2] * (threehalfs - (x2 * y[2] * y[2]));
}

void rsqrt_5f3759df(float number, float* y)
{
	long i;
	float x2;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y[0] = number;
	i = *(long*)&y[0];                       // evil floating point bit level hacking
	i = 0x5f3759df - (i >> 1);
	y[0] = *(float*)&i;
	y[1] = y[0] * (threehalfs - (x2 * y[0] * y[0]));   // 1st iteration
	y[2] = y[1] * (threehalfs - (x2 * y[1] * y[1]));   // 2nd iteration, this can be removed
	y[3] = y[2] * (threehalfs - (x2 * y[2] * y[2]));
}

void rsqrt_5f37642f(float number, float* y)
{
	long i;
	float x2;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y[0] = number;
	i = *(long*)&y[0];                       // evil floating point bit level hacking
	i = 0x5f37642f - (i >> 1);
	y[0] = *(float*)&i;
	y[1] = y[0] * (threehalfs - (x2 * y[0] * y[0]));   // 1st iteration
	y[2] = y[1] * (threehalfs - (x2 * y[1] * y[1]));   // 2nd iteration, this can be removed
	y[3] = y[2] * (threehalfs - (x2 * y[2] * y[2]));
}

void rsqrt_5f375a86(float number, float* y)
{
	long i;
	float x2;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y[0] = number;
	i = *(long*)&y[0];                       // evil floating point bit level hacking
	i = 0x5f375a86 - (i >> 1);
	y[0] = *(float*)&i;
	y[1] = y[0] * (threehalfs - (x2 * y[0] * y[0]));   // 1st iteration
	y[2] = y[1] * (threehalfs - (x2 * y[1] * y[1]));   // 2nd iteration, this can be removed
	y[3] = y[2] * (threehalfs - (x2 * y[2] * y[2]));
}

int main() {
	float x = 2.0f;
	printf("input x=%f\n", x);
	printf("ideal result=%.7f\n", 1 / sqrt(x));

	float y[4] = {0};
	rsqrt_5f400000(x,y);
	printf("calc with 5f400000 = %.7f,%.7f,%.7f,%.7f\n", y[0], y[1], y[2], y[3]);
	rsqrt_5f3759df(x, y);
	printf("calc with 5f3759df = %.7f,%.7f,%.7f,%.7f\n", y[0], y[1], y[2], y[3]);
	rsqrt_5f37642f(x, y);
	printf("calc with 5f37642f = %.7f,%.7f,%.7f,%.7f\n", y[0], y[1], y[2], y[3]);
	rsqrt_5f375a86(x, y);
	printf("calc with 5f375a86 = %.7f,%.7f,%.7f,%.7f\n", y[0], y[1], y[2], y[3]);
    return 0;
}

对 $1/\sqrt{2}$ 而言：

先看初值，“5f400000” 的精度明显要低于后面三者，但抛开精度不谈，这四个常数的近似结果也和理想值0.7071068都差的比较远。经过一次迭代后，由于“5f400000”初值的精度低，他一次迭代后的结果也很低。后面三者的精度也只停留在小数点后两位，其中，“5f3759df”的0.7069300是最接近0.7071068的。如果是按只迭代一次的原程序的话，计算到这里那还是“5f3759df”的精度是最高的，即，最接近正确答案0.7071068。

但经过两轮迭代之后却是"5f375a86"的0.7071068精度最高。最后，经过三次迭代，4个常数都得到了和正确答案一样的结果。

对 $1/\sqrt{3}$ 而言：

就初值来看依然是“5f400000”的精度最低，而且要比后三者低很多。这是符合预期的，毕竟这是最4者中最糙的一个近似。对后三者而言，“5f37642f”的0.5913724是最接近0.5773503。但整体看下来，三者的精度依然不理想。经过两次迭代后“5f3759df”的结果最接近标准答案，但整体的精度也仅能达到小数点后两位。

对 $1/\sqrt{4}$ 而言：

这是一个比较特殊的数，直接能够得到精确答案0.5。从程序运行的结果上看“5f400000”这个常数对这种计算最敏感，计算初值的时候就直接得到了标准答案，后面的多次迭代也是标准答案。

就后三者而言，“5f37642f”的初值0.4831862最接近0.5，但三者初值的精度依旧都不高。

小结：

经过我简单的实验，我得到如下结论：

1，四个常数初值的计算结果都不太理想(除了“5f400000”在计算结果为可整除的小数时)。因此，后续的NR-iteration是必不可少的。

2，由于源代码中只算了一次NR-iteration，单看这一次的结果，平均下来也就只能精确到小数点后两位。

3，实验中所用到的这四个常数，除了“5f400000”。基本上在迭代一次后都能达到差不多的效果，不存在那个比另一个精度高很多的情况。因此，Greg walsh在code中所使用的magic number---“5f3759df”能够出色的完成任务，即便是后续有“更好的常数”，在这个算法框架下相较于“5f3759df”也很难提高太多。

4，Greg walsh的这套平方根倒数快速算法固然快，也设计的非常精妙，同时用到了很多学科的知识点。但如果只用一次NR-iteration，计算结果的精度依然非常有限。

（全文完）

--- 作者，松下J27

参考文献(鸣谢)：

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

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

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

4，揭秘·变态的平方根倒数算法

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

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

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

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