快速求平方根

1. 前置知识

• 建议首先阅读我的另外一篇文章《雷神之锤 III 竞技场》快速求平方根倒数的计算探究》。
• 建议大家自己看过《雷神之锤 III 竞技场》快速求平方根倒数的计算探究》学会快速求平方根倒数算法后，不看我这篇文章，自己推导一篇快速求平方根的算法，这也是对自己是否彻底掌握的一个测试。
• 本篇文章是《雷神之锤 III 竞技场》快速求平方根倒数的计算探究》基础上的继续推导，写出快速计算平方根的函数，要计算的公式如下：
  $$\begin{gathered} y=\sqrt{x} \\ y=x^{\frac{1}{2}} \end{gathered}$$

2. 计算过程

2.1第一步：对等式两边同时取以2为底的对数

$$\begin{gathered} {\log }_2(y)={\log }_2(x^{\frac{1}{2}}) \\ {\log }_2(y)=\frac{1}{2}{\log }_2(x) \end{gathered}$$

2.2 第二步：浮点数取log的整数表达形式

《IEEE754 -浮点数的表示》 提到浮点数的计算公式为：$$(-1{)}^S\ast (1+\frac{F}{2^{23}})\ast 2^{(E-127)}$$

• 不考虑符号位。对其取对数
  $${\log }_2\left(\left(1+\frac{\mathrm{F}}{2^{23}}\right)\ast 2^{\mathrm{E}-127}\right)$$
  $$\approx \frac{1}{2^{23}}(F+2^{23}\ast E)+u-127$$

2.3第三步：将 第二步 带入 第一步

• 化简一下公式
  $$\begin{gathered} lo{g}_2(y)=\frac{1}{2}{\log }_2(x) \\ \frac{1}{2^{23}}({\mathrm{F}}_y+2^{23}\ast {\mathrm{E}}_y)+\mu -127=\frac{1}{2}\left(\frac{1}{2^{23}}({\mathrm{F}}_x+2^{23}\ast {\mathrm{E}}_x)+\mu -127\right) \\ ------ \\ ({\mathrm{F}}_y+2^{23}\ast {\mathrm{E}}_y)=\frac{1}{2}\ast 2^{23}(127-\mu )+\frac{1}{2}({\mathrm{F}}_x+2^{23}\ast {\mathrm{E}}_x) \end{gathered}$$
• 结果已经呼之欲出了，如下图所示,红框代表的是$y$的二进制值， 蓝框代表的是$x$的二进制值。
  

$$y_{二进制}=\frac{1}{2}\ast 2^{23}(127-\mu )+\frac{1}{2}(x_{二进制})$$

• 令$\mu =0.0450465679168701171875$(这个数字是《雷神之锤 III 竞技场》快速求平方根倒数的计算探究》中所给出的魔数反推出来的）,计算出 新的魔数 。
  $$y_{二进制}=0x1fbd1df5+\frac{1}{2}(x_{二进制})$$

3. 最后在使用牛顿迭代法再次逼近一次结果，提高精度

3.1 将$y=\sqrt{x}$变换为$f(x)=0$的形式

• 消除根号$$\begin{gathered} y=\sqrt{x} \\ {y}^2=x \end{gathered}$$
• 整理一下,得到$f(x)=0$ 的形式：$$y^2-x=0$$
• 带入牛顿迭代法的通用公式$$\begin{gathered} f(y)=y^2-x \\ {f}^{\prime }(y)=2y \\ ----------- \\ {y}_{n+1}=y_n-\frac{f(y_n)}{f^{\prime }(y_n)} \\ ----------- \\ {y}_{n+1}=y_n-\frac{y_n^2-x}{2y_n} \\ ----------------------- \\ {y}_{n+1}=\frac{1}{2}\left(y_n+\frac{x}{y_n}\right) \end{gathered}$$

4. 代码实现

float Q_sqrt(float number) {
    int i;
    float x, y;
    const float half = 0.5f;
    x = number;
    
    // 1. 位级黑魔法：生成初始近似值
    i = *(long*)&x;           // 将浮点位模式转为整数
    i = 0x1fbd1df5 + (i >> 1); // 魔法常数调整
    y = *(float*)&i;         // 转回浮点数
    
    // 2. 牛顿迭代法优化精度
    y = half * (y + x / y);  // 第一次迭代
    // y = half * (y + x / y);  // 可选：第二次迭代
    
    return y;
}