檢視反平方根快速演算法的原始碼

'''反平方根快速演算法'''（英語：Fast Inverse Square Root，亦常以「Fast InvSqrt ( )」抑是其使用的十六進位常數'''零 x 五 f 三千七百五十九 df'''代稱）是用佇咧快速計算 $ \ textstyle x ^ {-二分之一 } $（即 $ \ textstyle x $ 的平方根的倒算，在此 $ \ textstyle x $ 需要符合 IEEE 七百五十四標準格式的三十二位元浮點數）的一種演算法。這演算法的優勢佇咧減少求平方根倒算時浮點運算操作帶來的真大的運算了，啊若佇電腦圖學領域，若要求取照明佮投影的波動角度佮反射效果，就常在需要算平方根倒算。

此演算法首先接收一个三十二位元帶符浮點數，然後共之作為一个三十二位元整數看待，以共正面進行一擺邏輯移位的方式共取半，並且用佇咧浮點數規格代表 √ 二千一百二十七近似值的十六進位「魔術數字」零 x 五 f 三千七百五十九 df 減之，按呢就會當對輸入的浮點數的平方根上尾仔幾若值；後來重新共其做浮點數，用牛做法迵天，以求出閣較精確的近若值，一直到求出符合精確度要求的近似值。佇計算浮點數的平方根倒算的仝一精度的近來親像值時，若這个演算法比直接使用浮點數除法著愛緊四倍。

此演算法較早予人認為是由約翰 ・ 卡馬克於九零年進前 SGI Indigo 的開發中使用，後來就是一九九九年佇咧《雷神之槌 III 競技場》的原始碼中應用，毋過一直到二千空二－二空空三年才佇 Usenet 一類的公共論壇頂面出現。後來的調查顯示，該演算法佇這進前就佇電腦圖學的硬體佮軟體領域有所應用，如 SGI 佮三 dfx 就捌佇產品應用這馬算法，毋過到今猶是袂當確切知影演算法內底所使用的特殊常數的起源。

==演算法的切入點==

浮點數的平方根倒算捷用佇計算正規化向量。三 D 圖形程式需要使用正規化向量來實現光照佮投影效果，所以每一秒攏需要百萬擺平方根倒算來做，咧處理坐標轉換佮光源的專用硬體裝置出現進前，遮的計算攏是軟體完成，計算速度亦相當之慢；佇一九九空年代這段代碼開發出來之時，多數浮點數操作的速度閣較是遠遠無後就開始操作，因為重要針對正規化向量演算法的最佳化就顯得尤為重要。下跤陳述計算正規化向量的原理：

欲將一个向量標準化，就必須愛計算其歐幾里著愛算，以求得向量長度，為此便需要對向量的各分量的平方佮求平方根；所以求取著其長度，除了這向量的每一个分量了後，所得的新向量就是佮原向量仝向的單位向量，若用公式表示：


: $ \ | { \ boldsymbol { v } } \ |={ \ sqrt { v _ { 一 } ^ { 二 } + v _ { 二 } ^ { 二 } + v _ { 三 } ^ { 二 } } } $ 會當求甲向量 v 的歐幾里著範數，此演算法正類如對歐幾里得空間的兩點求取其歐幾里得距離，


: 而且 $ { \ boldsymbol { \ hat { v } } }={ \ boldsymbol { v } } / \ | { \ boldsymbol { v } } \ | $ 求得的就是標準化的向量，若以 $ { \ boldsymbol { x } } $ 代表 $ v _ { 一 } ^ { 二 } + v _ { 二 } ^ { 二 } + v _ { 三 } ^ { 二 } $，則有 $ { \ boldsymbol { \ hat { v } } }={ \ boldsymbol { v } } / { \ sqrt { x } } $，

可見標準化向量的時陣，對向量分算平方根倒算實為必需需，所以乎，對平方根倒算演算法的最佳化對計算正規化向量嘛有大有標益。

為著加速圖像處理單元計算，《 雷神之槌 III 競技場》使用矣反平方根快速演算法，尾仔採用現場的可程式化邏輯閘陣列的頂點著色器也應用這樣演算法。

==代碼概覽==

下列代碼是《雷神之槌 III 競技場》原始碼中反平方根快速演算法之實例。範例去除了 C 預處理器的指令，但是附上矣原有的注釋：

為計算平方根倒算的值，軟體首先愛先確定一个近似值，煞尾使用某寡數值的方法一直算修改較近的若值，一直到達到會當接受的精度。佇一九九空年代初（嘛即時演算法發明的大概時間）， 軟體開發的時陣通用的平方根算法濟是對走揣表中取得近似值，這段代碼取近來親像值了時比較短，達到精確度要求的速度嘛比通常使用的浮點除法計演算法咧欲四倍，雖然講這个演算法會損失一寡精度，但效能上誠大優勢已經誠以補償損失的精度。由代碼中對原資料的變數類型聲明為 float 看出講，這演算法是針對 IEEE 七百五十四標準格式的三十二位元浮點數設計的，毋過根據 Chris Lomont 佮後來的 Charles McEniry 的研究看，這演算法亦可套用於其他類型的浮點數上。

反平方根快速演算法佇速度上的優勢源自將浮點數轉化做長款以作整數看待，並用特定定數零 x 五 f 三千七百五十九 df 佮之相減。毋過對代碼閱讀者來講，𪜶煞歹隨領悟出使用這一常數的目的，所以佮其他佇代碼中出現的僫理解的常數仝款，這一常數亦被稱為「魔術數字」。 按呢共浮點數當做是整數先位徙後減法，所得的浮點數結果就是對輸入數字的平方根抑算了的粗略估計值，了後閣進行一擺牛做迵天代法，會使閣較精確了後，代碼即執行了畢。因為演算法所生的用於輸入牛頓法的頭一擺近若值已經誠精確，這个演算法所得近來若值的精度已經會當接受，若使用佮《雷神之槌 III 競技場》仝為一九九九年發布的 Pentium III 中的 SSE 指令 rsqrtss 計算，則計算平方根倒算的收斂速度閣較慢，精度嘛閣較低。

===共浮點數轉化做整數===

愛理解講這段代碼，首先需要了解浮點數的儲存格式。一个浮點數以三十二个位元表示一个有理數，這三十二位元對其意義分做三段：首先首位做符號位，如果若是零則為正數，反之為負數；紲落來的八位元表示經過偏移處理（這是為著使用表示-一百二十七喔－一百二十八）後的指數；上尾仔二十三位表示的就是有效數字內底除上懸位以外的數字其他的數字。欲共結構表示成公式即為 $ \ textstyle x=( 影一 ) ^ { \ mathrm { Si } } \ cdot ( 一 + m ) \ cdot 二 ^ { ( E-B ) } $，其中 $ \ textstyle m $ 表示有效數字的尾數（此處 $ \ textstyle 零 \ leq m < 一 $，偏徙量 $ \ textstyle B=一百二十七喔 $，指數的值 $ \ textstyle E-B $ 決定有效數字（佇咧 Lomont 和 McEniry 的論文中講做「尾數」（_ mantissa _）） 代表的是小數抑是整數。以上圖為例，欲來講帶入去有 $ \ textstyle m=一 \ times 二 ^ { 鋪二 }=空九二五空 $）， 而且 $ \ textstyle E-B=一百二十四孵一百二十七=ma三 $，著愛會當表示的浮點數為 $ \ textstyle x=( 一 + 空九二五空 ) \ cdot 二 ^ { ma三 }=空九一五六二五 $。

如何咧講，一个有符號正整數佇咧二補數系統中的表示中頭一个為零，若後壁的各位用表示其數值。共浮點數號別名儲存做整數時，該整數的數值即為 $ \ textstyle I=E \ times 二 ^ { 二十三 } + M $，其中 E 表示指數，M 表示有效數字；若以上圖為例，圖內底看出來若做浮點數看待有 $ \ textstyle E=百二四 $，$ M=一 \ cdot 二 ^ { 二十一 } $，好知影其轉化會著的整數型號數值為 $ I=百二四 \ times 二 ^ { 二十三 } + 二 ^ { 二十一 } $。因為平方根倒算函式干焦會當處理正數，因此浮點數的符號位（即如上的 Si）必為零，這就保證著轉換所得的有符號整數嘛必為正數。以上轉換就為後壁的計算帶來矣可行性，了後的頭一步操作（邏輯正徙一个）也就是該有幾若个形式被二所除了。

===「魔術數字」===

著猜測反平方根快速演算法的上頭構想來講，算頭擺近若值所使用的常數'''零 x 五 f 三千七百五十九 df'''嘛是重要的線索。為著確定程式設計師上早選此常數以近好親像求取平方根倒數的方法，Charles McEniry 代先檢定矣佇代碼內面選擇任意常數 R 所求取出的頭一擺近來若值的精度。回想頂節關於整數佮浮點數表示較：對仝款的三十二位元二進位數位，若為浮點數表示的時實際數值為 $ \ textstyle x=( 一 + m _ { x } ) 二 ^ { e _ { x } } $，若為著整數表示時實際數值則為 $ \ textstyle I _ { x }=E _ { x } L + M _ { x } $，其中 $ \ textstyle L=二 ^ { n 影一-b } $。以下等式引入一寡由指數和有效數字匯出的新元素：


: $ m _ { x }={ \ frac { M _ { x } } { L } } $


: $ e _ { x }=E _ { x }-B $，其中 $ B=二 ^ { b 影一 } 影一 $

閣繼續看 McEniry 二千空七里的進一步說明：


: $ y={ \ frac { 一 } { \ sqrt { x } } } $

著等式的兩爿取二進位對數（$ \ textstyle \ log _ { 二 } $，即函式 $ \ textstyle f ( n )=二 ^ { n } $ 的反函式）， 有


: $ \ log _ { 二 } { ( y ) }=-{ \ frac { 一 } { 二 } } \ log _ { 二 } { ( x ) } $


: $ \ log _ { 二 } ( 一 + m _ { y } ) + e _ { y }=-{ \ frac { 一 } { 二 } } \ log _ { 二 } { ( 一 + m _ { x } ) }-{ \ frac { 一 } { 二 } } e _ { x } $

以如方法，就會當共浮點數 x 和 y 的相關指數消去，對乘方運算化為加法運算。因為 $ \ textstyle \ log _ { 二 } { ( x ) } $ 佮 $ \ textstyle \ log _ { 二 } { ( x ^ {-二分之一 } ) } $ 線性相依，所以佇遮 $ \ textstyle x $ 佮 $ \ textstyle y _ { 零 } $（即輸入值佮頭擺近來親像值）間就會當線性組合的方式建立方程式。在此 McEniry 再度引入新數 $ \ sigma $ 是咧講 $ \ textstyle \ log _ { 二 } { ( 一 + x ) } $ 佮近似值 R 間的精差：因為 $ \ textstyle 零 \ leq x < 一 $，有 $ \ textstyle \ log _ { 二 } { ( 一 + x ) } \ approx { x } $，則在此可定義 $ \ sigma $ 佮 x 的關係為著 $ \ textstyle \ log _ { 二 } { ( 一 + x ) } \ cong x + \ sigma $，這一定義就會當提供二進位對數的頭一改的精度值（此處 $ 零 \ leq \ sigma \ leq { \ tfrac { 一 } { 三 } } $；當 R 為零 x 五 f 三千七百五十九 df 時，有 $ \ textstyle \ sigma=空八空四五空四六一八七五七九一六八七 $）。 由此將 $ \ textstyle \ log _ { 二 } { ( 一 + x ) }=x + \ sigma $ 代入上式，有：


: $ m _ { y } + \ sigma + e _ { y }=-{ \ frac { 一 } { 二 } } m _ { x }-{ \ frac { 一 } { 二 } } \ sigma-{ \ frac { 一 } { 二 } } e _ { x } $

參照首段等式代入 $ M _ { x } $，$ E _ { x } $，$ B $ 佮 $ L $，有：


: $ M _ { y } + ( E _ { y }-B ) L=-{ \ frac { 三 } { 二 } } \ sigma { L }-{ \ frac { 一 } { 二 } } M _ { x }-{ \ frac { 一 } { 二 } } ( E _ { x }-B ) L $

項項整理甲：


: $ E _ { y } L + M _ { y }={ \ frac { 三 } { 二 } } ( B-\ sigma ) L-{ \ frac { 一 } { 二 } } ( E _ { x } L + M _ { x } ) $

如何咧講，對以浮點規格儲存的正浮點數 x，若是共做長整型表示示示值為 $ \ textstyle I _ { x }=E _ { x } L + M _ { x } $，由此就會當根據 x 的整數表示匯出 y（在此 $ \ textstyle y={ \ frac { 一 } { \ sqrt { x } } } $，亦即 x 的平方根倒算的頭擺近來親像值）的整數表示值，嘛即：


: $ I _ { y }=E _ { y } L + M _ { y }=R-{ \ frac { 一 } { 二 } } ( E _ { x } L + M _ { x } )=R-{ \ frac { 一 } { 二 } } I _ { x } $，其中 $ R={ \ frac { 三 } { 二 } } ( B-\ sigma ) L $ .

上尾匯出的等式 $ \ textstyle I _ { y }=R-{ \ frac { 一 } { 二 } } I _{ x } $ 即佮上節代的碼中 i=零 x 五 f 三千七百五十九 df-( i > > 一 ) ; 一逝相契，由此可見，佇咧反平方根快速演算法內底，著浮點數進行一擺移位元運算和整數減法，就會當靠地輸出一个浮點數的對應近若像值。為止喔，McEniry 只證明矣，佇常數 R 的輔助之下，可近來親像求取浮點數的平方根倒算，但是猶未確定講代碼內底 R 值的選取方法。

關於做一擺移位佮減法操作用浮點數的指數被-二除的原理，Chris Lomont 的論文中亦有一个相對簡單的解說：以 $ \ textstyle 一孵=十 ^ { 四 } $ 做例，共伊指數除了-字會用得 $ \ textstyle 一孵 ^ {-二分之一 }=十 ^ { 鋪二 }=一百分之一 $；因為浮點表示的指數有進行過偏移處理，所以指數的真實值 e 應為 $ \ textstyle e=E 鋪百二七喔 $，所以會當知影講除法來操作的實際結果為著 $ \ textstyle-e / 二 + 一百二十七喔 $，這時用 R（在此即為「魔術數字」零 x 五 f 三千七百五十九 df）減之即可使指數的最低有效數位轉入有效數字域，了後重新轉換做浮點數的時，就通得著相當精確的平方根倒算近來若值。佇遮對常數 R 的選取亦有所講究，若會當選一个好的 R 值，就會當減少對指數進行除法佮對有效數字域來進行移位的時可能產生的錯誤。是因為這標準，零 xbe 即是上合適的 R 值，啊若無 xbe 正徙一位即可得著零 x 五 f，這是魔術數字 R 的頭一个位元組。

===精確度===

如何咧講，反平方根快速演算法所得的近似值驚人的精確，正圖抑是展示以上述代碼計算（用反平方根快速演算法計算後閣進行一擺牛頓法迵天代）所得近若值的誤差：做輸入去零交零一時，以 C 語言標準庫函式計算會當十影零，而且 InvSqrt ( ) 得為著九石九八二五八二二，其間誤差為空普遍空一七四七九，相對誤差是為百分之空石一七五，閣做輸入閣較大的數值的時陣，絕對誤差不多下降，相對的精差嘛一直控制佇咧一定的範圍內底。

===牛頓法提懸精度===

咧進行矣若上的整數操作了後，範例程式再度將被轉做長正型的浮點數轉做浮點數（對應 x=\ * ( float \ * ) & i ;）， 並對其進行一改浮點運算操作（對應 x=x \ * ( 一垺五 f-xhalf \ * x \ * x ) ;）， 遮的浮點運算操作就是嘿其進行一擺牛翻頭天頂，就這个例說明：


: $ y={ \ frac { 一 } { \ sqrt { x } } } $ 所求的是講 y 的平方根倒算，以伊起造以 y 為自變數的函式，有 $ f ( y )={ \ frac { 一 } { y ^ { 二 } } }-x=零 $，


: 共其代入去牛頓法的通用公式 $ y _ { n + 一 }=y _ { n }-{ \ frac { f ( y _ { n } ) } { f'( y _ { n } ) } } $（其中 $ \ , y _ { n } $ 為頭一擺近若值），


: 有 $ y _ { n + 一 }={ \ frac { y _ { n } ( 三-xy _ { n } ^ { 二 } ) } { 二 } } $，其中 $ f ( y )={ \ frac { 一 } { y ^ { 二 } } }-x $，$ f'( y )={ \ frac { 鋪二 } { y ^ { 三 } } } $。


: 整理有 $ \ , y _ { n + 一 }={ \ frac { y _ { n } ( 三-xy _ { n } ^ { 二 } ) } { 二 } }=y _ { n } ( 一垺五-{ \ frac { xy _ { n } ^ { 二 } } { 二 } } ) $，對應的代碼即為 y=y \ * ( threehalfs-xhalf \ * y \ * y ) ;。

佇以上的一節的整數操作產生頭擺近若值了後，程式會共頭擺近似值作為參數送入函式落尾兩句來進行精化處理，代碼內底迵天（迵過迵天的輸出（對應公式中的 $ y _ { n + 一 } $）作為迵天天代的輸入）正正是為著進一步提懸結果的精度，但是因為雷神之槌 III ia̋n-jín 的圖形算中並無需要傷懸的精度，所以代碼內底干焦進行迵天代，迵過迵天代的代碼是予人注意講。

==歷史佮考究==

《 雷神之槌 III》的代碼到甲 QuakeCon 兩千空五才正式放出講，毋過早佇二空空二年（抑是二空空三年）時，反平方根快速演算法的代碼就已經出現佇 Usenet 佮其他的論壇上矣。頭先人咧臆講是卡馬克寫這段代碼，但是伊在回覆講問伊的郵件的時陣是毋是定著這个觀點，並且咧臆可能是進前捌替 id Software 上好化雷神的摃槌仔的資深組譯程式設計師 Terje Mathisen 去寫這个代碼落去；啊若佇咧 Mathisen 的郵件內底，伊表示，佇一九九空年代初，伊干焦做過類似的實作，確切來講這段代碼亦非伊所作。這馬知影的上早實作是由 Gary Tarolli 佇咧 SGI Indigo 中實作的，但伊亦坦承伊干焦對常數 R 的取值做著一定的改進，實際上伊嘛毋是作者。咧向以發明 MATLAB 啊若有名的 Cleve Moler 查證後，Rys Sommefeldt 則認為原始的演算法是 Ardent Computer 公司的 Greg Walsh 所發明，但伊嘛無任何決定性的證據會當證明這點。

毋但該演算法的原作者不明，人嘛猶無法度確定講當初選擇這个「魔術數字」的方法。Chris Lomont 捌做一个研究：伊推算出一个函式以討論此速演算法的精差，並揣出著使誤差上細的上好 R 值零 x 五 f 三更七千六百四十二 f（佮代碼中使用的零 x 五 f 三千七百五十九 df 相當接近）， 但是以代入去演算法計算閣進行一擺牛做迵天了後，所得近若值之精度猶原略低代入散 x 五 f 三千七百五十九 df 的結果；所以 Lomont 將目標改為走揣佇咧進行一－二次牛做迵天以後會用得到上大的精度 R 值，佇咧暴力搜揣了後會出最佳 R 大部份攏 x 五 f 三百七十五 a 八十六，用這值代入演算法並進行牛做迵天，所得的結果攏比代入原始值（零 x 五 f 三千七百五十九 df）閣較精確，所以伊的論文尾仔講「若可能我誠想欲問原作者，此速演算法是以數學推導猶是以重複試毋著的方式求得？」。 佇論文中，Lomont 亦指出，六十四位箍的 IEEE 七仔五十四浮點數（即雙精度類型）所對應的魔術數字是零 x 五 fe 六 ec 八十五 e 七 de 三十 da，但後來的研究表明，代入零 x 五 fe 六 eb 五十 c 七 aa 十九 f 九的結果精確度閣較懸（McEniry 得出的結果則是零 x 五 FE 六 EB 五十 C 七 B 五百三十七 AA，精度介紹兩个之間，英文 wiki 予出的精度閣較懸的值是零 x 五 FE 六 EB 五十 C 七 B 五百三十七 A 九）。 佇咧 Charles McEniry 的論文中，伊使用一種的類似 Lomont 毋過閣較複雜的方法來最佳化 R 值：伊上開始使用散赤攑步來搜揣，所得結果佮 Lomont 相仝；了後伊試驗用帶權二分法走揣最佳值，所得結果恰是代碼中所使用的魔術數字零 x 五 f 三千七百五十九 df，所以，McEniry 認為講，這一常數上蓋起初凡勢是以「佇可容忍誤差範圍內使用二分法」的方式求得。

==注釋==

==參考==

===註跤===

===參考文獻===

==延伸閱讀==

* Kushner , David . The wizardry of Id . IEEE Spectrum . Aug 兩千空二 ,'''三十九'''( 八 ) : 四十二–四十七 . doi : 十席一一空九 / MSPEC . 二千空二四一空二一九四三 .
* 顧森（Matrix 六十七）. 鋪排一爿：神祕的零 x 五 f 三千七百五十九 df 不可思學的 Quake III 源孵 . 二千空七孵十一孵二十四 [二千空一十五孵八孵一] .（原始內容存檔佇兩千空一十五五分八配四）.

==外部連結==

* Origin of Quake 三's Fast InvSqrt ( ) , Beyond 三 D . com
* Quake III Arena source code , id Software
* Margolin , Tomer . Magical Square Root Implementation In Quake III . CodeMaestro . The Coding Experience . 二千空五孵八鋪二十七 .（原始內容存檔佇兩千空一十二抹四四）.

[[分類: 待校正]]