一個Sqrt函式引發的血案
我們平時經常會有一些資料運算的操作，需要呼叫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原始碼的下載地址：
我們知道，越底層的函式，呼叫越頻繁。3D引擎歸根到底還是數學運算。那麼找到最底層的數學運算函式（在game/code/q_math.c），必然是精心編寫的。裡面有很多有趣的函式，很多都令人驚奇，估計我們幾年時間都學不完。在game/code/q_math.c裡發現了這樣一段程式碼。它的作用是將一個數開平方並取倒，經測試這段程式碼比(float)(1.0/sqrt(x))快4倍：
float Q_rsqrt( float number )
{
longi;
floatx2, y;
constfloat threehalfs = 1.5F;
x2= number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what thefuck?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); //1st iteration
//y = y * ( threehalfs - ( x2 * y * y )); // 2nd iteration, this can be removed
returny;
}
函式返回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。
最後，給出最精簡的1/sqrt()函式：
float InvSqrt(float x)
{
floatxhalf = 0.5f*x;
inti = *(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
returnx;
}前兩天有一則新聞，大意是說 Ryszard Sommefeldt 很久以前看到這麼樣的一段code (可能出自 Quake III 的 source code)：
float InvSqrt (float x)
{
floatxhalf = 0.5f*x;
inti = *(int*)&x;
i= 0x5f3759df - (i>>1);
x= *(float*)&i;
x= x*(1.5f - xhalf*x*x);
returnx;
}PS. 這個 function 之所以重要，是因為求開根號倒數這個動作在 3D 運算 (向量運算的部份) 裡面常常會用到，如果你用最原始的 sqrt() 然後再倒數的話，速度比上面的這個版本大概慢了四倍吧… XD原始碼下載地址：
問題出在我簽入的來自卡馬克的求平方根函式程式碼。
double InvSqrt(double number)
{
__int64 i;
double x2, y;
const double threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(__int64 *)&y;
i = 0x5fe6ec85e7de30da - (i >> 1);
y = *( double *)&i;
y = y * (threehalfs - (x2 * y * y)); //1stiteration
y = y * (threehalfs - (x2 * y * y)); //2nditeration, this can be removed
return y;
}紅色部分程式碼在gcc開啟-fstrict-aliasing選項後將得到錯誤的程式碼。由於使用了type-punned pointer將打破strict-aliasing規則。
由於-fstrict-aliasing選項在-O2, -O3, -Os等優化模式下都將開啟（目前dev不帶優化，main帶-O3所以該問題只在main上出現）所以建議對linux編譯中產生
warning:dereferencing type-punned pointer will break strict-aliasing rule
警告的情況作為編譯失敗。以便防止出現類似問題。
上述程式碼應當使用聯合體重寫為：
double InvSqrt(double number)
{
double x2, y;
const double threehalfs = 1.5F;
union
{
double d;
__int64 i;
}d;
x2 = number * 0.5F;
y = number;
d.d = y;
d.i = 0x5fe6ec85e7de30da - (d.i >> 1);
y = d.d;
y = y * (threehalfs - (x2 * y * y)); //1stiteration
y = y * (threehalfs - (x2 * y * y)); //2nditeration, this can be removed
return y;
}這樣就不會打破該規則。
平方根倒數速演算法
平方根倒數速演算法（英語：Fast Inverse Square Root，亦常以“FastInvSqrt()”或其使用的十六進位制常數0x5f3759df代稱）是用於快速計算（即的平方根的倒數，在此需取符合標準格式的32位浮點數）的一種演算法。此演算法最早可能是於90年代前期由SGI所發明，後來則於1999年在《雷神之錘III競技場》的原始碼中應用，但直到2002－2003年間才在一類的公共論壇上出現。這一演算法的優勢在於減少了求平方根倒數時浮點運算操作帶來的巨大的運算耗費，而在領域，若要求取照明和投影的波動角度與反射效果，就常需計算平方根倒數。
此演算法首先接收一個32位帶符浮點數，然後將之作為一個32位整數看待，以將其向右進行一次邏輯移位的方式將之取半，並用十六進位制“魔術數字”0x5f3759df減之，如此即可得對輸入的浮點數的平方根倒數的首次近似值；而後重新將其作為浮點數，以牛頓法反覆迭代，以求出更精確的近似值，直至求出符合精確度要求的近似值。在計算浮點數的平方根倒數的同一精度的近似值時，此演算法比直接使用浮點數除法要快四倍。
1介紹
平方根倒數速演算法最早被認為是由約翰·卡馬克所發明，但後來的調查顯示，該演算法在這之前就於計算機圖形學的硬體與軟體領域有所應用，如SGI和3dfx就曾在產品中應用此演算法。而就現在所知，此演算法最早由Gary Tarolli在SGIIndigo的開發中使用。雖說在隨後的相關研究中也提出了一些可能的來源，但至今為止仍未能確切知曉此常數的起源。
2演算法的切入點
浮點數的平方根倒數常用於計算正規化向量。[1]3D圖形程式需要使用正規化向量來實現光照和投影效果，因此每秒都需做上百萬次平方根倒數運算，而在處理座標轉換與光源的專用硬體裝置出現前，這些計算都由軟體完成，計算速度亦相當之慢；在1990年代這段程式碼開發出來之時，多數浮點數操作的速度更是遠遠滯後於整數操作，因而針對正規化向量演算法的優化就顯得尤為重要。下面陳述計算正規化向量的原理：
要將一個向量標準化，就必須計算其歐幾里得範數以求得向量長度，而這時就需對向量的各分量的平方和求平方根；而當求取到其長度並以之除該向量的每個分量後，所得的新向量就是與原向量同向的單位向量，若以公式表示：
· 可求得向量v的歐幾里得範數，此演算法正類如對歐幾里得空間的兩點求取其歐幾里得距離，
· 而求得的就是標準化的向量，若以代表，則有，
可見標準化向量時需要用到對向量分量的平方根倒數計算，所以對平方根倒數計算演算法的優化對計算正規化向量也大有裨益。
為了加速影象處理單元計算，《雷神之錘III競技場》使用了平方根倒數速演算法，而後來採用現場可程式設計邏輯閘陣列的頂點著色器也應用了此演算法。[2]
3將浮點數轉化為整數
要理解這段程式碼，首先需瞭解浮點數的儲存格式。一個浮點數以32個二進位制位表示一個有理數，而這32位由其意義分為三段：首先首位為符號位，如若是0則為正數，反之為負數；接下來的8位表示經過偏移處理（這是為了使之能表示-127－128）後的指數；最後23位表示的則是有效數字中除最高位以外的其餘數字。將上述結構表示成公式即為，其中表示有效數字的尾數（此處，偏移量，而指數的值決定了有效數字（在Lomont和McEniry的論文中稱為“尾數”（mantissa））代表的是小數還是整數。以上圖為例，將描述帶入有），且，則可得其表示的浮點數為。
|
符號位 |
|||||||||
|
0 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
= |
127 |
|
0 |
0 |
0 |
0 |
0 |
0 |
1 |
0 |
= |
2 |
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
= |
1 |
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
= |
0 |
|
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
= |
−1 |
|
1 |
1 |
1 |
1 |
1 |
1 |
1 |
0 |
= |
−2 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
= |
−127 |
|
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
= |
−128 |
8位二進位制整數補碼示例
如上所述，一個有符號正整數在二進位制補碼系統中的表示中首位為0，而後面的各位則用於表示其數值。將浮點數取別名儲存為整數時，該整數的數值即為，其中E表示指數，M表示有效數字；若以上圖為例，圖中樣例若作為浮點數看待有，，則易知其轉化而得的整數型號數值為。由於平方根倒數函式僅能處理正數，因此浮點數的符號位（即如上的Si）必為0，而這就保證了轉換所得的有符號整數也必為正數。以上轉換就為後面的計算帶來了可行性，之後的第一步操作（邏輯右移一位）即是使該數的長整形式被2所除。[3]
4歷史與考究
id Software創始人約翰·卡馬克
《雷神之錘III》的程式碼直到QuakeCon 2005才正式放出，但早在2002年（或2003年）時平方根倒數速演算法的程式碼就已經出現在與其他論壇上了。最初人們猜測是卡馬克寫下了這段程式碼，但他在詢問郵件的回覆中否定了這個觀點，並猜測可能是先前曾幫id Software優化雷神之錘的資深彙編程式設計師Terje Mathisen寫下了這段程式碼；而在Mathisen的郵件裡他表示在1990年代初他只曾作過類似的實現，確切來說這段程式碼亦非他所作。現在所知的最早實現是由Gary Tarilli在SGIIndigo中實現的，但他亦坦承他僅對常數R的取值做了一定的改進，實際上他也不是作者。RysSommefeldt則在向以發明而聞名的CleveMoler查證後認為原始的演算法是Ardent Computer公司的GregWalsh所發明，但他也沒有任何決定性的證據能證明這一點。
目前不僅該演算法的原作者不明，人們也仍無法明確當初選擇這個“魔術數字”的方法。Chris Lomont在研究中曾做了個試驗：他編寫了一個函式，以在一個範圍內遍歷選取R值的方式將逼近誤差降到最小，以此方法他計算出了線性近似的最優R值0x5f37642f（與程式碼中使用的0x5f3759df相當接近），但以之代入演算法計算並進行一次牛頓迭代後，所得近似值與代入0x5f3759df的結果相比精度卻仍略微更低；而後Lomont將目標改為遍歷選取在進行1－2次牛頓迭代後能得到最大精度的R值，並由此算出最優R值為0x5f375a86，以此值代入演算法並進行牛頓迭代後所得的結果都比代入原始值（0x5f3759df）更精確，於是他的論文最後以“原始常數是以數學推導還是以反覆試錯的方式求得”的問題作結。在論文中Lomont亦指出64位的IEEE754浮點數（即雙精度型別）所對應的魔術數字是0x5fe6ec85e7de30da，但後來的研究表明代入0x5fe6eb50c7aa19f9的結果精確度更高（McEniry得出的結果則是0x5FE6EB50C7B537AA，精度介於兩者之間）。在Charles McEniry的論文中，他使用了一種類似Lomont但更復雜的方法來優化R值：他最開始使用窮舉搜尋法，所得結果與Lomont相同；而後他嘗試用帶權二分法尋找最優值，所得結果恰是程式碼中所使用的魔術數字0x5f3759df，因此McEniry確信這一常數或許最初便是以“在可容忍誤差範圍內使用二分法”的方式求得。[4]
5註釋
1. 由於現代計算機系統對長整型的定義有所差異，使用長整型會降低此段程式碼的可移植性。具體來說，由此段浮點轉換為長整型的定義可知，如若這段程式碼正常執行，所在系統的長整型長度應為4位元組（32位），否則重新轉為浮點數時可能會變成負數；而由於C99標準的廣泛應用，在現今多數64位計算機系統（除使用LLP64資料模型的Windows外）中，長整型的長度都是8位元組。[5]
2. 此處“浮點數”所指為標準化浮點數，也即有效數字部分必須滿足，可參見DavidGoldberg. What Every Computer Scientist Should Know About Floating-PointArithmetic. ACM Computing Surveys. 1991.March, 23 (1): 5–48. doi:10.1145/103162.103163.
3. Lomont 2003確定R的方式則有所不同，他先將R分解為與，其中與分別代表R的有效數字域和指數域。[6]
---平方根倒數演算法 [轉]
--------------------------------------------------------------------------------
快速平方根（平方根倒數）演算法
日前在書上看到一段使用多項式逼近計算平方根的程式碼，至今都沒搞明白作者是怎樣推算出那個公式的。但在嘗試解決問題的過程中，學到了不少東西，於是便有了這篇心得，寫出來和大家共享。其中有錯漏的地方，還請大家多多指教。
的確，正如許多人所說的那樣，現在有有FPU，有3DNow，有SIMD，討論軟體演算法好像不合時宜。關於sqrt的話題其實早在2003年便已在 GameDev.net上得到了廣泛的討論（可見我實在非常火星了，當然不排除還有其他尚在冥王星的人，嘿嘿）。而嘗試探究該話題則完全是出於本人的興趣和好奇心（換句話說就是無知）。
我只是個beginner，所以這種大是大非的問題我也說不清楚（在GameDev.net上也有很多類似的爭論）。但無論如何，Carmack在DOOM3中還是使用了軟體演算法，而多知道一點數學知識對3D程式設計來說也只有好處沒壞處。3D圖形程式設計其實就是數學，數學，還是數學。
=========================================================
在3D圖形程式設計中，經常要求平方根或平方根的倒數，例如：求向量的長度或將向量歸一化。C數學函式庫中的sqrt具有理想的精度，但對於3D遊戲程式來說速度太慢。我們希望能夠在保證足夠的精度的同時，進一步提高速度。
Carmack在QUAKE3中使用了下面的演算法，它第一次在公眾場合出現的時候，幾乎震住了所有的人。據說該演算法其實並不是Carmack發明的，它真正的作者是Nvidia的GaryTarolli（未經證實）。－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
//
// 計算引數x的平方根的倒數
//
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;
}－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－該演算法的本質其實就是牛頓迭代法（Newton-RaphsonMethod，簡稱NR），而NR的基礎則是泰勒級數（Taylor Series）。NR是一種求方程的近似根的方法。首先要估計一個與方程的根比較靠近的數值，然後根據公式推算下一個更加近似的數值，不斷重複直到可以獲得滿意的精度。其公式如下：－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－函式：y=f(x)
其一階導數為：y'=f'(x)
則方程：f(x)=0 的第n+1個近似根為
x[n+1] = x[n] - f(x[n]) / f'(x[n])
－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
NR最關鍵的地方在於估計第一個近似根。如果該近似根與真根足夠靠近的話，那麼只需要少數幾次迭代，就可以得到滿意的解。
現在回過頭來看看如何利用牛頓法來解決我們的問題。求平方根的倒數，實際就是求方程1/(x^2)-a=0的解。將該方程按牛頓迭代法的公式展開為：
x[n+1]=1/2*x[n]*(3-a*x[n]*x[n])
將1/2放到括號裡面，就得到了上面那個函式的倒數第二行。
接著，我們要設法估計第一個近似根。這也是上面的函式最神奇的地方。它通過某種方法算出了一個與真根非常接近的近似根，因此它只需要使用一次迭代過程就獲得了較滿意的解。它是怎樣做到的呢？所有的奧妙就在於這一行：
i = 0x5f3759df - (i >> 1); // 計算第一個近似根
超級莫名其妙的語句，不是嗎？但仔細想一下的話，還是可以理解的。我們知道，IEEE標準下，float型別的資料在32位系統上是這樣表示的（大體來說就是這樣，但省略了很多細節，有興趣可以GOOGLE）：－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
bits：31 30... 0
31：符號位
30-23：共8位，儲存指數（E）
22-0：共23位，儲存尾數（M）－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－所以，32位的浮點數用十進位制實數表示就是：M*2^E。開根然後倒數就是：M^(-1/2)*2^(-E/2)。現在就十分清晰了。語句i>
>1其工作就是將指數除以2，實現2^(E/2)的部分。而前面用一個常數減去它，目的就是得到M^(1/2)同時反轉所有指數的符號。
至於那個0x5f3759df，呃，我只能說，的確是一個超級的MagicNumber。
那個Magic Number是可以推匯出來的，但我並不打算在這裡討論，因為實在太繁瑣了。簡單來說，其原理如下：因為IEEE的浮點數中，尾數