高效率嵌入式系統(tǒng)開平方根.doc

開平方根開平方根目錄1. 開平方根22. 開平方根說明8181. 開平方根我們平時經(jīng)常會有一些數(shù)據(jù)運(yùn)算的操作，需要調(diào)用sqrt，exp，abs等函數(shù)，那么時候你有沒有想過：這個些函數(shù)系統(tǒng)是如何實(shí)現(xiàn)的？就拿最常用的sqrt函數(shù)來說吧，系統(tǒng)怎么來實(shí)現(xiàn)這個經(jīng)常調(diào)用的函數(shù)呢？雖然有可能你平時沒有想過這個問題，不過正所謂是“臨陣磨槍，不快也光”，你“眉頭一皺，計上心來”，這個不是太簡單了嘛，用二分的方法，在一個區(qū)間中，每次拿中間數(shù)的平方來試驗(yàn)，如果大了，就再試左區(qū)間的中間數(shù)；如果小了，就再拿右區(qū)間的中間數(shù)來試。比如求sqrt(16)的結(jié)果，你先試（0+16）/2=8，8*8=64，64比16大，然后就向左移，試（0+8）/2=4，4*4=16剛好，你得到了正確的結(jié)果sqrt(16)=4。然后你三下五除二就把程序?qū)懗鰜砹耍?用二分法 float SqrtByBisection(float n) /小于0的按照你需要的處理 if(n n)up=mid; else low=mid;last=mid;mid=(up+low)/2; /精度控制 while(abs(mid-last) eps);return mid; 然后看看和系統(tǒng)函數(shù)性能和精度的差別（其中時間單位不是秒也不是毫秒，而是CPU Tick，不管單位是什么，統(tǒng)一了就有可比性）。二分法和系統(tǒng)的方法結(jié)果上完全相同，但是性能上整整差了幾百倍。為什么會有這么大的區(qū)別呢？難道系統(tǒng)有什么更好的辦法？難道。哦，對了，回憶下我們曾經(jīng)的高數(shù)課，曾經(jīng)老師教過我們“牛頓迭代法快速尋找平方根”，或者這種方法可以幫助我們，具體步驟如下。求出根號a的近似值：首先隨便猜一個近似值x，然后不斷令x等于x和a/x的平均數(shù)，迭代個六七次后x的值就已經(jīng)相當(dāng)精確了。例如，我想求根號2等于多少。假如我猜測的結(jié)果為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)的切線來逼近方程x2-a=0的根。根號a實(shí)際上就是x2-a=0的一個正實(shí)根，這個函數(shù)的導(dǎo)數(shù)是2x。也就是說，函數(shù)上任一點(diǎn)(x,f(x)處的切線斜率是2x。那么，x-f(x)/(2x)就是一個比x更接近的近似值。代入 f(x)=x2-a得到x-(x2-a)/(2x)，也就是(x+a/x)/2。相關(guān)的代碼如下：float SqrtByNewton(float x)/ 最終float val = x; / 保存上一個計算的值float last;dolast = val;val =(val + x/val) / 2; while(abs(val-last) eps);return val;牛頓迭代法性能提高了很多，可是和系統(tǒng)函數(shù)相比，還是有這么大差距，這是為什么呀？想啊想啊，想了很久仍然百思不得其解。突然有一天，我在網(wǎng)上看到一個神奇的方法，于是就有了今天的這篇文章，廢話不多說，看代碼先：float InvSqrt(float x)float xhalf = 0.5f*x;int i = *(int*)&x; / get bits for floating VALUE i = 0x5f375a86- (i1); / gives initial guess y0x = *(float*)&i; / convert bits BACK to floatx = x*(1.5f-xhalf*x*x); / Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); / Newton step, repeating increases accuracyx = x*(1.5f-xhalf*x*x); / Newton step, repeating increases accuracyreturn 1/x;這次真的是質(zhì)變了，結(jié)果竟然比系統(tǒng)的還要好。到現(xiàn)在你是不是還不明白那個“鬼函數(shù)”，到底為什么速度那么快嗎？不急，先看看下面的故事吧：Quake-III Arena (雷神之錘3)是90年代的經(jīng)典游戲之一。該系列的游戲不但畫面和內(nèi)容不錯，而且即使計算機(jī)配置低，也能極其流暢地運(yùn)行。這要?dú)w功于它3D引擎的開發(fā)者約翰-卡馬克（John Carmack）。事實(shí)上早在90年代初DOS時代，只要能在PC上搞個小動畫都能讓人驚嘆一番的時候，John Carmack就推出了石破天驚的Castle Wolfstein, 然后再接再勵，doom, doomII, Quake.每次都把3-D技術(shù)推到極致。他的3D引擎代碼資極度高效，幾乎是在壓榨PC機(jī)的每條運(yùn)算指令。當(dāng)初MS的Direct3D也得聽取他的意見，修改了不少API。最近，QUAKE的開發(fā)商ID SOFTWARE 遵守GPL協(xié)議，公開了QUAKE-III的原代碼，讓世人有幸目睹Carmack傳奇的3D引擎的原碼。這是QUAKE-III原代碼的下載地址： /file.x?fid=7547。我們知道，越底層的函數(shù)，調(diào)用越頻繁。3D引擎歸根到底還是數(shù)學(xué)運(yùn)算。那么找到最底層的數(shù)學(xué)運(yùn)算函數(shù)（在game/code/q_math.c）， 必然是精心編寫的。里面有很多有趣的函數(shù)，很多都令人驚奇，估計我們幾年時間都學(xué)不完。在game/code/q_math.c里發(fā)現(xiàn)了這樣一段代碼。它的作用是將一個數(shù)開平方并取倒，經(jīng)測試這段代碼比(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 hackingi = 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#endifreturn y; 函數(shù)返回1/sqrt(x)，這個函數(shù)在圖像處理中比sqrt(x)更有用。注意到這個函數(shù)只用了一次疊代！（其實(shí)就是根本沒用疊代，直接運(yùn)算）。編譯，實(shí)驗(yàn)，這個函數(shù)不僅工作的很好，而且比標(biāo)準(zhǔn)的sqrt()函數(shù)快4倍！要知道，編譯器自帶的函數(shù)，可是經(jīng)過嚴(yán)格仔細(xì)的匯編優(yōu)化的??！ 這個簡潔的函數(shù)，最核心，也是最讓人費(fèi)解的，就是標(biāo)注了“what the fuck?”的一句：i = 0x5f3759df - ( i 1 );再加上y = y * ( threehalfs - ( x2 * y * y ) ); 兩句話就完成了開方運(yùn)算！而且注意到，核心那句是定點(diǎn)移位運(yùn)算，速度極快！特別在很多沒有乘法指令的RISC結(jié)構(gòu)CPU上，這樣做是極其高效的。算法的原理其實(shí)不復(fù)雜,就是牛頓迭代法,用x-f(x)/f(x)來不斷的逼近f(x)=a的根。沒錯，一般的求平方根都是這么循環(huán)迭代算的但是卡馬克(quake3作者)真正牛B的地方是他選擇了一個神秘的常數(shù)0x5f3759df 來計算那個猜測值，就是我們加注釋的那一行，那一行算出的值非常接近1/sqrt(n)，這樣我們只需要2次牛頓迭代就可以達(dá)到我們所需要的精度。好吧如果這個還不算NB，接著看：普渡大學(xué)的數(shù)學(xué)家Chris Lomont看了以后覺得有趣，決定要研究一下卡馬克弄出來的這個猜測值有什么奧秘。Lomont也是個牛人，在精心研究之后從理論上也推導(dǎo)出一個最佳猜測值，和卡馬克的數(shù)字非常接近, 0x5f37642f?？R克真牛，他是外星人嗎？傳奇并沒有在這里結(jié)束。Lomont計算出結(jié)果以后非常滿意，于是拿自己計算出的起始值和卡馬克的神秘數(shù)字做比賽，看看誰的數(shù)字能夠更快更精確的求得平方根。結(jié)果是卡馬克贏了. 誰也不知道卡馬克是怎么找到這個數(shù)字的。最后Lomont怒了，采用暴力方法一個數(shù)字一個數(shù)字試過來，終于找到一個比卡馬克數(shù)字要好上那么一丁點(diǎn)的數(shù)字，雖然實(shí)際上這兩個數(shù)字所產(chǎn)生的結(jié)果非常近似，這個暴力得出的數(shù)字是0x5f375a86。Lomont為此寫下一篇論文，F(xiàn)ast Inverse Square Root。 論文下載地址：/clomont/Math/Papers/2003/InvSqrt.pdf ，/data/InvSqrt.pdf。最后，給出最精簡的1/sqrt()函數(shù)：float InvSqrt(float x)float xhalf = 0.5f*x;int i = *(int*)&x; / get bits for floating VALUE i = 0x5f375a86- (i1); / gives initial guess y0x = *(float*)&i; / convert bits BACK to floatx = x*(1.5f-xhalf*x*x); / Newton step, repeating increases accuracyreturn x; 大家可以嘗試在PC機(jī)、51、AVR、430、ARM、上面編譯并實(shí)驗(yàn)，驚訝一下它的工作效率。前兩天有一則新聞，大意是說 Ryszard Sommefeldt 很久以前看到這么樣的一段 code (可能出自 Quake III 的 source code)：float InvSqrt (float x) float xhalf = 0.5f*x;int i = *(int*)&x;i = 0x5f3759df - (i1);x = *(float*)&i;x = x*(1.5f - xhalf*x*x);return x;他一看之下驚為天人，想要拜見這位前輩高人，但是一路追尋下去卻一直找不到人；同時間也有其他人在找，雖然也沒找到出處，但是 Chris Lomont 寫了一篇論文 (in PDF) 解析這段 code 的算法 (用的是 Newtons Method，牛頓法；比較重要的是后半段講到怎么找出神奇的 0x5f3759df 的)。 PS. 這個 function 之所以重要，是因?yàn)榍?開根號倒數(shù) 這個動作在 3D 運(yùn)算 (向量運(yùn)算的部份) 里面常常會用到，如果你用最原始的 sqrt() 然后再倒數(shù)的話，速度比上面的這個版本大概慢了四倍吧 XD PS2. 在他們追尋的過程中，有人提到一份叫做 MIT HACKMEM 的文件，這是 1970 年代的 MIT 強(qiáng)者們做的一些筆記 (hack memo)，大部份是 algorithm，有些 code 是 PDP-10 asm 寫的，另外有少數(shù)是 C code (有人整理了一份列表)。2. 開平方根說明人們很早就在Quake3源代碼中發(fā)現(xiàn)了類似如下的C代碼，它可以快速的求1/sqrt(x)，在3D圖形向量計算方面應(yīng)用很廣float invSqrt(float x)float xhalf = 0.5 * x;int i = *(int*)&x; / get bits for floating valuei = 0x5f3759df - (i 1); / gives initial guessx = *(float*)&i; / convert bits back to floatx = x * (1.5 - xhalf * x * x); / Newton stepreturn x;在分析這段代碼之前，先看看傳統(tǒng)方法是怎么求一個數(shù)的平方根的倒數(shù)的，一般采用牛頓迭代法，為描述方便，假設(shè)輸入數(shù)為a，顯然需要滿足a0，sqrt是C語言的求平方根函數(shù)，為方便起見，下面用sqrt(x)的形式代替x(1/2)求1/sqrt(a)，用迭代法，即求方程f(x)=x(-2)-a在f(x)=0時的解，選擇適當(dāng)?shù)某跏贾祒0，代入迭代式：x=x-f(x)/ f(x)化簡此式得：x=3x/2-ax3/2這實(shí)際上就是上面函數(shù)倒數(shù)第二行，從函數(shù)注釋也可以直接看出，這一步就是牛頓迭代，一般選擇一個合適的初始值開始迭代后，迭代次數(shù)越多越接近解，換句話說就是精度越高，誤差越小，當(dāng)誤差小于可接受值，即可獲得近似結(jié)果了，就這個問題而言，初始值的選擇一般要在區(qū)間(0, sqrt(3/a)，證明從略而invSqrt這個函數(shù)厲害的地方就在于，在正式迭代開始前的三行計算已經(jīng)得到了一個非常接近于解的數(shù)，因此只需一次迭代，即可得到近似值，經(jīng)測試，對于常用的浮點(diǎn)數(shù)范圍，invSqrt(x)與標(biāo)準(zhǔn)解1/sqrt(x)的最大相對誤差為1.75，平均相對誤差為0.95，這個精度在很多時候已經(jīng)滿足QuakeIII的基本要求1了，而invSqrt(x)的速度則比直接計算1/sqrt(x)快4倍，這對于Quake和CS這類游戲的性能是非常重要的，而且如果需要更高的精度，則將迭代那一行再重復(fù)一次就可以了，相對誤差會降到百萬分之一的級別，只不過速度會慢一些，接下來我們來分析下這三行代碼的原理最令人迷惑的是gives initial guess這行，即對i做了移位和減法運(yùn)算，不過熟悉C語言的人應(yīng)該能看出來，這個算法和float浮點(diǎn)數(shù)的內(nèi)部表示有關(guān)，分析應(yīng)該從這里入手正式開始前先輕松下，講些歷史故事，人們在QuakeIII源碼發(fā)現(xiàn)了這個函數(shù)，于是很自然的認(rèn)為這是卡馬克（John Carmack）的杰作，其中0x5f3759df這個數(shù)被稱為卡馬克密碼，我們在下面稱這個數(shù)為magic，Beyond3D.com的Ryszard Sommefeldt一直在想到底是哪個家伙寫了這些神奇的代碼，于是就開始找作者，John Carmack在郵件回復(fù)中明確表示不是他，也不是Michael。Terje Mathisen說他寫過類似的高效代碼，但上面的不是。后來猜測這個來自于一些早期黑客的算法筆記，作者究竟是誰自然也難以追查了，可以肯定的是這個家伙對計算機(jī)和高數(shù)知識都有較好理解，很聰明2003年普渡大學(xué)的數(shù)學(xué)家Chris Lomont寫了一篇文章對這段代碼進(jìn)行了分析。論文是英文的，地址在：/view/80b84d1fb7360b4c2e3f644b.html在這篇12頁的論文中，Lomont對這個算法做了分析，并從推導(dǎo)出了一個理論上最優(yōu)的magic 0x5f37642f，有意思的是，這個數(shù)居然沒有invSqrt里的0x5f3759df效果好，最大相對誤差達(dá)到1.78，Lomont一怒之下，用暴力搜索枚舉了所有可能的magic，終于找到一個最優(yōu)的magic 0x5f375a86，只比0x5f3759df效果好一點(diǎn)點(diǎn)，至于invSqrt的作者究竟如何找到0x5f3759df的，也就是個迷了開始正式分析，這三行代碼是把float在內(nèi)存中的表示作為一個整數(shù)i看待，然后對i進(jìn)行一次移位和減法，然后再將i的值作為一個float看待，所以我們先看看float在內(nèi)存中的表示，一般計算機(jī)的浮點(diǎn)數(shù)遵循IEEE754標(biāo)準(zhǔn)，采用以2為底數(shù)的科學(xué)計數(shù)法，例如二進(jìn)制的11010.11001記為1.101011001*10100，float和int都是32位，占4個字節(jié)（注意，這個函數(shù)早期的代碼中整數(shù)類型應(yīng)該是long，因?yàn)槟菚r候在dos下，int是16位的）：最高位d31：符號位，0表示非負(fù)，1表示非正，為什么不直接說正負(fù)呢，因?yàn)橛袛?shù)值0的存在，浮點(diǎn)數(shù)有+0和-0的區(qū)別，這個位用S表示d30d23：指數(shù)域，存放一個整數(shù)，表示127+E，E為指數(shù)，由于指數(shù)域的范圍是0255，因此理論上可以表示的指數(shù)范圍是-127128，不過0和255有特殊含義，所以范圍實(shí)際要稍微小一點(diǎn)，這個先按下不表，我們認(rèn)為常用浮點(diǎn)數(shù)不包括這兩種極端情況d22d0：有效數(shù)字域，只是小數(shù)部分，由于科學(xué)計數(shù)法的規(guī)定，整數(shù)部分肯定是1，就省略了，這樣可以避免不必要的精度浪費(fèi)，為描述方便，這個域所表示的小數(shù)設(shè)為F。當(dāng)然有人會問，那0怎么辦，+0和-0有自身的特殊表示法，S位表示符號，其他位都為0的時候是+0或-0于是除去0和IEEE754規(guī)定的特殊值，一個常用浮點(diǎn)數(shù)的表示可以看做：(-1)S*(1+F)*2E，具體到我們需要分析的問題，由于輸入是正數(shù)，S位肯定是0，就不做考慮了，簡化為：(1+F)*2E好，現(xiàn)在我們需要求(1+F)*2E的平方根的倒數(shù)，即求1/sqrt(1+F)*2E)，求得的結(jié)果當(dāng)然也要用這個浮點(diǎn)數(shù)表示法，有效數(shù)字必須在1, 2)，指數(shù)域?yàn)檎麛?shù)，則結(jié)果分兩種情況推導(dǎo)出結(jié)果：E為奇數(shù)：sqrt(2/(1+F)*2(-(E+1)/2)E為偶數(shù)：2/sqrt(1+F)*2(-E/2-1)先看指數(shù)，如果我們需要通過計算機(jī)的整數(shù)運(yùn)算（移位、位運(yùn)算和加減法等）來逼近解，首先要在數(shù)量級上盡量靠近，因?yàn)橹灰獢?shù)量級一樣，兩個數(shù)的誤差范圍是最小的，也就是說，我們需要：E為奇數(shù)時，將127+E變成127-(E+1)/2E為偶數(shù)時，將127+E變成127-E/2-1于是，通過右移一位來實(shí)現(xiàn)除以二，通過用一個數(shù)減去指數(shù)域來將E變成負(fù)的，這樣invSqrt中的那句就很容易理解了：E為奇數(shù)時，127+E為偶數(shù)，右移等于除以二，190-(127+E)/2 = 190-63-(E+1)/2 = 127-(E+1)/2E為偶數(shù)時，127+E為奇數(shù)，右移等于先減一再除以二，189-(127+E-1)/2 = 189-63-E/2 = 127-E/2-1然后我們把invSqrt中的magic 0x5f3759df用float的形式展開，則其指數(shù)域?yàn)?x5f1這個操作將指數(shù)域最后一位的1（即上面說的“先減一”）也向右移了一位，i右移后d22位為1，這樣一來只要被減數(shù)的d22這一位是0，就會因?yàn)椴粔驕p而產(chǎn)生借位，指數(shù)域被借了1，自然就變成189了細(xì)心的童鞋應(yīng)該發(fā)現(xiàn)了，當(dāng)E為偶數(shù)時，減法做完后，指數(shù)域一定是127-E/2-1，但是如果E為奇數(shù)，則不保證指數(shù)域是127-(E+1)/2，因?yàn)檫@時候被減數(shù)和減數(shù)的d22位都是0，但如果減數(shù)的d21到d0這個域的數(shù)字比被減數(shù)的大，就會產(chǎn)生借位而使得指數(shù)域比預(yù)期的127-(E+1)/2要小1，這個問題會導(dǎo)致一定的誤差，但是在最后的迭代中誤差會被縮小，這個誤差具體有多大，這里就不詳細(xì)討論了，有興趣的童鞋可以自己算算看我們還是先證明用這種方法得到的x0是落在上述區(qū)間(0, sqrt(3/a)的，這里a就是(1+F)*2E由于輸入是一個非負(fù)數(shù)，則S為0，指數(shù)域肯定不為0（右移后不可能剛好為190），因此x0肯定大于0，我們將sqrt(3/a)展開成期望的解的指數(shù)的乘法形式：E為奇數(shù)時，展開為sqrt(6/(1+F)*2(-(E+1)/2)，如果按上面說的那種情況產(chǎn)生借位，則展開為sqrt(24/(1+F)*2(-(E+3)/2)E為偶數(shù)時，展開為sqrt(12/(1+F)*2(-E/2-1)可以看到，無論是哪種情況，在指數(shù)相同的情況下，有效數(shù)字都大于2，反過來說，用上述算法得到的x0是小于sqrt(3/a)的，而且還小了很多，非常接近解，這時候只需要一次迭代，就得到了誤差很小的近似結(jié)果到這里，基本原理都清楚了，只要我們保證x0在這個區(qū)間中，再做迭代總是能進(jìn)一步接近解的，現(xiàn)在的問題就在于magic中d21到d0這個域的值應(yīng)該怎么取了，這個取值關(guān)系到每次迭代的誤差，比如說，我們?nèi)?，這樣也避免了E為奇數(shù)時的借位情況，這樣magic就是0x5f3fffff，用這個magic測試，結(jié)果最大誤差超過1%，平均誤差超過6，顯然效果太差了我們假設(shè)將magic作為float看時，小數(shù)部分的值是M，由于d22位已確定是0，則0=M1)的減法運(yùn)算中，d22到d0域的運(yùn)算可看做是定點(diǎn)小數(shù)減法，分三種情況討論：E為奇數(shù)時，減數(shù)的d22位為0，則小數(shù)部分的值為F/2，因?yàn)橛乙茖π?shù)來說也是除以二，d0位如果是1則會舍棄，這個因?yàn)樘《雎?，假設(shè)M=F/2，則不需要向指數(shù)域借位，計算結(jié)果的小數(shù)部分為M-F/2，指數(shù)域符合結(jié)果預(yù)期，此時的相對誤差為：rd1(M,F)=|1-(1+M-F/2)/sqrt(2/(1+F)|，0=F=2M假設(shè)MF/2，則需要向指數(shù)域借位，計算結(jié)果小數(shù)部分為1+M-F/2，指數(shù)位比預(yù)期低1，計算誤差的時候需要將指數(shù)域差值補(bǔ)回去，此時的相對誤差：rd2(M,F)=|1-(2+M-F/2)/2/sqrt(2/(1+F)|，2MF1E為偶數(shù)時，減數(shù)的d22位為1，則小數(shù)部分的值為F/2+1/2，此時必定借位，相對誤差為：rd3(M,F)=|1-(2+M-F/2-1/2)/2/sqrt(1+F)|，0=F 1;x = *(float *)&i;x = (x + a / x) * 0.5;return x;推導(dǎo)過程從略，有興趣的童鞋可以自行研究，注意這里是先加再右移，所以用unsigned int防止負(fù)數(shù)右移，這個mySqrt計算平方根速度比invSqrt(a)*a快，而且最大誤差只有0.6，那么，有沒有更快的算法呢，如果單純用計算，可能很難超越mySqrt了，想要更快，得從另外的方向想辦法比較直接的想法是，造一個巨大的table，儲存a到sqrt(a)的映射關(guān)系，這樣不需要計算，只要查表就行，而且由于