1/sqrt(x)的速算法(InvSqrt)

在網路上看到一篇

開根號倒數 (InvSqrt(), 1 / sqrt(x)) 速算法

是求開根號倒數的速數法，

雖然有誤差但是出奇的快，

比1/sqrt(x)快了約4倍，

Code 如下：

float InvSqrt (float x) {
	float xhalf = 0.5f*x;
	int i = *(int*)&x;
	i = 0x5f3759df - (i>>1);
	x = *(float*)&i;
	x = x*(1.5f - xhalf*x*x);
	return x;
}

看了真是驚嘆不已，

一時也看不懂為什麼…

好在 Chris Lomont 寫了一篇論文 (in PDF) 來解說，

這段code是由牛頓法而來，

複習一下牛頓法即Xn+1=Xn-(f(Xn)/f'(Xn))

此Xn+1將會比Xn更趨近f(x)=0之x值。

Let y=1/sqrt(x), f(y)=1/y^2-x,

(PS. 此邊和之後的 sqrt(x)非指sqrt function call，

指的是x開根號的數學意義)

give better one Yn+1 using

Yn+1=Yn-(f(Yn)/f'(Yn))=Yn-((1/(Yn)^2-x)/(-2)*(Yn)^(-3))

=Yn-(1/2)(-Yn+x(Yn)^3)=(3/2)Yn-(1/2)x(Yn)^3

=((1/2)Yn)(3-x(Yn)^2)

也就是原code中的 x = x*(1.5f - xhalf*x*x);

之後把float number 視為 0,E,M三個欄位，

此正是電腦儲存float number的格式，

由y=sqrt(x)=(1/sqrt(1+M))2^(-e/2),

Chris Lomont的論文說

For the general case we take a magic constant R,
and analyze y0 = R−(i>>1),

至於怎麼analyze的，到此好像沒說很清楚，

不過可能是把y=1/sqrt(x)，

找一條漸近線，取其導數得到-1/2為斜率的直線，

視做(-1/2)x+a，也就相當於(-1/2)i+R，

即y0 = R−(i>>1),

論文之後有一大串都在說R怎麼來的，

並且有計算其relative error，

分成E欄位為even和odd的情況討論之，

並實測all float number得到relative error最大值仍然很小，

所以此法不只神奇，最重要的是"有用的"。

最後論文下結論時，

找來三個magic number來做比較，第一個是code中的magic number，

第二個是前半論文算出的magic number，還有第三個magic number。

論文好像是說，第二個number沒考慮到牛頓法步驟的影響，

所以經過牛頓法一次後，並不能得到最小的maximal relative error，

考慮牛頓法後，去實測才得到第三個magic number，

也就是0x5f375a86，才是最佳的。

原論文摘錄如下：

The results of the theory
implied any good approximations should be near those few above, thus
limiting the range of the search.
Starting at the initial constant, and testing all constants above and
below until the maximal relative error exceeds 0.00176 gives the third
constant 0x5f375a86 as the best one

這邊讓人覺得有點雞助…

原論文也提到，

It is clear that y0 is a nonlinear approximation; however, by
being nonlinear, it actually fits better!

注意此approximation是 nonlinear的，

因此讓我感到極為神奇。

以上是小弟個人不才的一些見解，

如有錯誤，還請指教。