[Go] Fast Inverse Square Root

程式語言：Go

Package：

math

功能：快速計算出平方根的倒數

程式碼

Playground

package main
 
import (
    "fmt"
    "math"
)
 
func InvSqrt(x float32) float32 {
    var xhalf float32
    xhalf = 0.5 * x
    i := math.Float32bits(x)
    i = 0x5f3759df - (i >> 1) // 證明如下，直接求得接近的值
    x = math.Float32frombits(i)
    x = x * (1.5 - xhalf*x*x) // 牛頓-拉弗森方法
 
    return x
}
 
func main() {
    fmt.Println(InvSqrt(4))
}

牛頓-拉弗森方法
$$\begin{gathered} 0=(x-x_0)\cdot f^{\prime }(x_0)+f(x_0) \\ {x}_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)} \end{gathered}$$

證明

在單精度浮點數的定義下
$$x=(-1{)}^s\cdot 2^{e^{\prime }-127}\cdot (1+m)=(-1{)}^s\cdot 2^e\cdot (1+m)$$ 開根號，故只看正浮點數，$s=0$，且將對應的 4 bytes 轉成整數來看的話
$$\begin{array}{rl}X& =EL+M\\ E& =e+127=E+B\\ M& =m\cdot 2^{23}=mL\end{array}$$ $$\begin{array}{rl}lo{g}_{2}x& =lo{g}_2[2^e(1+m)]\\ & =e+lo{g}_2(1+m)\\ & \approx e+m+\delta \end{array}$$ 最大 $\delta =0.04504656791687011719$
$$\begin{array}{rl}lo{g}_{2}x& \approx e+m+\delta \\ \because e=E-B,m=\frac{M}{L}\therefore & =E-B+\frac{M}{L}+\delta \\ & =\frac{EL+M}{L}-(B-\delta )\\ \because X=EL+M\therefore & =\frac{X}{L}-(B-\delta )\end{array}$$ 可得
$$lo{g}_{2}x\approx \frac{X}{L}-(B-\delta )\cdots (1)$$ 當
$$x^{\prime }=\frac{1}{\sqrt{x}}\Rightarrow lo{g}_2{x}^{\prime }=-\frac{1}{2}lo{g}_{2}x$$ 代入$(1)$
$$\begin{array}{rl}\frac{X^{\prime }}{L}-(B-\delta )& =-\frac{1}{2}(\frac{X}{L}-(B-\delta ))\\ \frac{X^{\prime }}{L}& =-\frac{1}{2}\frac{X}{L}+\frac{3}{2}(B-\delta )\\ X^{\prime }& =-\frac{1}{2}X+\frac{3}{2}L(B-\delta )\\ X^{\prime }& =\frac{3}{2}L(B-\delta )-\frac{1}{2}X\end{array}$$ 將 $L=2^{23},B=127,\delta =0.04504656791687011719$ 的值代入，可得證

i = 0x5f3759df - (i >> 1) 

依此想法，可推廣至其他的數學運算，快速求一解
像是立方根、平方根、倒數立方根等…
但 $L=2^{23},B=127,s=0$ 可能需再調整，此處的條件為 32-bit 的單精度浮點數，求平方根倒數

參考

0x5f3759df 这个快速开方中的常数的数学依据和原理
FAST INVERSE SQUARE ROOT
如何通俗易懂地讲解牛顿迭代法？