我們平時經常會有一些數據運算的操作，需要調用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。然後你三下五除二就把程序寫出來了：
float SqrtByBisection(float n) //用二分法
{
if(n<0) //小於0的按照你需要的處理
return n;
float mid,last;
float low,up;
low=0,up=n;
mid=(low+up)/2;
do
{
if(mid*mid>n)
up=mid;
else
low=mid;
last=mid;
mid=(up+low)/2;
}while(abs(mid-last) > eps);//精度控制
return mid;
} 然後看看和系統函數性能和精度的差別（其中時間單位不是秒也不是毫秒，而是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。
相關的代碼如下：
float SqrtByNewton(float x)
{
float val = x;//最終
float last;//保存上一個計算的值
do
{
last = val;
val =(val + x/val) / 2;
}while(abs(val-last) > eps);
return val;
}然後我們再來看下性能測試：
哇塞，性能提高了很多，可是和系統函數相比，還是有這麼大差距，這是為什麼呀？想啊想啊，想了很久仍然百思不得其解。突然有一天，我在網上看到一個神奇的方法，於是就有了今天的這篇文章，廢話不多說，看代碼先：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return 1/x;
}然後我們最後一次來看下性能測試：
這次真的是質變了，結果竟然比系統的還要好。。。哥真的是震驚了！！！哥吐血了！！！一個函數引發了血案！！！血案，血案。。。
到現在你是不是還不明白那個“鬼函數”，到底為什麼速度那麼快嗎？不急，先看看下面的故事吧：
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
(下面是官方的下載網址，搜索 “quake3-1.32b-source.zip” 可以找到一大堆中文網頁的。ftp://ftp.idsoftware.com/idstuff/source/quake3-1.32b-source.zip)
我們知道，越底層的函數，調用越頻繁。3D引擎歸根到底還是數學運算。那麼找到最底層的數學運算函數（在game/code/q_math.c）， 必然是精心編寫的。裡面有很多有趣的函數，很多都令人驚奇，估計我們幾年時間都學不完。在game/code/q_math.c裡發現了這樣一段代碼。它的作用是將一個數開平方並取倒，經測試這段代碼比(float)(1.0/sqrt(x))快4倍：
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the fuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
#ifndef Q3_VM
#ifdef __linux__
assert( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
} 函數返回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
參考：<IEEE Standard 754 for Binary Floating-Point Arithmetic><FAST INVERSE SQUARE ROOT>
最後，給出最精簡的1/sqrt()函數：
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating VALUE
i = 0x5f375a86- (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits BACK to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
} 大家可以嘗試在PC機、51、AVR、430、ARM、上面編譯並實驗，驚訝一下它的工作效率。
前兩天有一則新聞，大意是說 Ryszard Sommefeldt 很久以前看到這麼樣的一段 code (可能出自 Quake III 的 source 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;
}
