- Quake 3で有名になった fast inverse square root は、
1 / sqrt(x) を float のビット再解釈と Newton-Raphson 補正によって高速に近似する、当時の性能上の解決策である
- 核心は、IEEE-754 32ビット float の整数ビットパターンを スケールおよびシフトされた
log2(x) 近似 のように扱える点にある
0x5f3759df - (i >> 1) は log2(x^-0.5) = -0.5 * log2(x) を整数シフトと減算に移した形であり、マジック定数は 3/2 * 2^23 * (127 - σ) から導かれる- その後
y = y * (1.5 - 0.5x * y * y) を1回適用して Newton-Raphson 補正 を行い、Quake のコードでは2回目の反復はコメントアウトされている
- 1999年には照明計算や 3D ベクトル正規化のために逆平方根が毎秒数百〜数千回必要だったが、現代のハードウェアでは専用の浮動小数点処理のおかげで同じトリックの実用性は低下している
Quake のコードがしていること
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long*)&y;
i = 0x5f3759df - ( i >> 1 );
y = *(float*)&i;
y = y * ( threehalfs - ( x2 * y * y ) );
return y;
}
- この関数は
number に対する 逆平方根 1 / sqrt(number) の近似値を計算する
- 最も有名な部分は、float 値を
long として解釈した後に 0x5f3759df - (i >> 1) を実行するビット操作である
- Quake 3 が発売された 1999 年には逆平方根は遅く高価な演算であり、照明方程式や正規化が必要な 3D ベクトル計算で毎秒数百〜数千回必要とされていた
- 現代のハードウェアでは、この種の計算は CPU で実行されないことも多く、CPU 上で実行される場合でも発達した専用浮動小数点ハードウェアのおかげで高速である
IEEE-754 32ビット float 表現
- 32ビット float は3つの部分で構成される
- Sign: 1ビットで、正負を表す
- Exponent: 8ビットで、値が属する範囲を決める
- Mantissa: 23ビットで、その範囲内での位置を線形に表す
- 一般的な値は次の形で解釈される
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
B = 127 は biased exponent に使われるバイアス値で、実際の指数は e = E - B である- mantissa は単純に
m を掛けるのではなく 1 + m の形で使われる
m = 0 なら 2^em = 1 に近づくと、次の指数範囲である 2^(e+1) の直前までを表現する
- 指数がすべて 0 の場合は sub-normal 数であり、式が異なる
N = (-1)^S * 2^-126 * m
- sub-normal は 0 および 0 に非常に近い小さな数を表現するために必要である
- 指数がすべて 1 の場合は特殊値として扱われる
E = 255, M = 0 なら Infinity または -InfinityM != 0 なら NaN
float のビットを整数として見ると現れる対数関係
- float の内部表現を 32ビット整数として見ると、次の式で表せる
I_x = 2^31 S + 2^23 E + M
- 逆平方根は正の入力を対象とするため、
S = 0 とすると式は単純になる
L = 2^23
I_x = L E + M
- 同じ exponent 範囲内では mantissa が位置を線形に表すが、exponent が大きくなるほど同じ数の mantissa ステップがより広い数直線区間をカバーする
E = 127、つまり e = 0 はおおよそ [1, 2) の範囲E = 128、つまり e = 1 はおおよそ [2, 4) の範囲- どちらの範囲も mantissa ステップ数は同じだが、後者の範囲は2倍広い
- この構造のため、float の生ビットパターンを整数として見ると 対数的な関係 が現れる
生ビットは log2(x) の近似
- float のビットパターンを整数
I_x として解釈すると、log2(x) に対する区分線形近似のように見ることができる
- この関係は次の近似式で表される
log2(x) ≈ I_x / L - B
- 生ビット整数を mantissa サイズ
L = 2^23 で割り、exponent bias B = 127 を引くと、log2(x) に近い値になる
- mantissa 区間内の対数は線形近似として扱う
log2(1 + x) ≈ x + σ
σ は近似を調整する チューニングパラメータ であり、x は [0, 1] の範囲で exponent 区間内の位置を表す
逆平方根を対数恒等式に変換する
y = 1 / sqrt(x)
y = x^-0.5
- 対数恒等式を適用すると、逆平方根の計算は次の関係になる
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
- float のビットが
log2(x) 近似のように動作することを利用すると、x の整数ビット表現 I_x から y の整数ビット表現 I_y を直接近似できる
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
- この式が Quake のコードにある核心の1行につながる
i = 0x5f3759df - ( i >> 1 );
i >> 1 は整数ビットを右に1ビットシフトして 1/2 を掛ける役割を果たす- 手前の定数
0x5f3759df は 1.5 * L * (B - σ) に相当する
0x5f3759df 定数の正体
1.5 * 2^23 * 127 = 1598029824
- この値の16進表現は
0x5f400000 である
- Quake の実際の定数
0x5f3759df とは 566817 だけ差がある
- この差から Quake のコードに対応する
σ の値を計算すると次のようになる
σ = 377878 / 2^23
σ = 0.04504656
int32_t compute_magic(void) {
double sigma = 0.0450465;
double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
int32_t i = expression;
return i;
}
// -> 0x5f3759df
- ここでは
double を使用しており、整数変換はビット再解釈ではなく通常のキャストである
- この
σ の値は近似を最適化するために選ばれたが、実際の最適値ではなく、誰が作ったのかもはっきりしていない
単なるハックではない理由
0x5f3759df - (i >> 1) は、float の生ビットが対数近似であるという事実を利用して 逆平方根の初期値 を作る式である- 複雑な数学的関係に基づいているが、実行段階ではシフトと減算のような高速な演算だけを使う
- 当時は高価な演算を毎秒数千回処理する必要があったため、この方式はハードウェア制約に適した工学的設計だった
- ただし、このアルゴリズムは normal float でのみ動作する
- sub-normal 値では
log2(1 + x) ≈ x + σ という近似の前提が成り立たない
- sub-normal では実際には
0 + x に近い形になるため、近似が破綻する
Newton-Raphson 補正で誤差を減らす
- ビット操作で得た初期値はかなり良いが、測定可能な誤差は残る
- 次の行が近似を大きく改善する
y = y * ( threehalfs - ( x2 * y * y ) );
- この行は Newton-Raphson method を適用した形である
- 逆平方根の問題を Newton method に合わせるため、次の関数の根を求める問題に変換する
f(y) = 1 / y^2 - x = 0
- Newton method は現在の近似値
y_n から、より良い近似値 y_(n+1) を次のように作る
y_(n+1) = y_n - f(y_n) / f'(y_n)
f(y) = y^-2 - x の導関数は次のとおりである
f'(y) = -2y^-3 = -2 / y^3
除算のない Newton 補正式
- Newton の式をそのまま使うと、浮動小数点除算が何度も入る
- このアルゴリズムが高速な理由の1つは 浮動小数点除算の回避 にある
- 代数的に整理すると、除算なしで乗算だけを使う形になる
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
- Quake のコードでは
x2 = number * 0.5F で 0.5x をあらかじめ計算し、これを次の行で使っている
y = y * ( threehalfs - ( x2 * y * y ) );
- この1回の反復の後、最大絶対誤差は 0.175% であり、多くの場合はこれより小さい
- 元のコードには2回目の Newton 反復もあるが、コメントアウトされている
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
起源と関連アルゴリズム
- このアルゴリズムは John Carmack が発明したものではなく、正確な起源は 100% 確実ではない
- 関連情報として Beyond3D の記事がリンクされている: The truth is the exact origin is not 100% certain
- Chris Lomont は、対数近似段階で 最適な sigma 値 を探す論文を書いている: InvSqrt.pdf
- CORDIC は浮動小数点なしで加算とビットシフトだけを使って sine と cosine を計算するアルゴリズムであり、fast inverse square root とは詳細な方式が大きく異なる
- 両アルゴリズムに共通するのは、数学的観察を当時のハードウェア制約に合わせて効率的に適用している点である
2件のコメント
忘れた頃にまた上がってくる不思議なコード…w
Hacker News のコメント
1999年以降に作られたコンピュータなら、たいてい SSE 命令セットをサポートしており、その中には逆平方根を一度に4つ、より高速に計算する
_mm_rsqrt_psがある: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...それでも、ここで扱われている手法がまだ完全に無意味になったわけではない。float/int 変換は速いが、
rsqrt、sqrt、pow、log命令を持たないハードウェアも今なお存在し、こうした演算はこのトリックで近似できるsqrt(x)をx * 1/sqrt(x)として計算するほうが速いGPU 命令セット、ARM、RISC-V、AVR、PIC、8051、FPGA などでは近似逆平方根演算が組み込まれていることも多いが、おそらくこうしたアルゴリズムで実装されている可能性が高い
記事に少し難癖をつけるなら、こうした計算が最近の CPU では行われないかのような説明は正しくない。ゲームや浮動小数点演算の多いアプリが、すべての浮動小数点演算を GPU に渡したがっているというのはよくある誤解だ
実際に GPU に渡すのが妥当なのは、大きく均一な処理だけだ。あるオブジェクトを別のオブジェクトの方に向けるための回転行列を構成するような 一回限りのベクトル正規化なら、CPU 側に残したほうが速い。GPU への転送時間を除いて見ても、単一の浮動小数点演算は CPU のほうが速い。GPU は通常クロックが低く、高い FLOP 数を並列性で達成しているからだ
MMIX 実装を書いてみた。元の入力値が
2^-1021より大きいという仮定を置いている興味があるなら、Wikipedia にもこの関数と歴史についての良い説明がある: https://en.wikipedia.org/wiki/Fast_inverse_square_root
こういうものをいくつか集めている: https://github.com/ncruces/fastmath/blob/main/fast.go
関連する StackOverflow の投稿もある: https://stackoverflow.com/questions/32042673/optimized-low-a...
fastmathパッケージの ベンチマークも見てみたい難癖をつける時間だ。float の式にタイプミスがあり、
-1^Sではなく(-1)^Sであるべきだ。前者は常に-1になる生のビットパターンを解釈することが対数の区分的線形近似だという説明も正確ではない。青いグラフのデータ点同士を結ぶ線は実際には存在せず、ビットを半分だけ 1 に設定することはできない。むしろ対数の 離散版に近く、実際に存在するデータ点、つまり赤い線と青い線が交わる点は、文字どおりスケール・移動された対数と同じだ。それ以外は良い記事だ
[010000, 010111]の区間は 2、2.25、2.5、2.75、3、3.25、3.5、3.75 を含むしかし、これらの数の底 2 の対数が示唆する仮数はそれぞれ
.0000000、.0010101、.0101001、.0111010、.1001010、.1011001、.1100111、.1110100であり、最初を除けば float の001、010などとは一致しない。[2,4)区間の float は線形間隔だが、対応する対数はそうではないため、記事の言うように float は対数の 区分的線形近似と見なせるグラフ全体であれば、区分的線形パターンの中に
2^32個の選択肢があるはずだが、元記事が描いているのはそうした全体グラフではない。記事が 32 ビット整数と IEEE-754 32 ビット float 演算を扱っている以上、説明で「離散」を省略しても問題ないと思う興味深い概念をたくさん説明している良い記事だが、あるセクションの 代数展開が驚くほどひどい
「最初の形からこの形へ至る正確なステップは多いが、完全性のためにすべて含めた」の後の展開には不要なステップが多く、互いに打ち消し合う符号の誤りもいくつもある。特に2行目から3行目に行くところで、負号が正しく分配されていない。2行目以降なら
y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2)から始めて、y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n)までずっと短く進められるし、中間ステップも正しい。代数を理解している人には自明なステップばかりだと思う有名なコード片の マジックナンバーは最適な定数ではない。別の定数を使えば、相対誤差をたぶん 0.5% ほどさらに減らせる
当時は絶対的な最適値を見つけるのが難しかったかもしれないが、今では比較的簡単だ。私もかつてこの沼にはまり、
(1/x^2)と(1/x)の最適なマジックナンバーを探す Jupyter Notebook を持っているこの記事で一番興味深かったのは “How Java's Floating-Point Hurts Everyone Everywhere” へのリンクだった: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
著者は “Old Man of Floating-Point” としても知られる William Kahan だ: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
JAVAhurtの PDF を読み始めたところ、組版がひどい。単語間隔を過度に、しかも不揃いに広げる TeX パッケージを使ったような感じで、別の文書を OCR したあと余分な空白が入ったようにも見える等幅フォント部分にも妙な追加スペースがある。読むのに集中するのが本当に難しく、実際には違うと分かっているが、ほとんど科学オタクのマニフェストのように感じられた
以前見たこの動画が本当に良かった: https://www.youtube.com/watch?v=p8u_k2LIZyo