平方根倒数速算法

近似求解平方根倒数
最近更新:2019-08-14

【参考维基百科】平方根倒数速算法

源码

1
2
3
4
5
6
7
8
9
10
float invSqrt(float x) {
float halfx = 0.5f * x;
float y = x;
long i = *(long*)&y;
i = 0x5f3759df - (i>>1);
y = *(float*)&i;
y = y * (1.5f - (halfx * y * y));
y = y * (1.5f - (halfx * y * y));
return y;
}

该代码是先求出平方根倒数的近似值,然后使用牛顿迭代法进行精确。


浮点数

浮点数的存储格式

  • 首位符号位
  • 8 位表示经过偏移处理后的指数,表示范围为 $[-127, \, 128]$
  • 最后 23 位表示有效数字中除最高位外的其余数字
  • m 表示有效数字的尾数($0 \le m \lt 1$)
  • 偏移量 $B = 127$
  • 以上图为例
    • $m = 2^{-2} = 0.25$
    • $E - B = 4 + 8 + 16 + 32 + 64 - 127 = -3$
    • $x = (1 +0.25) \cdot 2^{-3} = 0.15625$


将浮点数取别名存储为整数

该整数的数值为

  • $E$ 表示指数
  • $M$ 表示有效数字
  • 以上图为例
    • $E = 124$,$M = 2^{21}$
    • $I = 124 \times 2^{23} + 2^{21}$

魔术数字

Sign(符号) Exponent(指数) Mantissa(尾数)
1 位 $b$ 位 $(n - 1 -b)$ 位
  • 一共 $n$ 位
  • 浮点数表示实际数值时
    • $x = (1 + m_x) \cdot 2^{e_x}$
  • 若为整数表示时
    • $I_x = E_x L + M_x$
    • $L = 2^{n - 1 -b}$
  • 导出
    • $m_x = \dfrac{M_x}{L}$
    • $e_x = E_x - B$,其中 $B = 2^{b-1} - 1$

对等式的两边取二进制对数有

  • $\log_2 (y) = - \dfrac{1}{2} \log_2 (x)$
  • 带入 $x = (1 + m_x) \cdot 2^{e_x}$
    • $\log_2 (1 + m_y) + e_y = - \dfrac{1}{2} \log_2 (1 + m_x) - \dfrac{1}{2} e_x$
    • 因为 $0 \lt m_y \lt 1$,$0 \lt m_x \lt 1$
      • 且 $0 \lt x \lt 1$ 时,$\log_2(1 + x) \approx x$
      • 设 $0 \le \sigma \le \dfrac{1}{3}$,$\log_2(1 + x) \approx x +\sigma$
      • 当近似值 $R = {\rm 0x5f3759df}$ 时,$\sigma \approx 0.04504618758$
    • 将 $\log_2(1 + x) \approx x +\sigma$ 代入上式
      • $m_y + \sigma + e_y = - \dfrac{1}{2} m_x - \dfrac{1}{2} \sigma - \dfrac{1}{2} e_x$
    • 带入 $M_x$,$E_x$,$B$ 和 $L$
      • $M_x + (E_y - B) L = - \dfrac{1}{2} M_x - \dfrac{1}{2} (E_x - B) L - \dfrac{3}{2} \sigma L$
      • $E_y L +M_y = \dfrac{3}{2} (B - \sigma) L - \dfrac{1}{2} (E_x L + M_x)$
      • $\begin{array}{l} I_y &= \dfrac{3}{2} (B - \sigma) L - \dfrac{1}{2} I_x \\ &= R - \dfrac{1}{2} I_x \end{array}$
1
R = 0x5f3759df = 0 10111110 01101110101100111011111 = 1597463007;

$\dfrac{3}{2} (B - \sigma) L = \dfrac{3}{2} B L - \dfrac{3}{2} \sigma L$

  • $\begin{array}{l} \dfrac{3}{2} (B - \sigma) L &= \dfrac{3}{2} B L - \dfrac{3}{2} \sigma L \\\\ &= \dfrac{3}{2} ( 2^{b-1} - 1 ) 2^{n-1-b} - \dfrac{3}{2} \sigma L \\\\ &= 3 \cdot 2^{n-3} - \dfrac{3}{2} L - \dfrac{3}{2} \sigma L \\\\ &= 3 \cdot 2^{n-3} - \dfrac{3}{2} (1 + \sigma) L \\\\ \end{array}$ 且 $n = 32$,$b = 8$

    • $3 \cdot 2^{n-3} = 3 \cdot 2^{29} = 1610612736$
    • $L = 2^{b-1} -1 = 127$
      • $190.5 \le \dfrac{3}{2} (1 + \sigma) L \le 254$
    • $1610612545.5 \le 3 \cdot 2^{n-3} - \dfrac{3}{2} (1 + \sigma) L \le 1610612482$
  • 所以 $\dfrac{3}{2} (B - \sigma) L \approx R$


牛顿迭代法

1
近似求解方程的方法
  • 假设求解方程 $f(x) = 0$,且解为 $x = x_0$
  • 在 $x_0$ 旁取一近似点 $x_1$
    • 过 $(x_1, \, f(x_1))$ 作 $f(x)$ 的切线,切线方程为 $y = f(x_1) + f’(x_1)(x - x_1)$
    • 切线与 x 轴的交点为 $x_2 = x_1 - \dfrac{ f(x_1) }{ f’(x_1) }$
  • 将 $x_2$ 作为新的近似点进行迭代,直到 $x_{n+1} - x_n < \text{阈值}$

案例

设 $f(x) = \dfrac{1}{ \sqrt{x} }$,求解 $x = x’$ 时,$f(x’) \approx ?$

  • 即求解 $\dfrac{1}{y^2} - x = 0$ 的解,即 $f(y) = \dfrac{1}{y^2} - x = 0$
  • 设求出某一近似点为 $y_n$,则更精确的值为 $y_{n+1} = y_n - \dfrac{ f(y_n) }{ f’(y_n) }$
    • $f(y_n) = \dfrac{1}{y_n^2} - x’$,$f’(y_{n+1}) = - \dfrac{2}{y_n^3}$
  • 所以 $y_{n+1} = \dfrac{ y_n (3 - x’ y_n^2) }{2} = y_n \left(1.5 - \dfrac{1}{2} x’ y_n^2 \right)$

true