平方根倒数快速算法
平方根常出现在游戏的图形计算中，尤其是求一个向量的基向量时
约翰卡马克的代码
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;
}