日學一算法---平方根倒數速算法【轉】

本文來源：http://blog.csdn.net/xiaoguohaha/article/details/21652643

文章太牛了，裏面的算法更是讓我膜拜萬分，不得不轉過來收藏一下。

 

 

---

 

     咱們平時常常會有一些數據運算的操做，須要調用sqrt，exp，abs等函數，那麼時候你有沒有想過：這個些函數系統是如何實現的？就拿最經常使用的sqrt函數來講吧，系統怎麼來實現這個常常調用的函數呢？

     雖然有可能你平時沒有想過這個問題，不過正所謂是「臨陣磨槍，不快也光」，你「眉頭一皺，計上心來」，這個不是太簡單了嘛，用二分的方法，在一個區間中，每次拿中間數的平方來試驗，若是大了，就再試左區間的中間數；若是小了，就再拿右區間的中間數來試。好比求sqrt(16)的結果，你先試（0+16）/2=8，8*8=64，64比16大，而後就向左移，試（0+8）/2=4，4*4=16恰好，你獲得了正確的結果sqrt(16)=4。而後你三下五除二就把程序寫出來了：

 1 //用二分法
 2 float SqrtByBisection(float n)
 3 {
 4     //小於0的按照你須要的處理
 5     if(n < 0)
 6         return n;
 7     float mid,last;
 8     float low,up;
 9     low=0,up=n;
10     mid=(low+up)/2;
11     do
12     {
13         if(mid*mid>n)
14             up=mid;
15         else
16             low=mid;
17         last=mid;
18         mid=(up+low)/2;
19     }
20     //精度控制
21     while(abs(mid-last) > eps);
22     return mid;
23 }

View Code

     而後看看和系統函數性能和精度的差異（其中時間單位不是秒也不是毫秒，而是CPU Tick，無論單位是什麼，統一了就有可比性）。二分法和系統的方法結果上徹底相同，可是性能上整整差了幾百倍。爲何會有這麼大的區別呢？難道系統有什麼更好的辦法？難道。。。。哦，對了，回憶下咱們曾經的高數課，曾經老師教過咱們「牛頓迭代法快速尋找平方根」，或者這種方法能夠幫助咱們，具體步驟以下。

求出根號a的近似值：首先隨便猜一個近似值x，而後不斷令x等於x和a/x的平均數，迭代個六七次後x的值就已經至關精確了。例如，我想求根號2等於多少。假如我猜想的結果爲4，雖然錯的離譜，但你能夠看到使用牛頓迭代法後這個值很快就趨近於根號2了：

(       4  + 2/4        ) / 2 = 2.25

(     2.25 + 2/2.25     ) / 2 = 1.56944..

( 1.56944..+ 2/1.56944..) / 2 = 1.42189..

( 1.42189..+ 2/1.42189..) / 2 = 1.41423..

....

這種算法的原理很簡單，咱們僅僅是不斷用(x,f(x))的切線來逼近方程x^2-a=0的根。根號a實際上就是x^2-a=0的一個正實根，這個函數的導數是2x。也就是說，函數上任一點(x,f(x))處的切線斜率是2x。那麼，x-f(x)/(2x)就是一個比x更接近的近似值。代入 f(x)=x^2-a獲得x-(x^2-a)/(2x)，也就是(x+a/x)/2。

相關的代碼以下： 

 1 float SqrtByNewton(float x)
 2 {
 3     // 最終
 4     float val = x;
 5     // 保存上一個計算的值
 6     float last;
 7     do
 8     {
 9         last = val;
10         val =(val + x/val) / 2;
11     }
12     while(abs(val-last) > eps);
13     return val;
14 }

View Code

     牛頓迭代法性能提升了不少，但是和系統函數相比，仍是有這麼大差距，這是爲何呀？想啊想啊，想了好久仍然百思不得其解。忽然有一天，我在網上看到一個神奇的方法，因而就有了今天的這篇文章，廢話很少說，看代碼先：

 1 float InvSqrt(float x)
 2 {
 3     float xhalf = 0.5f*x;
 4     int i = *(int*)&x; // get bits for floating VALUE
 5     i = 0x5f375a86- (i>>1); // gives initial guess y0
 6     x = *(float*)&i; // convert bits BACK to float
 7     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
 8     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
 9     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
10 
11     return 1/x;
12 }

View Code

此次真的是質變了，結果居然比系統的還要好。到如今你是否是還不明白那個「鬼函數」，到底爲何速度那麼快嗎？不急，先看看下面的故事吧：

Quake-III Arena (雷神之錘3)是90年代的經典遊戲之一。該系列的遊戲不但畫面和內容不錯，並且即便計算機配置低，也能極其流暢地運行。這要歸功於它3D引擎的開發者約翰-卡馬克（John Carmack）。事實上早在90年代初DOS時代，只要能在PC上搞個小動畫都能讓人驚歎一番的時候，John Carmack就推出了石破天驚的Castle Wolfstein, 而後再接再勵，doom, doomII, Quake...每次都把3-D技術推到極致。他的3D引擎代碼資極度高效，幾乎是在壓榨PC機的每條運算指令。當初MS的Direct3D也得聽取他的意見，修改了很多API。

最近，QUAKE的開發商ID SOFTWARE 遵照GPL協議，公開了QUAKE-III的原代碼，讓世人有幸目擊Carmack傳奇的3D引擎的原碼。這是QUAKE-III原代碼的下載地址： http://www.fileshack.com/file.x?fid=7547。咱們知道，越底層的函數，調用越頻繁。3D引擎歸根到底仍是數學運算。那麼找到最底層的數學運算函數（在game/code/q_math.c）， 必然是精心編寫的。裏面有不少有趣的函數，不少都使人驚奇，估計咱們幾年時間都學不完。在game/code/q_math.c裏發現了這樣一段代碼。它的做用是將一個數開平方並取倒，經測試這段代碼比(float)(1.0/sqrt(x))快4倍：

 1 float Q_rsqrt( float number )
 2 {
 3     long i;
 4     float x2, y;
 5     const float threehalfs = 1.5F;
 6 
 7     x2 = number * 0.5F;
 8     y   = number;
 9     i   = * ( long * ) &y;   // evil floating point bit level hacking
10     i   = 0x5f3759df - ( i >> 1 ); // what the fuck?
11     y   = * ( float * ) &i;
12     y   = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
13     // y   = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
14 
15     #ifndef Q3_VM
16     #ifdef __linux__
17          assert( !isnan(y) ); // bk010122 - FPE?
18     #endif
19     #endif
20     return y;
21 } 

View Code

函數返回1/sqrt(x)，這個函數在圖像處理中比sqrt(x)更有用。注意到這個函數只用了一次疊代！（其實就是根本沒用疊代，直接運算）。編譯，實驗，這個函數不只工做的很好，並且比標準的sqrt()函數快4倍！要知道，編譯器自帶的函數，但是通過嚴格仔細的彙編優化的啊！

這個簡潔的函數，最核心，也是最讓人費解的，就是標註了「what the fuck?」的一句：i = 0x5f3759df - ( i >> 1 );

再加上y = y * ( threehalfs - ( x2 * y * y ) );

兩句話就完成了開方運算！並且注意到，核心那句是定點移位運算，速度極快！特別在不少沒有乘法指令的RISC結構CPU上，這樣作是極其高效的。

算法的原理其實不復雜,就是牛頓迭代法,用x-f(x)/f'(x)來不斷的逼近f(x)=a的根。

沒錯，通常的求平方根都是這麼循環迭代算的可是卡馬克(quake3做者)真正牛B的地方是他選擇了一個神祕的常數0x5f3759df 來計算那個猜想值，就是咱們加註釋的那一行，那一行算出的值很是接近1/sqrt(n)，這樣咱們只須要2次牛頓迭代就能夠達到咱們所須要的精度。好吧若是這個還不算NB，接着看：

普渡大學的數學家Chris Lomont看了之後以爲有趣，決定要研究一下卡馬克弄出來的這個猜想值有什麼奧祕。Lomont也是個牛人，在精心研究以後從理論上也推導出一個最佳猜想值，和卡馬克的數字很是接近, 0x5f37642f。卡馬克真牛，他是外星人嗎？

傳奇並無在這裏結束。Lomont計算出結果之後很是滿意，因而拿本身計算出的起始值和卡馬克的神祕數字作比賽，看看誰的數字可以更快更精確的求得平方根。結果是卡馬克贏了... 誰也不知道卡馬克是怎麼找到這個數字的。

最後Lomont怒了，採用暴力方法一個數字一個數字試過來，終於找到一個比卡馬克數字要好上那麼一丁點的數字，雖然實際上這兩個數字所產生的結果很是近似，這個暴力得出的數字是0x5f375a86。

Lomont爲此寫下一篇論文，"Fast Inverse Square Root"。 論文下載地址：http://www.math.purdue.edu/~clomont/Math/Papers/2003/InvSqrt.pdf ，http://www.matrix67.com/data/InvSqrt.pdf。

最後，給出最精簡的1/sqrt()函數：

1 float InvSqrt(float x)
2 {
3     float xhalf = 0.5f*x;
4     int i = *(int*)&x; // get bits for floating VALUE
5     i = 0x5f375a86- (i>>1); // gives initial guess y0
6     x = *(float*)&i; // convert bits BACK to float
7     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
8     return x;
9 }

View Code

你們能夠嘗試在PC機、5一、AVR、430、ARM、上面編譯並實驗，驚訝一下它的工做效率。

前兩天有一則新聞，大意是說 Ryszard Sommefeldt 好久之前看到這麼樣的一段 code (可能出自 Quake III 的 source code)：

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

View Code

他一看之下驚爲天人，想要拜見這位前輩高人，可是一路追尋下去卻一直找不到人；同時間也有其餘人在找，雖然也沒找到出處，可是 Chris Lomont 寫了一篇論文 (in PDF) 解析這段 code 的算法 (用的是 Newton’s Method，牛頓法；比較重要的是後半段講到怎麼找出神奇的 0x5f3759df 的)。

PS. 這個 function 之因此重要，是由於求 開根號倒數 這個動做在 3D 運算 (向量運算的部份) 裏面經常會用到，若是你用最原始的 sqrt() 而後再倒數的話，速度比上面的這個版本大概慢了四倍吧… XD

PS2. 在他們追尋的過程當中，有人提到一份叫作 MIT HACKMEM 的文件，這是 1970 年代的 MIT 強者們作的一些筆記 (hack memo)，大部份是 algorithm，有些 code 是 PDP-10 asm 寫的，另外有少數是 C code (有人整理了一份列表)。

好了，故事就到這裏結束了，但願你們能有有收穫:)

原文轉自：http://www.nowamagic.net/algorithm/algorithm_EfficacyOfFunctionSqrt.php

 

說實話，通過測試，本文所提到的計算開平方的算法在精度方面跟sqrt相比仍是有差距的。

測試時，我把精度定爲1e-7，分別用兩個算法計算了1！10000的開平方根，發現不少數

經二者開根號結果是不同的。

下面是測試代碼：

 

 1 #include<stdio.h>
 2 #include<math.h>
 3 double delta=1e-7;
 4 float InvSqrt(float x);
 5 int comparDoubleNum(double x,double y);
 6 int main()
 7 {
 8     int i;
 9     double ans1,ans2;
10     freopen("out.txt","w",stdout);
11     for(i=1;i<=10000;i++)
12     {
13         ans1=sqrt(i);
14         ans2=InvSqrt(i);
15         printf("%d---%lf---%lf---  %d\n",i,ans1,ans2,comparDoubleNum(ans1,ans2));
16     }
17     return 0;
18 }
19 float InvSqrt(float x)
20 {
21     float xhalf = 0.5f*x;
22     int i = *(int*)&x; // get bits for floating VALUE
23     i = 0x5f375a86- (i>>1); // gives initial guess y0
24     x = *(float*)&i; // convert bits BACK to float
25     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
26     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
27     x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
28 
29     return 1/x;
30 }
31 int comparDoubleNum(double x,double y)//當x與y之差的絕對值小於0.00001（即：1e-5）時 認爲x等於y 
32 {
33     if((x-y)>delta) return 1;
34     else if((x-y)<-delta) return -1;
35     else return 0;
36 }

View Code