平方根倒數速演算法

平方根倒數速演算法的一種演算法。這一演算法的優勢在於減少了求平方根倒數時浮點運算操作帶來的巨大的運算耗費，而在計算機圖形學領域，若要求取照明和投影的波動角度與反射效果，就常需計算平方根倒數。此演算法首先接收一個32位帶符浮點數，然後將之作為一個32位整數看待，以將其向右進行一次邏輯移位的方式將之取半，並用在浮點數規格代表formula_4近似值的十六進位制“-{zh-hans:魔術數字;zh-hant:魔術數字}-”0x5f3759df減之，如此即可得對輸入的浮點數的平方根倒數的首次近似值；而後重新將其作為浮點數，以牛頓法反覆迭代，以求出更精確的近似值，直至求出符合精確度要求的近似值。在計算浮點數的平方根倒數的同一精度的近似值時，此演算法比直接使用浮點數除法要快四倍。此演算法最早被認為是由約翰·卡馬克於90年代前期在的開發中使用，後來則於1999年在《雷神之錘III競技場》的原始碼中應用，但直到2002－2003年間才在Usenet一類的公共論壇上出現。後來的調查顯示，該演算法在這之前就於計算機圖形學的硬體與軟體領域有所應用，如SGI和3dfx就曾在產品中應用此演算法，但至今為止仍未能確切知曉演算法中所使用的特殊常數的起源。演算法的切入點.浮點數的平方根倒數常用於計算正規化向量。3D圖形程式需要使用正規化向量來實現光照和投影效果，因此每秒都需做上百萬次平方根倒數運算，而在處理座標轉換與光源的專用硬體裝置出現前，這些計算都由軟體完成，計算速度亦相當之慢；在1990年代這段程式碼開發出來之時，多數浮點數操作的速度更是遠遠滯後於整數操作，因而針對正規化向量演算法的最佳化就顯得尤為重要。下面陳述計算正規化向量的原理：要將一個向量標準化，就必須計算其歐幾里得範數，以求得向量長度，為此便需對向量的各分量的平方和求平方根；而當求取到其長度，並以之除該向量的每個分量後，所得的新向量就是與原向量同向的單位向量，若以公式表示：formula_5可求得向量v的歐幾里得範數，此演算法正類如對歐幾里得空間的兩點求取其歐幾里得距離，而formula_6求得的就是標準化的向量，若以formula_7代表formula_8，則有formula_9，可見標準化向量時，對向量分量計算平方根倒數實為必需，所以，對平方根倒數計算演算法的最佳化對計算正規化向量也大有裨益。為了加速影象處理單元計算，《雷神之錘III競技場》使用了平方根倒數速演算法，而後來採用現場可程式設計邏輯閘陣列的頂點著色器也應用了此演算法。程式碼概覽.下列程式碼是《雷神之錘III競技場》原始碼中平方根倒數速演算法之例項。示例去除了C預處理器的指令，但附上了原有的註釋：float Q_rsqrt( float number )long i;float x2, y;const float threehalfs = 1.5F;x2 = number * 0.5F;y = number;i = * ( long * ) &amp;y; // evil floating point bit level hackingi = 0x5f3759df - ( i » 1 ); // what the fuck?y = * ( float * ) &amp;i;y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration // y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removedreturn y;為計算平方根倒數的值，軟體首先要先確定一個近似值，而後則使用某些數值方法不斷計算修改近似值，直至達到可接受的精度。在1990年代初，軟體開發時通用的多是從查詢表中取得近似值，而這段程式碼取近似值耗時比之更短，達到精確度要求的速度也比通常使用的浮點除法計算法快四倍，雖然此演算法會損失一些精度，但效能上的巨大優勢已足以補償損失的精度。由程式碼中對原資料的變數型別宣告為float可看出，這一演算法是針對IEEE 754標準格式的32位浮點數設計的，不過據Chris Lomont和後來的Charles McEniry的研究看，這一演算法亦可套用於其他型別的浮點數上。平方根倒數速演算法在速度上的優勢源自將浮點數轉化為長整型以作整數看待，並用特定常數0x5f3759df與之相減。然而對於程式碼閱讀者來說，他們卻難以立即領悟出使用這一常數的目的，因此和其它在程式碼中出現的難以理解的常數一樣，這一常數亦被稱為“-{zh-hans:魔術數字;zh-hant:魔術數字}-”。如此將浮點數當作整數先位移後減法，所得的浮點數結果即是對輸入數字的平方根倒數的粗略估計值，而後再進行一次牛頓迭代法，以使之更精確後，程式碼即執行完畢。由於演算法所生成的用於輸入牛頓法的首次近似值已經相當精確，此演算法所得近似值的精度已可接受，而若使用與《雷神之錘III競技場》同為1999年釋出的Pentium III中的SSE指令rsqrtss計算，則計算平方根倒數的收斂速度更慢，精度也更低。將浮點數轉化為整數.要理解這段程式碼，首先需瞭解浮點數的儲存格式。一個浮點數以32個二進位制位表示一個有理數，而這32位由其意義分為三段：首先首位為符號位，如若是0則為正數，反之為負數；接下來的8位表示經過偏移處理，且formula_16，則可得其表示的浮點數為formula_17。8位二進位制補碼示例如上所述，一個有符號正整數在二進位制補碼系統中的表示中首位為0，而後面的各位則用於表示其數值。將浮點數取儲存為整數時，該整數的數值即為formula_18，其中E表示指數，M表示有效數字；若以上圖為例，圖中樣例若作為浮點數看待有formula_19，formula_20，則易知其轉化而得的整數型號數值為formula_21。由於平方根倒數函式僅能處理正數，因此浮點數的符號位即是使該數的長整形式被2所除。“-{zh-hans:魔術數字;zh-hant:魔術數字}-”.對猜測平方根倒數速演算法的最初構想來說，計算首次近似值所使用的常數0x5f3759df也是重要的線索。為確定程式設計師最初選此常數以近似求取平方根倒數的方法，Charles McEniry首先檢驗了在程式碼中選擇任意常數R所求取出的首次近似值的精度。回想上一節關於整數和浮點數表示的比較：對於同樣的32位二進位制數碼，若為浮點數表示時實際數值為formula_22，而若為整數表示時實際數值則為formula_23，其中formula_24。以下等式引入了一些由指數和有效數字匯出的新元素：formula_25formula_26，其中formula_27再繼續看McEniry 2007裡的進一步說明：formula_28對等式的兩邊取二進位制對數，有formula_31formula_32以如上方法，就能將浮點數x和y的相關指數消去，從而將乘方運算化為加法運算。而由於formula_33與formula_34線性相關，因此在formula_2與formula_36。由此將formula_45代入上式，有：formula_46參照首段等式代入formula_47，formula_48，formula_49與formula_50，有：formula_51移項整理得：formula_52如上所述，對於以浮點規格儲存的正浮點數x，若將其作為長整型表示則示值為formula_23，由此即可根據x的整數表示匯出y的整數表示值，也即：formula_55，其中formula_56.最後匯出的等式formula_57即與上節程式碼中#重定向 一行相契合，由此可見，在平方根倒數速演算法中，對浮點數進行一次移位操作與整數減法，就可以可靠地輸出一個浮點數的對應近似值。到此為止，McEniry只證明了，在常數R的輔助下，可近似求取浮點數的平方根倒數，但仍未能確定程式碼中的R值的選取方法。關於作一次移位與減法操作以使浮點數的指數被-2除的原理，Chris Lomont的論文中亦有個相對簡單的解釋：以formula_58為例，將其指數除-2可得formula_59；而由於浮點表示的指數有進行過偏移處理，所以指數的真實值e應為formula_60，因此可知除法操作的實際結果為formula_61，這時用R減之即可使指數的最低有效數位轉入有效數字域，之後重新轉換為浮點數時，就能得到相當精確的平方根倒數近似值。在這裡對常數R的選取亦有所講究，若能選取一個好的R值，便可減少對指數進行除法與對有效數字域進行移位時可能產生的錯誤。基於這一標準，0xbe即是最合適的R值，而0xbe右移一位即可得到0x5f，這恰是-{zh-hans:魔術數字;zh-hant:魔術數字}-R的第一個位元組。精確度.如上所述，平方根倒數速演算法所得的近似值驚人的精確，右圖亦展示了以上述程式碼計算所得近似值的誤差：當輸入0.01時，以C語言標準庫函式計算可得10.0，而InvSqrt()得值為9.9825822，其間誤差為0.017479，相對誤差則為0.175%，且當輸入更大的數值時，絕對誤差不斷下降，相對誤差也一直控制在一定的範圍之內。牛頓法提高精度.在進行了如上的整數操作之後，示例程式再度將被轉為長整型的浮點數迴轉為浮點數，這裡的浮點運算操作就是對其進行一次牛頓法迭代，若以此例說明：formula_28所求的是y的平方根倒數，以之構造以y為自變數的函式，有formula_63，將其代入牛頓法的通用公式formula_64，有formula_66，其中formula_67，formula_68。整理有formula_69，對應的程式碼即為#重定向 。在以上一節的整數操作產生首次近似值後，程式會將首次近似值作為引數送入函式最後兩句進行精化處理，程式碼中的兩次迭代正是為了進一步提高結果的精度，但由於雷神之錘III引擎的圖形計算中並不需要太高的精度，所以程式碼中只進行了一次迭代，二次迭代的程式碼則被註釋。歷史與考究.《雷神之錘III》的程式碼直到QuakeCon 2005才正式放出，但早在2002年時，平方根倒數速演算法的程式碼就已經出現在Usenet與其他論壇上了。最初人們猜測是卡馬克寫下了這段程式碼，但他在回覆詢問他的郵件時否定了這個觀點，並猜測可能是先前曾幫id Software最佳化雷神之錘的資深彙編程式設計師Terje Mathisen寫下了這段程式碼；而在Mathisen的郵件裡，他表示，在1990年代初，他只曾作過類似的實作，確切來說這段程式碼亦非他所作。現在所知的最早實作是由Gary Tarolli在SGI Indigo中實作的，但他亦坦承他僅對常數R的取值做了一定的改進，實際上他也不是作者。在向以發明MATLAB而聞名的Cleve Moler查證後，Rys Sommefeldt則認為原始的演算法是公司的Greg Walsh所發明，但他也沒有任何決定性的證據能證明這一點。不僅該演算法的原作者不明，人們也仍無法確定當初選擇這個“-{zh-hans:魔術數字;zh-hant:魔術數字}-”的方法。Chris Lomont曾做了個研究：他推算出了一個函式以討論此速演算法的誤差，並找出了使誤差最小的最佳R值0x5f37642f。在Charles McEniry的論文中，他使用了一種類似Lomont但更復雜的方法來最佳化R值：他最開始使用窮舉搜尋，所得結果與Lomont相同；而後他嘗試用帶權二分法尋找最優值，所得結果恰是程式碼中所使用的-{zh-hans:魔術數字;zh-hant:魔術數字}-0x5f3759df，因此，McEniry認為，這一常數最初或許便是以“在可容忍誤差範圍內使用二分法”的方式求得。參考.