平方根倒数快速算法

平方根倒数快速算法

平方根常出现在游戏的图形计算中，尤其是求一个向量的基向量时

约翰卡马克的代码

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）

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

浮点数

有效数字是从小数点后第一位开始的，小数点前的数字总为1，即

\[对任意浮点数x，总有x=(-1)^{Signature}\cdot(1+m)\cdot2^{E-B}\\ 其中const\ B = 127 \]

上图中可以得知符号部分0指数01111100数字部分1.01000000000000000000000，即0.15625

求平方根转化为求对数

\[y = Sqrt(x)，即y = x^{\frac{1}{2}}\\ \log_2y = \log_2x^{\frac{1}{2}}=\frac{1}{2}\log_2x\\ 令z=\log_2y\\ z=\frac{1}{2}\log_2x\\ \\ 先尝试对\log_2x进行变换\\ 把x=(-1)^{Signature}\cdot(1+m)\cdot2^{E-B}且Signature=0带入\log_2x，得\\ \log_2x=\log_2[(1+m)\cdot2^{E-B}]\\ 用m_{bit}表示m的二进制\\ _{（这里由于m表示的是一个小数，而m的二进制本身可以被视为整数，所以需要区分开来）}\\ \log_2x=\log_2[(1+\frac{m_{bit}}{2^{23}})\cdot2^{E-B}]\\ \log_2x=\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_2x\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_2x=\frac{1}{2^{23}}x_{bit}-B\qquad\qquad(1)\\ \\ 把(1)式代回z=\frac{1}{2}\cdot\log_2x，得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\\ \]

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

直接求平方根倒数

运算除法的开销比运算乘法高，能否直接求平方根倒数后再进行乘法呢？

\[y=x^{-\frac{1}{2}}\\ \log_2y = -\frac{1}{2}\log_2x\\ 把(1)式代入，得\\ \frac{1}{2^{23}}\cdot y_{bit}-B=-\frac{1}{2}\cdot(\frac{1}{2^{23}}\cdot x_{bit}-B)\\ 整理，得y_{bit}=3B\cdot2^{22}-\frac{1}{2}x_{bit}\\ 到这里就是最终形式了 \]

代入数据

float

\(y_{bit}=381\cdot2^{22}-\frac{1}{2}x_{bit}\)

对应的值为0x5F40 0000

double

推广到双精度浮点型变量

\(y_{bit}=3069\cdot2^{51}-\frac{1}{2}x_{bit}\)

对应的值为0x5FE8 0000 0000 0000

调整

为什么上述计算得到的值为0x5F40 0000，而不是0x5F37 59DF呢

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

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

后续处理

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

牛顿迭代法

原理

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

平方根

表达式已经在上述的计算中得到了：\(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的差值，得到0x1FB759DF

代码

constexpr float INV_MAGIC = 0x5F3759DF;
constexpr float MAGIC = 0x1FB759DF;

constexpr double INV_DMAGIC = 0x5FE6800000000000;
constexpr double DMAGIC = 0x1FF6800000000000;

float FastInvSqrt(float input)
{
    // y = input^{-\frac{1}{2}}
    const float threehalfs = 1.5F, halfInput = 0.5 * input;
    int32_t i = *(int32_t*)&input;
    float y;
    
    i = INV_MAGIC - (i >> 1);
    y = *(float*)&i;
    y = y * (threehalfs - halfInput * y * y);
    return y;
}
double FastInvSqrt(double input)
{
    // y = input^{-\frac{1}{2}}
    const double threehalfs = 1.5, halfInput = 0.5 * input;
    int64_t i = *(int64_t*)&input;
    double y;
    
    i = INV_DMAGIC - (i >> 1);
    y = *(double*)&i;
    y = y * (threehalfs - halfInput * y * y);
    y = y * (threehalfs - halfInput * y * y);
    return y;
}

float FastSqrt(float input)
{
    // y = input^{\frac{1}{2}}
    int32_t i = *(int32_t*)&input;
    float y;

    i = MAGIC + (i >> 1);
    y = *(float*)&i;
    y = 0.5 * (y + input / y);
    y = 0.5 * (y + input / y);
    return y;
}
double FastSqrt(double input)
{
    // y = input^{\frac{1}{2}}
    int64_t i = *(int64_t*)&input;
    double y;
    
    i = DMAGIC + (i >> 1);
    y = *(double*)&i;
    y = 0.5 * (y + input / y);
    y = 0.5 * (y + input / y);
    y = 0.5 * (y + input / y);
    return y;
}