平方根（倒数）快速算法

平方根倒数快速算法

平方根常出现在游戏的图形计算中，尤其是求一个向量的基向量时，还会以倒数的形式出现，所以如何快速求解平方根一度成为了游戏引擎优化中的一个难题。

本文会先介绍原理，最后放出实现的代码，包括手动实现和使用内联汇编一步到位。
如果对汇编语言感兴趣，欢迎移步2万字汇编语言学习笔记和指令集和现代x86汇编

引入：约翰卡马克的代码

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

原理

当知道一个数的值的时候，你就知道这个数的值；当知道一个数的精确近似值时，你就知道这个数的精确近似值。（原地tp）

牛顿迭代法需要迭代的次数取决于给定的近似值的近似程度，所以需要想办法给出更加精确的近似值。

浮点数

IEEE754规格化浮点数

有效数字是从小数点后第一位开始的，小数点前隐含一个1，即

$$\begin{gathered} 对任意浮点数x，总有x=(-1{)}^{Signature}\cdot (1+m)\cdot 2^{E-B} \\ 其中constB=127 \end{gathered}$$
上图中可以得知符号部分0指数01111100数字部分(1+0.01000000000000000000000)，即0.15625

求平方根转化为求对数

$$\begin{gathered} y=Sqrt(x)，即y=x^{\frac{1}{2}} \\ {\log }_{2}y={\log }_2{x}^{\frac{1}{2}}=\frac{1}{2}{\log }_{2}x \\ 先尝试对{\log }_{2}x进行变换 \\ 把x=(-1{)}^{Signature}\cdot (1+m)\cdot 2^{E-B}且Signature=0带入{\log }_{2}x，得 \\ {\log }_{2}x={\log }_2[(1+m)\cdot 2^{E-B}] \\ 用m_{bit}表示m的二进制{ \\ }_{（这里由于m表示的是一个小数，而m的二进制本身可以被视为整数，所以需要区分开来）} \\ {\log }_{2}x={\log }_2[(1+\frac{m_{bit}}{2^{23}})\cdot 2^{E-B}] \\ {\log }_{2}x={\log }_2(1+\frac{m_{bit}}{2^{23}})+E-B \\ 令t=\frac{m_{bit}}{2^{23}}，则t\in [0,1] \\ f(t)={\log }_2(1+t)在t\in [0,1]时可近似视为f(t)=t \\ 故{\log }_{2}x\approx \frac{m_{bit}}{2^{23}}+E-B=\frac{1}{2^{23}}\cdot (2^{23}\cdot E+m_{bit})-B \\ 此时凑出来的2^{23}非常巧妙，可以视为给E补23个0{ \\ }_{这里不用把E与E的二进制值区分开是因为E的意义与其二进制值表示的意义都是整数，是相同的} \\ 同时，m_{bit}刚好是23位，两者相加，可以视为整个x的二进制表示 \\ {\log }_{2}x=\frac{1}{2^{23}}{x}_{bit}-B(1) \\ 把(1)式代回z=\frac{1}{2}\cdot {\log }_{2}x，得z=\frac{1}{2}\cdot (\frac{1}{2^{23}}{x}_{bit}-B) \\ 即\frac{1}{2^{23}}{y}_{bit}-B=\frac{1}{2}\cdot (\frac{1}{2^{23}}{x}_{bit}-B) \\ 整理，得y_{bit}=\frac{1}{2}{x}_{bit}+2^{22}\cdot B \end{gathered}$$

此时得到了平方根，再在运算中直接除以这个数就可以了

直接求平方根倒数

运算除法的开销比运算乘法高，能否直接求平方根倒数后再进行乘法呢？
$$\begin{gathered} y=x^{-\frac{1}{2}} \\ {\log }_{2}y=-\frac{1}{2}{\log }_{2}x \\ 把(1)式代入，得 \\ \frac{1}{2^{23}}\cdot y_{bit}-B=-\frac{1}{2}\cdot (\frac{1}{2^{23}}\cdot x_{bit}-B) \\ 整理，得y_{bit}=3B\cdot 2^{22}-\frac{1}{2}{x}_{bit} \\ 到这里就是最终形式了 \end{gathered}$$

代入数据

float

$y_{bit}=381\cdot 2^{22}-\frac{1}{2}{x}_{bit}$

对应的值为0x5F40 0000

double

推广到双精度浮点型变量

$y_{bit}=3069\cdot 2^{51}-\frac{1}{2}{x}_{bit}$

对应的值为0x5FE8 0000 0000 0000

调整

为什么上述计算得到的值为0x5F40 0000，而不是之前提到的代码中的0x5F37 59DF呢

因为在把${\log }_2$近似为一条直线时，产生的误差仅在两端较小，中间较大，为了减小这一误差，便进行了一定的调整，使得结果更贴近真实值

double类型的值有待精确，但是通过单精度浮点型计算得出的结果已经够用了，0x5FE6 8000 0000 0000，这个值的出来的结果很可能不够精确，需要多迭代一次

后续处理

在上述计算过程中使用了近似处理，误差是不可避免的，因此一般还会再进行一次牛顿迭代法求解

牛顿迭代法

$$\begin{gathered} y=f(x)上一点(x_n,f(x_n)) \\ 在该点处作切线，其斜率k=f^{\prime }(x_n) \\ 点斜式:y-f(x_n)=f^{\prime }(x_n)(x-x_n) \\ 令y=0，得 \\ {x}_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)} \\ 在求平方根这一实例中，f(x)=x_n^2-input，f^{\prime }(x)=2x_n \\ {x}_{n+1}=x_n-\frac{x_n^2-input}{2x_n} \\ {x}_{n+1}=\frac{1}{2}(x_n+\frac{input}{x_n}) \\ 在求平方根倒数这一实例中，f(x)=\frac{1}{x^2}-input,f^{\prime }(x)=-2x^{-3} \\ {x}_{n+1}=x_n-\frac{\frac{1}{x_n^2}-input}{-\frac{2}{x_n^3}}=x_n+\frac{x_n-input\cdot x_n^3}{2} \\ {x}_{n+1}=x_n\cdot (\frac{3}{2}-\frac{1}{2}input\cdot x_n^2) \end{gathered}$$

平方根

表达式已经在上述的计算中得到了：$y_{bit}=\frac{1}{2}{x}_{bit}+2^{22}\cdot B$

代入数据，得

$y_{bit}=\frac{1}{2}{x}_{bit}+0x1FC00000$(float)

$y_{bit}=\frac{1}{2}{x}_{bit}+0x1FF8000000000000$(double)

没有严格的计算误差，直接参考了0x5F40 0000与0x5F37 59DF的差值，得到0x1FB7 59DF

普渡大学的数学家Chris Lomont看了以后觉得有趣，决定要研究一下卡马克弄出来的这个猜测值有什么奥秘。Lomont也是个牛人，在精心研究之后从理论上也推导出一个最佳猜 测值，和卡马克的数字非常接近, 0×5f37642f。卡马克真牛，他是外星人吗？传奇并没有在这里结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是卡马克赢了… 谁也不知道卡马克是怎么找到这个数字的。最后Lomont采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数字要好上那么一丁点的数字， 虽然实际上这两个数字所产生的结果非常近似，这个暴力得出的数字是0×5f375a86。

扩展

对数实现

泰勒展开

$$\begin{gathered} \ln x=\ln [(1+m)\cdot 2^{E-B}]=\ln (1+m)+(E-B)\ln 2 \\ \ln (1+m)=\ln [\sqrt{2}\times \frac{\sqrt{2}}{2}\times (1+m)]=\frac{1}{2}\ln 2+\ln [\frac{\sqrt{2}}{2}\times (1+m)] \\ [\frac{\sqrt{2}}{2}\times (1+m)]\in [\frac{\sqrt{2}}{2},\sqrt{2}] \\ 泰勒展开 \\ \ln (1+m)=\frac{1}{2}\ln 2+{\Sigma }_{i=1}^{\infty }[(-1{)}^{i+1}\cdot \frac{1}{i}\cdot (t-1{)}^i],其中t=\frac{\sqrt{2}}{2}\times (1+m) \\ \ln x=\frac{1}{2}\ln 2+{\Sigma }_{i=1}^{\infty }[(-1{)}^{i+1}\cdot \frac{1}{i}\cdot (t-1{)}^i]+(E-B)\ln 2 \end{gathered}$$

但由于泰勒展开的时间复杂度稍稍搞了一点，所以不是在此讨论的对象

快速获得近似值

$$\begin{gathered} {\log }_{2}x={\log }_2(1+m)+E-B \\ 如果1+m越趋近于1，那么结果越准确，越方便估算 \\ 不妨在此引入根号，反正基于上面的讨论，我们能尝试快速获得根号值 \\ {\log }_{2}x=2{\log }_2\sqrt{1+m}+E-B \\ 此时再进行近似操作，有 \\ {\log }_{2}x\approx 2(\sqrt{1+m}-1)+E-B=2\sqrt{1+m}+E-B-2 \\ {\log }_{2}x=4{\log }_2\sqrt{\sqrt{1+m}}+E-B \\ {\log }_{2}x\approx 4(\sqrt{\sqrt{1+m}}-1)+E-B=4\sqrt{\sqrt{1+m}}+E-B-4 \end{gathered}$$

关于近似操作，当$x\in [1,2]$时，总有${\log }_{2}x\ge x-1$，也就是这一操作会得到偏小的值

在实际计算时，获得的根号值偏小，结果会比预计的值偏小，再次开根号获得的值依旧偏小，结果会比开一次根号的值更小，而再开一次根号的开销虽然依旧被缩小了，但还是比较大，所以实际计算中使用开一次根号的值即可

综上，这种方式只能获得整数部分准确的近似值，小数部分从第一位开始就不够准确了，但是我们可以尝试使用迭代方法将结果变准确

代码

实现

#pragma once

namespace math
{
    inline float __fastcall sqrt(float x) {
        const int magic = 0x1FB75A86;
        int bits = *(int*)&x;
        float y;

        bits = magic + (bits >> 1);
        y = *(float*)&bits;
        y = 0.5 * (y + x / y);
        y = 0.5 * (y + x / y);
        return y;
    }
    inline double __fastcall sqrt(double x) {
        // 这个值可以参考下方的值，来变得更精确
        const long long magic = 0x1FF75FFF00000000;
        long long bits = *(long long*)&x;
        double y;

        bits = magic + (bits >> 1);
        y = *(double*)&bits;
        y = 0.5 * (y + x / y);
        y = 0.5 * (y + x / y);
        y = 0.5 * (y + x / y);
        return y;
    }
    inline float __fastcall fastInvSqrt(float x) {
        const int magic = 0x5F375A86;
        const float threehalfs = 1.5, halfX = x * 0.5;
        unsigned int bits = *(unsigned int*)&x;
        float y;

        bits = magic - (bits >> 1);
        y = *(float*)&bits;
        y = y * (threehalfs - (halfX * y * y));
        return y;
    }
    inline double __fastcall rsqrt(double x) {
        const long long magic = 0x5FE6EB50C7B537AA;
        const double threehalfs = 1.5, halfX = x * 0.5;
        unsigned long long bits = *(unsigned long long*) & x;
        double y;

        bits = magic - (bits >> 1);
        y = *(float*)&bits;
        y = y * (threehalfs - (halfX * y * y));
        y = y * (threehalfs - (halfX * y * y));
        return y;
    }
    inline float __fastcall approximateLog2(float x) {
        unsigned int bits = *(unsigned int*)&x; // get bits
        unsigned char E = bits >> 23; // get E
        bits = ((bits << 9) >> 9) + 0x3F800000; // get 1 + m
        float t = *(float*)&bits; // put 1+m in float
        return 2.0f * sqrt(t) + E - 129; // \log_2x\approx2(\sqrt{1+m}-1)+E-B=2\sqrt{1+m}+E-B-2

        // no need to do
        // \log_2x\approx4(\sqrt{\sqrt{1+m}}-1)+E-B=4\sqrt{\sqrt{1+m}}+E-B-4
        // return 4.0f * fastSqrt(fastSqrt(t)) + E - 131;
    }
    inline double __fastcall approximateLog2(double x) {
        unsigned long long bits = *(unsigned long long*) & x;
        unsigned short E = bits >> 52;
        bits = ((bits << 12) >> 12) + 0x3FF0000000000000;
        double t = *(double*)&bits;
        return 2.0 * sqrt(t) + E - 1025;
    }
}

最快实现

最坏的情况下花费时间和上面差不多

// intel语法实现（Clang、ICC）
// GNU要使用AT&T语法，我不是很熟悉AT&T语法

namespace math {
    // C++哪有汇编快啊[doge]
    inline float __fastcall regSqrt(float x) {
    	__asm {
        	vsqrtss xmm0, xmm0, xmm0
        	vmovss x, xmm0
    	}
    	return x;
	}
    inline double __fastcall regSqrt(double x) {
    	__asm {
        	vsqrtsd xmm0, xmm0, xmm0
        	vmovsd x, xmm0
    	}
    	return x;
	}
    inline float regRSqrt(float x) {
        // SSE, AVX
        __asm {
            rsqrtss xmm0, xmm0
            vmovss x, xmm0
        }
        return x;
    }
    // 没有双精度浮点平方根的寄存器实现
}

不知道能不能直接使用内联汇编进行计算，感觉会有寄存器相关的竞争，保守一点还是调用函数）