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

'''反平方根快速演算法'''（英語：Fast Inverse Square Root，乜輒常以「Fast InvSqrt ( )」或者使用个十六進位常數'''零 x 五 f 三千七百五十九 df'''代表）係用在快速計算 $ \ textstyle x ^ { - 二分之一 } $（就係 $ \ textstyle x $ 个平方根个倒算，在這 $ \ textstyle x $ 需要符合資格 IEEE七百五十四標準格式个三十二個銀浮點數）个一種演算法。這一演算法个優勢係在減少咧求平方根倒數个時節浮點運算操作帶來个當大个運算耗費，在電腦圖學个領域，若講要求取照明還有投影个動角度摎反射效果，就常需要計算平方根倒數。

這擺演算法首先接收一個三十二個元帶符浮點數，過仔呢就會分佢做一個三十二個銀整數來看，佢就用向右片進行一擺邏輯移位个方式摎這往年个方式取半，還用在浮點數規格代表 √ 二十一二十七近似值个十六進位「奇術數字」零 x 五 f 三千七百五十九 df 減之，恁樣就做得拿著輸入个浮點數个平方根橫忒幾下擺第一擺近似值；以後重新摎佢作為浮點數，以牛頓法反覆迭代，以求出更加精確个近似值，一直到求出符合精確度要求个近似值。在計算浮點數个平方根倒數个共精度个近似值時，這擺演算法比直接使用浮點數除法愛遽四倍。

這擺演算法最早分人認為係由約翰 ・ 卡馬克在九零年代前期在 SGIIndigo 个開發當中使用，過後係一九九年在嗬《雷神之鎚 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 } } $，

做得看著標準化向量時，對向量分量計算平方根橫忒係定著愛个，故所，對平方根項算演算法最佳化對計算當在規化向量乜大有737益。

為著加速圖像處理單元个計算，《 雷神之鎚 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=一百二十四八點一百二十七=三步八步九百五十二 $，係做得表示个浮點數係 $ \ textstyle x=( 一 + 零嗒二五零 ) \ cdot 二 ^ { 三步八步九百五十二 }=零嗒一五六二五 $。

比將講嗬，一個有符號正在整數在二補數系統裡背个表示第一位係零，後背个各位用來表示。將浮點數取別名儲存為整數个時節，這個整數值就係 $ \ 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 176一百二十七 $，所以做得知除法操作个實際結果係 $ \ 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 引擎个圖形計算當中並無需要昶高个精度，所以代碼肚項單淨進行了一擺迭代，第二擺迭代个代碼係分人注釋。==歷史同考究==

《 雷神之鎚 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 . 二千空五千八百五十二十七 .（頭擺內容存檔係在兩千空一十二pa24四千五百空四）.

[[分類: 待校正]]