平方根倒数速算法（Fast Inverse Square Root）

最新推荐文章于 2024-07-18 22:34:36 发布

winkemoji

最新推荐文章于 2024-07-18 22:34:36 发布

阅读量2.1k

点赞数 5

文章标签：算法

本文链接：https://blog.csdn.net/weixin_43899202/article/details/120540104

版权

本文深入解析了平方根倒数速算法（FastInverseSquareRoot），该算法源于90年代，最初用于《雷神之锤III竞技场》。通过浮点数位操作和牛顿迭代法，实现了高效计算32位浮点数的平方根倒数。核心步骤包括浮点数的二进制表示与对数关系，以及利用位移运算进行快速除法。此外，还介绍了算法的优化迭代过程，以提高精度。

摘要由CSDN通过智能技术生成

平方根倒数速算法（Fast Inverse Square Root）

一、背景

平方根倒数速算法是适用于快速计算积的算术平方根的倒数（在此需取符合IEEE 754标准格式的32位浮点数）的一种算法。

此算法最早可能是于90年代前期由SGI所发明，后来于1999年在《雷神之锤III竞技场》的源代码中应用，但直到2002－2003年间才在Usenet一类的公共论坛上出现。这一算法的优势在于减少了求平方根倒数时浮点运算操作带来的巨大的运算耗费，而在计算机图形学领域，若要求取照明和投影的波动角度与反射效果，就常需计算平方根倒数。
🔔 以上内容来自百度百科

源代码 source code

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 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 ) );		// 2st iteration, can be removed
    
    return y;
}

二、需要知道的知识😢

1. IEEE 754 单精度浮点数标准

2. 牛顿迭代法

3. C语言中移位运算符(>>)

4. 当x在区间[0,1]内，有 $\log_2(1+x) \approx x$

在这里插入图片描述

根据上图可知，当 $\in [0,1)$ ,使用 $y = x$ 近似的表示 $log_2(1+x)$ 好一点，所以有公式(1),其中 $\mu$ 为校正系数。
$\log_2(1+x) = x + \mu \tag{1}$

三、平方根倒数速算法推导

源代码中求平方根倒数主要步骤如下：

    i = * ( long* ) &y;							// 浮点数不能做位操作，&y获取浮点数y的地址，( long* )&y将该地址转换为长整型，最终取到该值
    i = 0x5f3759df - (i >> 1);					// 公式(6)，（i >> 1 右移将该数减半，相比于做除法速度更快）
    y = * ( float * ) &i;						// 将长整型地址转换为浮点数地址
    y = y * (threehalfs - ( x2 * y * y ) );		// 公式(7)

假设开平方根倒数为 $a$ (a>0), 结果为 $\tau$ 则有：
$\tau = \frac{1}{\sqrt{a}} \tag{2}$
根据IEEE754浮点数标准，浮点数 $x$ 其对应的二进制b表示（除去最高符号位）与其指数 $E$ （8个二进制位）, 尾数M（23个二进制位）之间的关系为:
$2^{23}E + M \tag{3}$

$x$ 与指数 $E$ , 尾数 $M$ 关系为：
$\frac{M}{2^{23} }) \cdot 2^{(E - 127)} \tag{4}$
公式 $(4)$ 两端取对数如下：

$\begin{aligned} \log_2{x} &= \log_2 {\{(1 + \frac{M}{2^{23} }) \cdot 2^{(E - 127)}\}}\\ &=\log_2 {(1 + \frac{M}{2^{23}})} +E - 127 \\ &= \frac{M}{2^{23} } + \mu + E -127 \\ &= \frac{1}{2^{23}}(M + 2^{23}E) + \mu -127 \\ &=\frac{1}{2^{23}} b + \mu -127\tag{5} \end{aligned}$

由公式(5)可知浮点数的二进制表示 $b$ 在某种意义上是浮点数 $x$ 的对数，所以相比于求 $\frac{1}{\sqrt{a}}$ , 直接对其取对数并计算出结果 $\tau$ 的二进制表示,对式(2)两端取对数见式(6)。

$\begin{aligned} \log_2 \tau &= \log_2 (\frac{1}{\sqrt{a}}) \\ \log_2 \tau&=-\frac{1}{2}\log_2 {a} \\ \frac{1}{2^{23}} b_{\tau} + \mu -127 &= -\frac{1}{2}( \frac{1}{2^{23}} b_a + \mu -127 ) \\ \frac{1}{2^{23}}b_{\tau} &=- \frac{1}{2^{24}} b_a - \frac{3}{2} \mu + \frac{3 \times 127}{2} \\ b_\tau &={3}\times 2^{22} (127-\mu) - \frac{1}{2} b_a \tag{6} \end{aligned}$
公式(6)得到 $b_\tau ={3}\times 2^{22} (127-\mu) - \frac{1}{2} b_a$ 其中 ${3}\times 2^{22} (127-\mu)$ 就是神秘数字 0x5f3759df 😆 ，校正系数 $\mu$ 可以取值0.043。

将 $b_{\tau}$ 转换为对应浮点数 $\tau$ ，同时为进一步减少误差，使用牛顿迭代法继续提高精度,假设关于 $\tau$ 的方程 $f(\tau)$ 并使其为0：
$f(\tau) = \tau - \frac{1}{\sqrt{a}} = 0 \Longrightarrow \frac{1}{\tau^2} - a = 0$

根据牛顿迭代法可得：
$\begin{aligned} \tau_{new} &= \tau - \frac{f(\tau)}{f^{'}{(\tau)}} \\ &=\tau - \frac{\frac{1}{\tau^2} - a}{-2 \frac{1}{\tau^3}} \\ &=\tau + \frac{1}{2} (1- a\tau^2)\tau \\ &= \tau(\frac{3}{2}-\frac{1}{2} a \tau^2) \end{aligned} \tag{7}$
证明完毕。👏