[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))
}

牛頓-拉弗森方法
$$ 0=(x-x_0) \cdot {f}'(x_0)+f(x_0) \\ x_{n+1}=x_n-\frac{f(x_n)}{{f}'(x_n)} $$

證明

在單精度浮點數的定義下
$$ x=(-1)^s\cdot 2^{{e}'-127}\cdot (1+m)=(-1)^s\cdot 2^e\cdot (1+m) $$ 開根號，故只看正浮點數，\(s=0\)，且將對應的 4 bytes 轉成整數來看的話
$$ \begin{align*} X &= EL+M \\ E &= e+127=E+B \\ M &= m\cdot 2^{23}=mL \\ \end{align*} $$ $$ \begin{align*} log_2\ x&=log_2\left [2^e(1+m) \right ]\\ &=e+log_2(1+m)\\ &\approx e+m+\delta \end{align*} $$ 最大 \(\delta=0.04504656791687011719\)
$$ \begin{align*} log_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{align*} $$ 可得
$$ log_2\ x\approx \frac{X}{L}-(B-\delta)\cdots (1) \\ $$ 當
$$ {x}'=\frac{1}{\sqrt{x}} \Rightarrow log_2\ {x}'=-\frac{1}{2}log_2\ x $$ 代入\((1)\)
$$ \begin{align*} \frac{{X}'}{L}-(B-\delta)&=-\frac{1}{2}\left (\frac{X}{L}-(B-\delta) \right )\\ \frac{{X}'}{L} &= -\frac{1}{2}\frac{X}{L}+\frac{3}{2}(B-\delta)\\ {X}' &= -\frac{1}{2}X+\frac{3}{2}L(B-\delta)\\ {X}' &= \frac{3}{2}L(B-\delta)-\frac{1}{2}X\\ \end{align*} $$ 將 \(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
如何通俗易懂地讲解牛顿迭代法？