算法学习笔记（二）：平方根倒数速算法

大家好，又见面了，我是你们的朋友全栈君。

序

这是一个神奇的算法！

一、介绍

起源于一篇《改变计算技术的伟大算法》文章，知道这个算法，然后google一下，维基讲的还不错，本文权当自己理清下思路。先贴源代码，为《雷神之锤III竞技场》源代码中的应用实例，剥离了C语言预处理器的指令，并附上了原有的注释。

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（对浮点数的邪恶位级hack）
	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;
}

算法的思路：

1、计算出该浮点数的平方根倒数的首次近似值，见源代码中的8-11行；

2、利用牛顿法迭代以加强精度，得到要求精度内的值（迭代次数根据精度要求调整，源代码中一次迭代就满足精度要求）。

算法的巧妙之处在于代码中的四行蓝色代码，大致过程是：将表示浮点数的字节序列用来表示整数（i = * ( long * ) &y;），然后通过一个巧妙的整数运算得到一个新的整数（i = 0x5f3759df – ( i >> 1 );），最后将表示整数的字节序列换回表示浮点数。所以弄清算法关键障碍是：在计算机中是如何表示浮点数和整数的、整数运算又怎能算出浮点数的平方根倒数的近似值、0x5f3759df怎么来的。

二、浮点数和整数的表示方法

浮点数和整数存储位数一定是相同的。这里浮点数和整数都占4字节，32位。

一个浮点数是由32位二进制位表示的有理数，分为三部分。其中符号占1位，表示正负，记为Si；指数占接下来的8位，表示经过偏移处理后的指数，即实际表示E（如图中为124），需要偏移B（图中为2的8次方减1,127。B为一个固定值），最后得指数值为E-B；有效数字（除最高位以外）占剩下的23位，记为m（0<m<1）,图中的

\scriptstyle m=1\times 2^{-2}=0.250

。

所以浮点数的结构公式为：

\scriptstyle x=(-1)^{\mathrm{Si}}\cdot(1+m)\cdot 2^{(E-B)}

， 图中

\scriptstyle x=(1+0.250)\cdot 2^{-3}=0.15625

算法学习笔记（二）：平方根倒数速算法

整数的表示相对简单，符号占1位，数值占剩下的31位。如果用上图的浮点数字节序列来表示整数，那么

\scriptstyle I=E\times 2^{23}+M

，即

I=124\times 2^{23} + 2^{21}

.平方根倒数函数仅能处理正数，所以符号位均为0。

小结：对于同样的32位二进制数码，若为浮点数表示时实际数值为

\scriptstyle x=(1+m_x)2^{e_x}

，而若为整数表示时实际数值则为

\scriptstyle I_x=E_xL+M_x

，其中

\scriptstyle L=2^{n-1-b}

，这里n=32，b=8。式子中引入的新变量为：

m_x=\frac{M_x}{L}

——-——-—–—–—–—–—–—–等式1

e_x=E_x-B

，其中

B=2^{b-1}-1

————等式2

三、浮点数的平方根倒数近似值

理解浮点数和整数的表示后，下面开始推导。

平方根倒数方程为：

y=\frac{1}{\sqrt{x}}

两边取对数有：

\log_2{(y)}=-\frac{1}{2}\log_2{(x)}

因为浮点数可表示为：

\scriptstyle x=(1+m_x)2^{e_x}

，所以也有

算法学习笔记（二）：平方根倒数速算法

，代入上式有：

\log_2(1+m_y)+e_y=-\frac{1}{2}\log_2{(1+m_x)}-\frac{1}{2}e_x

再度引入新数

\sigma

描述

\scriptstyle \log_2{(1+x)}

与近似值R间的误差：由于

\scriptstyle 0 \le x < 1

，有

\scriptstyle \log_2{(1+x)}\approx {x}

，则在此可定义

\sigma

与x的关系为

\scriptstyle \log_2{(1+x)}\cong x+\sigma

，其中

\sigma

介于0到1/3，所以将

\scriptstyle \log_2{(1+x)}= x+\sigma

代入上式得：

m_y+\sigma+e_y=-\frac{1}{2}m_x-\frac{1}{2}\sigma-\frac{1}{2}e_x

将第二部分小结中等式1，等式2代入到上述方程中，有：

M_y+(E_y-B)L=-\frac{3}{2}\sigma{L}-\frac{1}{2}M_x-\frac{1}{2}(E_x-B)L

移项整理得：

E_yL+M_y=\frac{3}{2}(B-\sigma)L-\frac{1}{2}(E_xL+M_x)

又因为浮点规格存储的正浮点数x，若将其作为长整型表示，则示值为：

\scriptstyle I_x=E_xL+M_x

，所以x的平方根倒数的首次近似值的整数表示值为：

I_y=E_yL+M_y=R-\frac{1}{2}(E_xL+M_x)=R-\frac{1}{2}I_x

，其中

R=\frac{3}{2}(B-\sigma)L

，

B=2^{b-1}-1

，

\scriptstyle L=2^{n-1-b}

，

n=32，b=8。

这个式子对应着源代码中的这一行：i = 0x5f3759df – ( i >> 1 );，然后将整数表示值换回表示浮点数：y = * ( float * ) &i;。这样就得到了浮点数的平方根倒数的近似值。

四、神秘的0x5f3759df

由第三部分可知：0x5f3759df 对应着R，即3/2(B-

\sigma

)L.当R为0x5f3759df时，有

\scriptstyle \sigma=0.0450461875791687011756

.

“现在不仅该算法的原作者不明，人们也仍无法明确当初选择这个“魔术数字”的方法。Chris Lomont在研究中曾做了个试验：他编写了一个函数，以在一个范围内遍历选取R值的方式将逼近误差降到最小，以此方法他计算出了线性近似的最优R值0x5f37642f（与代码中使用的0x5f3759df相当接近），但以之代入算法计算并进行一次牛顿迭代后，所得近似值与代入0x5f3759df的结果相比精度却仍略微更低。……在Charles McEniry的论文中，他使用了一种类似Lomont但更复杂的方法来优化R值：他最开始使用穷举搜索，所得结果与Lomont相同；而后他尝试用带权二分法寻找最优值，所得结果恰是代码中所使用的魔术数字0x5f3759df”—维基百科

五、结束

至于最后一步牛顿法，自行google。平方根倒数速算法的神奇之处在于：1、充分利用了浮点数和整数在计算机中的表示，然后以两次转换表示和一次整数运算替换复杂的浮点数计算，最后通过牛顿法加强精度；2、R的取值。

参考资料

http://zh.wikipedia.org/wiki/%E5%B9%B3%E6%96%B9%E6%A0%B9%E5%80%92%E6%95%B0%E9%80%9F%E7%AE%97%E6%B3%95

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/140041.html原文链接：https://javaforall.cn