我們平時經常會有一些資料運算的操作，需要呼叫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

我們知道，越底層的函式，呼叫越頻繁。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。

參考：<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;
}

他一看之下驚為天人，想要拜見這位前輩高人，但是一路追尋下去卻一直找不到人；同時間也有其他人在找，雖然也沒找到出處，但是 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 (有人整理了一份列表)

好了，故事就到這裡結束了，希望大家能有有收穫:)，我把原始碼也提供下載了，有興趣的朋友們可以自己執行下試試看。

求平方根倒數的演算法

下面這個求  的函式號稱比直接呼叫sqrt庫函式快4倍，來自遊戲Quake III的原始碼。

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

我們這裡分析一下它的原理（指程式的正確性，而不是解釋為何快）。

分析程式之前，我們必須解釋一下float資料在計算機裡的表示方式。一般而言，一個float資料  共32個bit，和int資料一樣。其中前23位為有效數字  ，後面接著一個8位資料  表示指數，最後一位表示符號，由於這裡被開方的數總是大於0，所以我們暫不考慮最後一個符號位。此時

如果我們把計算機內的浮點數  看做一個整數  ，那麼

現在開始逐步分析函式。這個函式的主體有四個語句，分別的功能是：

int i = *(int*)&x; 這條語句把  轉成  。

i = 0x5f3759df - (i>>1); 這條語句從  計算  。

y = *(float*)&i; 這條語句將  轉換為  。

y = y*(1.5f - xhalf*y*y); 這時候的y是近似解；此步就是經典的牛頓迭代法。迭代次數越多越準確。

關鍵是第二步 i = 0x5f3759df - (i>>1); 這條語句從  計算  ，原理:

令  ，用  和  帶入之後兩邊取對數，再利用近似表示  ，算一算就得到

若取  ，  就是程式裡所用的常量0x5f3759df。至於為何選擇這個 ，則應該是曲線擬合實驗的結果。