高速逆平方根

例えば、その数字は${\displaystyle x=0.15625}$計算に使用できます$\textstyle {\frac {1}{\sqrt {x}}}\approx 2.52982}$アルゴリズムの最初のステップを以下に示します。

0011_1110_0010_0000_0000_0000_0000_0000 xとiの両方のビットパターン 0001_1111_0001_0000_0000_0000_0000_0000 右に1つシフトします: (i > > 1) 0101_1111_0011_0111_0101_1001_1101_1111 魔法の数字 0x5F3759DF 0100_0000_0010_0111_0101_1001_1101_1111 0x5F3759DF の結果 - (i > > 1)

IEEE 32ビット表現として解釈します。

0_01111100_01000000000000000000000 1.25 × 2 −3 0_00111110_00100000000000000000000 1.125 × 2 −65 0_10111110_01101110101100111011111 1.432430... × 2 63 0_10000000_01001110101100111011111 1.307430... × 2 1

この最後のビットパターンを浮動小数点数として再解釈すると、近似値が得られます。${\displaystyle y=2.61486}$これは約3.4%の誤差があります。ニュートン法を1回繰り返した後の最終結果は${\displaystyle y=2.52549}$誤差はわずか0.17%です。

C 標準によれば、浮動小数点値をキャストしてそのポインタを逆参照することで整数として再解釈すると、指定されたアーキテクチャで整数と浮動小数点のサイズが一致しない場合に未定義の動作を引き起こす可能性があります。 ^{[ 23 ]これは、C の共用体[ 24 ]}やC++20の^{[ 25 ]}などの代替の型変換技術を使用することで回避できます。std::bit_cast

このアルゴリズムは単精度浮動小数点数のビットレベル表現に大きく依存しているため、ここではその表現について簡単に概説します。ゼロ以外の実数をエンコードするには${\displaystyle x}$単精度浮動小数点数として、最初のステップは${\displaystyle x}$正規化されたバイナリ数として: ^{[ 26 ]}

${\displaystyle {\begin{aligned}x&=\pm 1.b_{1}b_{2}b_{3}\ldots \times 2^{e_{x}}\end{aligned}}}$

指数${\textstyle e_{x}}$は整数であり、${\textstyle 1.b_{1}b_{2}b_{3}\ldots }$は仮数の二進数表現です。仮数の小数点前の1ビットは常に1なので、保存する必要はありません。この式は次のように書き換えることができます。

${\displaystyle {\begin{aligned}x&=(-1)^{S_{x}}\cdot 2^{e_{x}}(1+m_{x})\end{aligned}}}$

どこ${\textstyle m_{x}}$手段${\textstyle 0.b_{1}b_{2}b_{3}\ldots }$、 それで${\textstyle m_{x}\in [0,1)}$この形式から、3 つの符号なし整数が計算されます。^{[ 27 ]}

• ${\textstyle S_{x}}$「符号ビット」は${\displaystyle 0}$もし${\displaystyle x}$ポジティブで${\displaystyle 1}$負またはゼロ（1ビット）
• ${\textstyle E_{x}=e_{x}+B}$は「バイアス指数」であり、${\displaystyle B=127}$「指数バイアス」^{[注2 ]}（8ビット）
• ${\textstyle M_{x}=m_{x}\times L}$、 どこ${\textstyle L=2^{23}}$^{[注3 ]} (23ビット)

したがって：${\textstyle e_{x}=E_{x}-B}$そして${\textstyle m_{x}={\frac {M_{x}}{L}}}$。

これらのフィールドは、左から右に32ビットコンテナにパックされます。^{[ 28 ]}

例として、もう一度数字を考えてみましょう。${\displaystyle x=0.15625=0.00101_{2}}$正規化${\displaystyle x}$収量：

${\textstyle x=(-1)^{0}\cdot 2^{-3}(1+0.25)=+2^{-3}(1+0.25)}$

したがって、符号なし整数フィールドは次の3つです。

• ${\displaystyle S=0}$
• ${\displaystyle E=-3+127=124=0111\ 1100_{2}}$
• ${\displaystyle M=0.25\times 2^{23}=2\ 097\ 152=0010\ 0000\ 0000\ 0000\ 0000\ 0000_{2}}$

これらのフィールドは、下図に示すように配置されます。

この数値はバイナリで次のように表されます。${\displaystyle I_{x}=S_{x}\cdot 2^{31}+E_{x}L+M_{x}}$

また、このアルゴリズムは実数で動作するため、${\textstyle {\sqrt {x}}}$のみ定義されています${\textstyle x\geq 0}$コードは、${\textstyle x\geq 0}$そして${\displaystyle S_{x}=0}$。

平方根を計算するために与えられた数値は、次のように書き換えることができます。

${\textstyle I_{x}=E_{x}L+M_{x}}$

${\textstyle x=(-1)^{0}\cdot 2^{e_{x}}(1+m_{x})=+2^{e_{x}}(1+m_{x})}$

もし${\textstyle {\frac {1}{\sqrt {x}}}}$コンピュータや電卓を使わずに計算する場合、対数表と恒等式が役立つだろう。${\textstyle \log _{b}\left({\frac {1}{\sqrt {x}}}\right)=\log _{b}\left(x^{-{\frac {1}{2}}}\right)=-{\frac {1}{2}}\log _{b}(x)}$これはすべての基数に有効です${\displaystyle b}$高速逆平方根は、この恒等式と、float32 を整数にエイリアスするとその対数の近似値が得られるという事実に基づいています。その方法は次のとおりです。

もし${\displaystyle x}$は正の正規数です。

${\displaystyle x=2^{e_{x}}(1+m_{x})}$

それから

${\displaystyle \log _{2}(x)=e_{x}+\log _{2}(1+m_{x})}$

そしてそれ以来${\displaystyle m_{x}\in [0,1)}$右辺の対数は^{[ 29 ]で近似できる。}

${\displaystyle \log _{2}(1+m_{x})\approx m_{x}+\sigma }$

どこ${\displaystyle \sigma }$は近似を調整するために使用される自由パラメータです。たとえば、${\displaystyle \sigma =0}$区間の両端で正確な結果が得られる一方、${\textstyle \sigma ={\frac {1}{2}}-{\frac {1+\ln(\ln(2))}{2\ln(2)}}\approx 0.0430357}$これは最適な近似値（誤差の一様ノルムの意味で最良の値）をもたらします。ただし、この値は後続のステップを考慮しないため、アルゴリズムでは使用されません。

したがって、近似値が存在する

${\displaystyle \log _{2}(x)\approx e_{x}+m_{x}+\sigma .}$

浮動小数点ビットパターンの解釈${\displaystyle x}$整数として${\displaystyle I_{x}}$収量^{[注4 ]}

${\displaystyle {\begin{aligned}I_{x}&=E_{x}L+M_{x}\\&=L(e_{x}+B+m_{x})\\&=L(e_{x}+m_{x}+\sigma +B-\sigma )\\&\approx L\log _{2}(x)+L(B-\sigma ).\end{aligned}}}$

すると、${\displaystyle I_{x}}$は、スケーリングおよびシフトされた区分的線形近似である。${\displaystyle \log _{2}(x)}$右の図に示すように。言い換えれば、${\displaystyle \log _{2}(x)}$近似値は

${\displaystyle \log _{2}(x)\approx {\frac {I_{x}}{L}}-(B-\sigma ).}$

計算${\textstyle y={\frac {1}{\sqrt {x}}}}$アイデンティティに基づいています

${\displaystyle \log _{2}(y)=-{\tfrac {1}{2}}\log _{2}(x)}$

上記の対数の近似を両方に適用すると${\displaystyle x}$そして${\displaystyle y}$上記の式から、次の式が得られます。

${\displaystyle {\frac {I_{y}}{L}}-(B-\sigma )\approx -{\frac {1}{2}}\left({\frac {I_{x}}{L}}-(B-\sigma )\right)}$

したがって、${\displaystyle I_{y}}$は：

${\displaystyle I_{y}\approx {\tfrac {3}{2}}L(B-\sigma )-{\tfrac {1}{2}}I_{x}}$

コードには次のように書かれています。

i = 0x5f3759df - ( i >> 1 );

上記の最初の項は魔法の数字です

${\displaystyle {\tfrac {3}{2}}L(B-\sigma )={\mathtt {0x5f3759df}}}$

そこから推測できることは${\displaystyle \sigma \approx 0.0450466}$2番目の項は、${\displaystyle {\frac {1}{2}}I_{x}}$は、ビットをシフトすることによって計算されます。${\displaystyle I_{x}}$右に1つ位置。^{[ 30 ]}

ニュートン法による根探索法を0回、1回、2回、3回、4回繰り返した場合の、直接計算と高速逆平方根法との相対誤差。なお、倍精度浮動小数点数を使用しており、2つの倍精度浮動小数点数間の最小表現可能差は4回の繰り返し後に達成される。

番号${\displaystyle y={\tfrac {1}{\sqrt {x}}}}$は方程式の解である${\displaystyle {\tfrac {1}{y^{2}}}-x=0}$前のステップで得られた近似値は、関数の零点を見つける方法である根探索法を使用して改良できます。このアルゴリズムはニュートン法を使用します。近似値がある場合、${\displaystyle y_{n}}$のために${\displaystyle y}$すると、より良い近似値が得られます${\displaystyle y_{n+1}}$計算するには、${\displaystyle y_{n}-{\tfrac {f(y_{n})}{f'(y_{n})}}}$、 どこ${\displaystyle f'(y_{n})}$は、​${\displaystyle f(y)}$で${\displaystyle y_{n}}$[ ^{31 ]}方程式に^適用${\displaystyle {\tfrac {1}{y^{2}}}-x=0}$ニュートン法では

${\displaystyle {\begin{aligned}f(y)&={\frac {1}{y^{2}}}-x\\f'(y)&=-{\frac {2}{y^{3}}}\\y_{n+1}&=y_{n}-{\frac {f(y_{n})}{f'(y_{n})}}\\&=y_{n}+{\frac {y_{n}^{3}}{2}}\left({\frac {1}{y_{n}^{2}}}-x\right)\\&=y_{n}\left({\frac {3}{2}}-{\frac {x}{2}}y_{n}^{2}\right),\end{aligned}}}$

これはコードでは次のように書かれています。y=y*(threehalfs-(x2*y*y));

この手順を繰り返すことで、関数の出力（${\displaystyle y_{n+1}}$）を次の反復の入力として、アルゴリズムは${\displaystyle y}$逆平方根に収束する。^{[ 32 ]} Quake IIIエンジンの目的のために、1 つの反復のみが使用されました。2 番目の反復はコード内に残っていましたが、コメントアウトされました。^{[ 21 ]}

魔法数の正確な値がどのように決定されたかは正確には分かっていません。クリス・ロモントは、魔法数を選択することで近似誤差を最小化する関数を開発しました。${\displaystyle R}$範囲にわたって。彼はまず、線形近似ステップの最適定数を0x5F37642Fと計算しました。これは0x5F3759DFに近い値ですが、この新しい定数はニュートン法の 1 回反復後にわずかに精度が低下しました。^{[ 33 ]}ロモントは次に、ニュートン法の 1 回および 2 回反復後でも最適な定数を探し、すべての反復段階で元の値よりも精度が高い0x5F375A86 を見つけました。 ^{[ 33 ]}彼は最後に、元の定数の正確な値が導出によって選択されたのか、試行錯誤によって選択されたのかを尋ねました。^{[ 34 ]} ロモントは、64 ビット IEEE754 サイズ型の double のマジック ナンバーは0x5FE6EC85E7DE30DAであると述べましたが、後に Matthew Robertson によって正確には0x5FE6EB50C7B537A9であることが示されました。^{[ 35 ]}

Jan Kadlec は、単一のニュートン法反復における定数も調整することで相対誤差をさらに 2.7 倍に減らし、 ^{[ 36 ]}徹底的な探索の末に

conv.i = 0x5F1FFFF9 - ( conv.i >> 1 ) ; conv.f * = 0.703952253f * ( 2.38924456f - x * conv.f * conv.f ) ; return conv.f ;​​​​​​​​​

単精度浮動小数点数については、マジックナンバーを決定するための完全な数学的解析が利用可能になった。^{[ 37 ] [ 38 ]}