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

平方根倒数快速算法 --- 向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;
}

相应的输出为：