快速平方根倒数算法

3 minute read

Published:

快速平方根倒数算法

引入

计算机图形学中经常需要计算一个数的平方根的倒数,例如求向量 $(x,y,z)$ 的单位向量: \(x_0 = \frac{1}{\sqrt{x^2+y^2+z^2}} \cdot x \\ y_0 = \frac{1}{\sqrt{x^2+y^2+z^2}} \cdot y \\ z_0 = \frac{1}{\sqrt{x^2+y^2+z^2}} \cdot z \\\) 如今的英特尔芯片以及英伟达的 CUDA 库都已经有了硬件集成好的算法可以调用,而在上世纪九十年代,该计算过程一般用牛顿迭代法实现,即循环迭代去逼近解: \(y = f(x) = \frac{1}{\sqrt{x}} \ \ \ \Rightarrow \ \ \ x = \frac{1}{y^2} \\ 令 g(y) = \frac{1}{y^2}-x,则 g^{'}(y) = -2\cdot y^{-3} \\ 牛顿迭代法: y_{n+1} = y_n - \frac{g(y_n)}{g^{'}(y_n)} = y_n + \frac{y_n-xy^3_{n}}{2} \\ \Rightarrow y_{n+1} = y_n(\frac{3}{2}-\frac{1}{2}xy^2_n)\)

float rsqrt(float x) {
    float y = 1.0f / x;
    float err = 1e-5;
    float yy;
    do {
        yy = y;
        y = y * (1.5f - 0.5f * x * y * y);
    } while(fabs(y-yy) >= err);
    return y;
}

循环迭代计算的计算速度并不出色。2005 年,人们在游戏 Quake 3 开源代码中发现了一个快速计算平方根倒数的算法,不过经过一些考证之后(详见维基百科)发现此算法在上世纪九十年代就已经出现,对于此算法的最初发明者,目前也没有确凿的定论。以下是《雷神之锤3》中的相关代码原文:

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;
}

IEEE 754

C 语言中的 float 类型存储遵循 IEEE 754 标准,如下图所示:

IEEE754

  • float 类型变量大小为 32 位,其中第 1 位用来表示浮点数的符号,记为 S0 表示正数, 1 表示负数。
  • 之后的 8 位是浮点数的阶码,记为 e,存储的是移码,在这里,一个数的移码就是这个数加上 127
  • 最后的 23 位是浮点数的尾数,采用规格化二进制小数表示,记为 m,简单来讲,任何一个二进制小数都可以写成科学计数法的形式,例如 $1.1001001 \times 2^7 $,其中,整数位的那个 1 我们不存进内存,因为科学计数法的形式下,整数位肯定是 1 (如果不是,就一直右移小数点直至整数位为 1),这样就可以多存一位尾数了,提高了存储精度。

综上,按照 IEEE 754 标准,这个浮点数 x 可以写成: \(x = (-1)^S\cdot(1+m)\cdot 2^e\) 这里的小写字母表示的数都可以按照数学上的数字直接理解。

假如我们不考虑浮点数的格式之类的,就把内存中的 01 串当做整数来对待,那么 e 所在内存表示为整数就是 \(E = e + 127 \\ 例如 e = 3 = (0000 \ 0011)_2 \\ 则 E = e + 127 = (0000 \ 0011)_2 + (0111 \ 1111)_2 = (1000 \ 0010)_2\) 即把实际的值 e 加上 127 存储为移码的格式。

m 所在内存表示为整数就是 \(M = m \times 2^{23} \\ 例如 m = 0.11110000111100001111010 \\ 则 M = 11110000111100001111010 , 就是去掉小数点\\\) 另外,在这个算法中,因为是算平方根的倒数,所以这里的数都是正数,因此 S = 0

综上所述,我们有了: \(x = (1+m)\cdot 2^e \\ M = m\cdot L, L = 2^{23} \\ E = e + B, B = 127 \\\)

对数化简和函数近似

对于一个单精度浮点数 x,我们要计算的是 \(y = \frac{1}{\sqrt{x}}=x^{-\frac{1}{2}} \\ 两边同时取以 2 为底的对数 \\ \log_{2}y = -\frac{1}{2}\log_{2}x & (1)\\ 而 \begin{cases} x = (1+m_x)\cdot 2^{e_x} \\ y = (1+m_y)\cdot 2^{e_y} \\ \end{cases} & (2) \\ (2)代入(1)所以有 \\ \log_{2}((1+m_y)\cdot2^{e_y}) = -\frac{1}{2}\log_{2}((1+m_x)\cdot2^{e_x}) \\ 化简为 \\ \log_{2}(1+m_y)+e_y = -\frac{1}{2}(\log_{2}(1+m_x)+e_x) & (*)\) (*)式中含有对数,不好处理,考虑到 $m_x,m_y$都是$(0,1)$之间的小数,我们可以考虑用一个线性函数 $f(z)=z+b$去近似表示$g(z)=\log_{2}(1+z)$,只要 b 选择足够合适,就可以有较高的近似精度。

log

现在将(*)式中的对数函数用近似的线性函数替换: \(\log_{2}(1+m_y)+e_y = -\frac{1}{2}(\log_{2}(1+m_x)+e_x) & (*) \\ \bbox[#DB7D74, 5px]{\log_{2}(1+m_y)}+e_y = \bbox[#576690, 5px]{m_y+b} + e_y \\ -\frac{1}{2}( \bbox[#DB7D74, 5px]{\log_{2}(1+m_y)}+e_y) = -\frac{1}{2}(\bbox[#576690, 5px]{m_x+b} + e_y) \\ 所以m_y+b+e_y = -\frac{1}{2}(m_x+b+e_x) & (**)\) 接下来,我们将(**)式子中的小写字母变量改为用大写字母变量来表示: \(M_x = Lm_x, \ M_y = Lm_y \\ E_x = e_x + B, \ E_y = e_y + B \\ L = 2^{23}, \ B = 127\) 代入 (**) 式: \(\frac{M_y}{L} + b + E_y - B = -\frac{1}{2}(\frac{M_x}{L} + b + E_x - B) \\\) 化简为: \(LE_y+M_y = \frac{3}{2}L(B-b)-\frac{1}{2}(LE_x+M_x) \tag{***}\)

神之一手

观察(***)式的形式,我们会发现两边都有 $LE+M$ 的形式。

如果不考虑浮点数的存储格式,将一个 32 位的 float 变量看做一个 32 位的整数来看待,那么会发现 $LE+M$ 就是将其看作整数之后的那个整数的值(这里S=0故不用考虑)。

现在有两个 float 变量,如果按照 32 位单精度浮点数看待,值就是 $x$ 和 $y$;如果按照 32 位的整型看待,就是 $I_x,I_y$: \(I_x=LE_x+M_x \\ I_y = LE_y+M_y \\\) 代入(***)式,得到: \(I_y = \frac{3}{2}L(B-b)-\frac{1}{2}I_x\)

The Fucking Number 魔数

我们已经推导到了 \(I_y = \frac{3}{2}L(B-b)-\frac{1}{2}I_x\) 这里 LB 都是确定的常数,因此只要确定了 合适的 b 的值,就能在很小的误差之下求出平方根的倒数。为了表示方便,我们将这几个常数合并为 \(magic = \frac{3}{2}L(B-b) = \frac{3}{2}\cdot 2^{23}\cdot (127-b)\) 那么就有 \(I_y = magic \ - \ \frac{1}{2}I_x\) 在已知浮点数 $x$ 的前提下,将其按照整数形式读出为 $I_x$,根据这个公式就可以计算出 $I_y$,再将 $I_y$以单精度浮点数的形式读出,即可求得结果 $y$。

只要找到一个合适的 magic 值,就能计算出倒数平方根。在游戏《雷神之锤3》的源代码中,给出了这个神奇的数字:

0x5f3759df

后来,人们找到了一个更加精确的值:

0x5f375a86

代码

float fast_rsqrt(float x) {
	unsigned int Ix = *(unsigned int*)(&x);
	Ix  = Ix >> 1;
    const unsigned int MAGIC = 0x5f3759df;
	unsigned int Iy = MAGIC - Ix;
	float y = *(float*)(&Iy);
	y = y * (1.5f - 0.5*x*y*y);
	y = y * (1.5f - 0.5*x*y*y);
	y = y * (1.5f - 0.5*x*y*y);
	return y;
}