y_{1} = y_{0} (\frac{3}{2} - \frac{x}{2} y_{0}^{2})

計算の裏側

Quake IIIの高速逆平方根とニュートン法

Fast inverse square root and Newton's method in Quake III

定数から半分を引くだけで、1/√xの見当が誤差3.44%以内でつく。

ビットの並びを整数として読むとlog₂xの近似になる。0x5f3759dfから半分を引く見当の誤差3.44%が、ニュートン法1回で0.175%に縮む仕組みを、単精度の正規化数すべてで確かめる。

種類: 計算の裏側
公式: 3 本
人物: ジョセフ・ラフソン、アイザック・ニュートン、ウィリアム・カハン、クリス・ロモン

テイラー展開対数の法則三平方の定理

場面

割り算も平方根も使わない1/√x

$x = 4$ として、 $1/ x$ を割り算も平方根も使わずに求めてみる。4を単精度の浮動小数点数にした32ビットを、そのまま整数として読むと1,082,130,432になる。右へ1ビットずらして半分の541,065,216にし、0x5f3759df、10進では1,597,463,007から引くと、1,056,397,791。この整数を、今度は浮動小数点数のビットとして読み直すと0.4831になる。

本当の値0.5より、3.4%小さいだけだ。

仕上げは掛け算と引き算だけで済む。0.4831に、1.5から $2 \times 0.483 1^{2}$ を引いた1.0332を掛けると0.49915。誤差は0.17%まで減った。

この手順は、Quake III Arenaのソースにある関数Q_rsqrtそのものである。id Softwareが2005年8月20日にGPLで公開したソースの、数学の関数を集めたファイルq_math.cに入っている。浮動小数点数のビットの並びを整数として読むと $lo g_{2} x$ の近似になるので、半分にして定数0x5f3759dfから引くだけで、誤差3.44%以内の見当がつく。そのあとのニュートン法1回で、誤差は0.18%以内に縮む。

では、なぜビットの並びが対数になるのか。そして、定数の下の桁を変えると、答えのどこが変わるのか。

考え方

長さ1のベクトルを作る

3Dのゲームで面に光を当てるとき、面が光のほうをまっすぐ向いていれば明るく、斜めになるほど暗く塗る計算がある。明るさを、面の向きと光の向きがなす角 $θ$ のコサインに比例させるやり方だ。2つの向きを長さ1のベクトルにそろえておけば、 $cos θ$ は内積、つまり成分どうしの積を3つ足すだけで求まる。

ベクトル $(x, y, z)$ の長さは三平方の定理を2回使って $x^{2} + y^{2} + z^{2}$ になり、長さを1にするには各成分をこの長さで割ればよい。ただ、割り算を3回するより、 $1/ x^{2} + y^{2} + z^{2}$ を1回求めて3回掛けるほうが手間は少ない。Quake IIIのVectorNormalizeFastという関数は、まさにこの形をしている。成分の2乗の和をQ_rsqrtに渡し、返ってきた値を3つの成分に掛ける。関数の前のコメントによれば、長さが0かどうかは調べていない。

では、精度はどれだけあれば足りるのか。明るさを0から255の256段で表すとき、ニュートン法1回のあとの誤差0.18%は、いちばん明るい255でも0.45段にあたる。表示のために丸めても、ずれは1段までに収まる。最初の見当のままなら3.44%で、8.8段ずれてしまう。

考え方

整数として読むと対数になる

単精度の浮動小数点数では、正の数を $x = (1 + m) \times 2^{e}$ と表し、32ビットのうち8ビットに指数 $e + 127$ を、下の23ビットに仮数 $m$ を入れる。 $m$ は0以上1未満の数である。ビットの並びの詳しい作りは0.1 + 0.2が0.3にならないのはなぜかで見たので、ここでは並べ方だけを使う。

この32ビットを整数として読むと、指数の部分は $2^{23}$ の位から上に、仮数の部分は下の23桁に入っている。整数の値は $2^{23} (e + 127 + m)$ になる。2の何乗かを表す $e$ に、その間の位置を0から1の直線で表す $m$ を足した数だ。

ここで対数と比べる。 $lo g_{2} x$ は $e + lo g_{2} (1 + m)$ なので、 $lo g_{2} (1 + m)$ が $m$ と等しければ、整数を $2^{23}$ で割って127を引いた値は $lo g_{2} x$ そのものになる。実際に一致するのは区間の両端だけで、その間では曲線 $lo g_{2} (1 + m)$ が直線 $m$ より少し上を通る。差がいちばん開くのは $m = 1/ ln 2 - 1$ 、約0.443のときで、0.086である。

整数として読んだビットの並びは、2の累乗ごとに折れる折れ線で、 $lo g_{2} x$ をなぞっている。

計算機は数を持った時点で、対数の近似も持っていることになる。

図

折れ線と曲線の差

x = (1 + m)·2ᵉの1つの区間で、2を底とする対数の曲線（太線）と、ビットを整数として読んだ値にあたる直線m（細い実線）を比べた。網の部分が差で、m = 0.443で最大0.086になる。点線は直線を一律にσ = 0.045だけ上げたもの。

直線を $σ$ だけ上げると、曲線は直線から、上へ最大0.041、下へ最大0.045だけ離れる。次の章で見るように、0x5f3759dfはこの $σ$ を0.0450にした定数になっている。

式

0x5f3759dfを組み立てる

正の単精度の数を $x = (1 + m) 2^{e}$ と書くと、ビットの並びを整数として読んだ値 $I_{x}$ は、指数と仮数をつないだ数になる。
$I_{x} = 2^{23} (e + 127 + m)$
$0 \leq m < 1$ で、 $lo g_{2} (1 + m)$ を直線 $m + σ$ で近似する。 $σ$ は小さな定数である。
$lo g_{2} (1 + m) \approx m + σ$
$lo g_{2} x = e + lo g_{2} (1 + m)$ なので、手順1の $e + m$ を $lo g_{2} x - σ$ に置き換えられる。
$I_{x} \approx 2^{23} (lo g_{2} x + 127 - σ)$
求めたい $y = 1/ x$ は、対数では $lo g_{2} x$ の $- 1/2$ 倍になる。
$lo g_{2} y = - \frac{1}{2} lo g_{2} x$
$y$ のビットの並びにも手順3が使えるので、 $I_{y}$ を $I_{x}$ で書ける。
$I_{y} \approx \frac{3}{2} \cdot 2^{23} (127 - σ) - \frac{1}{2} I_{x}$
第1項は定数で、第2項はビットを1つ右へずらせば求まる。 $σ = 0.0450466$ とおくと、定数は0x5f3759dfになる。
$\frac{3}{2} \cdot 2^{23} (127 - 0.0450466) \approx 1, 597, 463, 007$

$σ$ を0にしたときの定数は0x5f400000で、16進の上の2桁5fはここで決まる。残りの桁が、 $σ$ の選び方にあたる。ソースの中では、次の3行が手順1から6の計算で、4行目が次の章のニュートン法にあたる。行末に付いている短いコメントは省いた。

i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) );

引用元はQuake III Arenaのソースのcode/game/q_math.cで、id SoftwareがGPL第2版以降で公開している。x2は $x /2$ 、threehalfsは1.5を表す。

式

割り算のないニュートン法

関数 $f$ の根の近くの点 $y_{n}$ で、 $f$ をテイラー展開の1次までで近似し、その直線が0になる点を次の $y_{n + 1}$ にする。これがニュートン法である。
$f (y_{n}) + f^{'} (y_{n}) (y_{n + 1} - y_{n}) = 0 \Rightarrow y_{n + 1} = y_{n} - \frac{f ( y _{n} )}{f ^{'} ( y _{n} )}$
$1/ x$ を根に持つ関数として $f (y) = 1/ y^{2} - x$ を選ぶ。微分すると次のとおり。
$f (y) = \frac{1}{y ^{2}} - x, f^{'} (y) = - \frac{2}{y ^{3}}$
手順1に入れて整理すると、割り算が消える。
$y_{n + 1} = y_{n} + \frac{y _{n} - x y _{n}^{3}}{2} = y_{n} (\frac{3}{2} - \frac{x}{2} y_{n}^{2})$
$y_{n}$ の相対誤差を $δ$ とし、 $y_{n} = (1 + δ) / x$ とおく。 $x y_{n}^{2} = (1 + δ)^{2}$ なので、次の値は次の形になる。
$y_{n + 1} = \frac{( 1 + δ ) ( 3 - ( 1 + δ ) ^{2} )}{2 x} = \frac{1 - \frac{3}{2} δ ^{2} - \frac{1}{2} δ ^{3}}{x}$
1回あとの相対誤差 $δ^{'}$ は $δ$ の2乗に比例する。 $δ > - 3$ なら0以下なので、式の上では、1回かけた値が本当の値を超えることはない。
$δ^{'} = - \frac{δ ^{2} ( 3 + δ )}{2}$
0x5f3759dfの見当がいちばん大きく外れる $δ = - 0.034376$ を入れると、1回目で約 $- 0.001752$ 、2回目で約 $- 0.0000046$ になる。
$- 0.034376 \to - 0.001752 \to - 0.0000046$

x = 1のとき、f(y) = 1/y² − xのグラフに、y₀ = 0.7での接線を引くと、横軸と交わる点がy₁になる。同じことをくり返すと根1/√xに近づき、相対誤差は−30%、−12%、−2.1%と縮む。

手順2で $f (y) = y^{2} - 1/ x$ や $f (y) = x y^{2} - 1$ を選ぶと、手順3の式に $x$ や $y$ で割る計算が残る。 $1/ y^{2} - x$ を選んだのは、掛け算と引き算だけにするためだ。

手順5の式はニュートン法の性質をそのまま表している。根の近くから始めれば、1回ごとに誤差はほぼ「2乗の1.5倍」になる。最初の見当が3.4%なら、式の上では1回で0.18%、2回で0.00046%。

最初の見当を少し良くするより、ニュートン法を1回かけるほうがずっと効く。

シミュレーション

3つの定数を全部の数で試す

開いたときと「数え直す」を押したときは、1以上4未満の単精度の数16,777,216個をその場で順に試し、最初の推定とニュートン法1回のあとの相対誤差の最大を数える。定数を切り替えたときは、同じ計算で数え終えた最大がすぐに出る。 $x$ を4倍すると整数は $2 \cdot 2^{23}$ 増え、半分にして引いた値は $2^{23}$ 減るので、推定はちょうど半分になって誤差は変わらず、ここで数えた最大は正規化数の全体での最大と同じになる。 $x$ のつまみかグラフを横になぞると、丸い拡大鏡の中に、縦の向きだけを拡大した最初の推定（■）とニュートン法1回の値（●）が並ぶ。

考え方

最良の見当が負ける理由

パデュー大学のクリス・ロモンは2003年の論文で、ずらしたときに指数の最後のビットが仮数へ入る場合と入らない場合を分けて、最初の推定の誤差を式に書いた。そのうえで、最大の誤差がいちばん小さくなる定数を求め、0x5f37642fを得た。

ところが、全部の単精度の数で試すと、ニュートン法を1回かけたあとは元の0x5f3759dfのほうが誤差が小さかった。ロモンは近くの定数を1つずつ試して0x5f375a86を見つけ、なぜ最良の見当が1回あとの最良にならないのかを、論文の最後に問いとして残している。

筆者はQuake IIIの0x5f3759dfと、ロモンが式から求めた0x5f37642f、探し当てた0x5f375a86を、正の正規化数2,130,706,432個すべてで数え直した。ニュートン法も、C言語のfloatと同じく、掛け算と引き算を1回ごとに単精度で丸めている。表の「推定」は最初の推定が上と下に外れる最大、「1回」と「2回」はニュートン法をかけたあとの最大の相対誤差である。

定数	推定（上）	推定（下）	1回	2回
0x5f3759df	3.3960%	3.4376%	0.17523%	0.00047%
0x5f37642f	3.4213%	3.4213%	0.17759%	0.00049%
0x5f375a86	3.3976%	3.4365%	0.17513%	0.00047%

ロモンの表とも、ほぼ同じ値になった¹。

違いは、推定の上と下の欄に出ている。

手順5の式で、上に外れた見当（ $δ > 0$ ）には $(3 + δ) /2$ 、つまり1.5より大きい係数が掛かり、同じ大きさだけ下に外れた見当より、1回あとの誤差が大きくなる。0x5f37642fでは上下とも3.4213%にそろっているので、上側が1回あとに0.1776%まで残る。0x5f3759dfでは、下に3.4376%と大きく外れる代わりに、上は3.3960%で止まっている。1回あとの最大は、下側から来る0.1752%で済む。

0x5f375a86は、その釣り合いをもう一段詰めた定数だ。上側の3.3976%と下側の3.4365%を手順5の式に入れると、どちらも0.17512%になる。2016年にはモロズらが、各段階の誤差を式で表して最適な定数を求め直し、ニュートン法1回のあとの相対誤差を最小にする定数として、同じ0x5f375a86と最大誤差 $1.75118 \times 1 0^{- 3}$ を得ている。

意外なのは、Quake IIIの定数がこの最適な値とほとんど変わらないことだ。2つの定数の差は167で、1回あとの最大誤差の差は0.0001%しかない。最初の見当ではなく、ニュートン法をかけたあとの誤差を見ながら選んだ定数なのかもしれない。ただし、これは表から読める推測にとどまる。

歴史

王立協会の議事録とfdlibmの草稿

ニュートン法は、ニュートン=ラフソン法とも呼ばれる。セント・アンドルーズ大学のMacTutorにある伝記によれば、ジョセフ・ラフソンは1690年の著書『方程式の普遍解析』でこの方法を示した。ニュートンも『流率法』に同じ方法を書き、 $x^{3} - 2 x - 5 = 0$ の2と3の間の根を求めている。ただ、この本は1671年に書かれたものの、出版は1736年だった。

同じ伝記は、1690年12月17日の王立協会の議事録を引いている。議事録によれば、ラフソンの方法は1回の操作ごとに「すでに分かっている根の桁数を倍にする」もので、3、4回くり返せば8桁から10桁まで正しくなる。上の表で誤差が3.4%、0.18%、0.00047%と進むのは、正しい桁がおよそ1.5桁、2.8桁、5.3桁と、ほぼ倍ずつ増えていく様子にあたる。

ビットの並びから逆平方根の見当をつける考えは、Quake IIIより前の文書にも残っている。サン・マイクロシステムズの数学ライブラリfdlibmにある平方根のソースe_sqrt.cには、1986年5月にウィリアム・カハンとK. C. Ngが書いた草稿の写しが付いている。草稿によれば、倍精度の数の上位32ビットを整数として扱い、0x5fe80000から半分を引いて、さらに表にある補正の値を引くと、 $1/ x$ をほぼ7.8ビットの精度で「魔法のように（magically）」近似できる。続く反復はy := y*(1.5-0.5*x*y*y)で、Quake IIIの4行目と同じ形だ。草稿はこの反復を、割り算を避けるために使うと説明している。

Quake IIIの定数を、だれが最初に考えたのかははっきりしない。

ロモンは2003年2月の論文の冒頭で、ゲーム開発者向けのサイトにある数学のフォーラムで、このコードに出会ったと記している。作者については、名前の挙がった人はいるが確かめられないとして、知っている人はいないかと問いかけ、論文の終わりでもこの問いは開いたままだとしている。筆者も、定数を最初に書いた人を示す一次資料は見つけられなかった。

公開されたソースでは、同じQ_rsqrtが、PowerPCの機械向けには別の形で書かれている。最初の見当にはfrsqrteというCPUの命令を使い、そのあとに同じニュートン法の1行を置く。

見当をつける部分だけを、機械に合わせて取り替えていたわけだ。

使いどころ

今のCPUとGPUで1/√xを求める

今の計算機には、逆平方根の見当を1つの命令で返す仕組みがある。インテルのIntrinsics Guideによれば、SSEの命令RSQRTSSの相対誤差は $1.5 \times 2^{- 12}$ 、約0.037%より小さい。Quake IIIの見当にニュートン法を1回かけた0.175%と比べると、最大の誤差は5分の1ほどだ。AVX-512の命令VRSQRT14PSでは、 $2^{- 14}$ 、約0.0061%より小さくなる。NVIDIAのGPU向けの命令セットPTXの文書によれば、命令rsqrt.approx.f32の最大の相対誤差は、正の有限の範囲全体で $2^{- 22.9}$ 、約 $1.3 \times 1 0^{- 7}$ で、単精度の丸めの大きさに近い。

速さは機械とコンパイラーで変わる。ロモンは2003年の論文で、当時のVisual C++ .NETでは(float)(1.0/sqrt(x))のおよそ4倍速かったと書いているが、今もそのまま成り立つとは限らない。逆平方根の命令がある機械では、命令で見当をつけ、要る精度に合わせてニュートン法を0回か1回足せばよい。Quake IIIのPowerPC向けのコードと同じ組み立てである。

それでも自分でビットを読み替えるなら、元のコードはそのまま写さないほうがよい。1つは型の大きさで、Q_rsqrtはlongで4バイトを読むが、longの大きさは機械によって違う。筆者が64ビットのLinuxでgcc 15.2を使って調べると、longは8バイトだった。もう1つはC言語の規格で、ある型で置いた値を別の型として読める場合は限られていて、floatをlongとして読むのはそこに入らない。この決まりを破ったときの動作は未定義とされている²。4バイトをmemcpyでuint32_tに写すか、C++20ならstd::bit_castを使えば、この問題は起きない。

JavaScriptなら、同じ4バイトをFloat32ArrayとUint32Arrayで共有すれば、このページのシミュレーションと同じようにビットの並びを読み替えられる。ブラウザーの開発者ツールで試すなら、次のとおり。

const f = new Float32Array(1), u = new Uint32Array(f.buffer);と打つ。
f[0] = 4; u[0] = 0x5f3759df - (u[0] >>> 1);のあと、f[0]を表示すると0.48310753…になる。
let y = f[0]; y = y * (1.5 - 2 * y * y);で、yが0.4991535…になる。 $x$ を変えるときは、2を $x /2$ に替える。
倍精度に広げるなら、Float64ArrayとBigUint64Arrayを使い、ロモンが倍精度向けに挙げた定数0x5fe6ec85e7de30daから半分を引く。筆者が1から4の間の $2^{26}$ 点で試すと、最初の推定の最大誤差は3.4213%、ニュートン法1回のあとは0.1776%だった。

そもそも長さが要らず、2つの距離を比べるだけなら、平方根を取らずに2乗のまま比べればよい。物理エンジンでのこの使い方は、衝突判定とすり抜けにある。

q_math.cで2回目のニュートン法を書いた行は、今も行頭の2本の斜線で無効にされたまま残っている。0.18%で足りると見た判断の跡なのだろう。

参考文献

id Software, Quake III Arena GPL source release, code/game/q_math.c（Q_rsqrt・VectorNormalizeFast）とcode/game/q_shared.h（PowerPC向けのQ_rsqrt）. GitHub id-Software/Quake-III-Arena（GPL第2版以降。最初のコミットの説明に、2005年8月20日にGPLで公開したソースとある）
C. Lomont, "Fast Inverse Square Root", Purdue University（2003年2月）
L. V. Moroz, C. J. Walczyk, A. Hrynchyshyn, V. Holimath, J. L. Cieśliński, "Fast calculation of inverse square root with the use of magic constant - analytical approach", arXiv:1603.04483 (2016)
W. Kahan, K. C. Ng, 1986年5月の草稿の写し（fdlibm e_sqrt.c, Sun Microsystems, 1993。netlibで公開）
J. J. O'Connor, E. F. Robertson, "Joseph Raphson", MacTutor History of Mathematics, University of St Andrews（2017年5月更新）
T. J. Ypma, "Historical Development of the Newton-Raphson Method", SIAM Review 37(4) (1995), 531-551. DOI 10.1137/1037125（書誌はCrossrefで確認）
Intel, Intel Intrinsics Guide（データ3.6.9、2024年7月12日）の_mm_rsqrt_ssと_mm512_rsqrt14_ps
NVIDIA, Parallel Thread Execution ISA Version 9.4, 9.7.3.16 "Floating Point Instructions: rsqrt"
ISO/IEC JTC1/SC22/WG14 N1570（C11の委員会草案、2011年4月12日）、WG21 N4861（C++20の作業草案、2020年4月1日）

注

ロモンの表は、0x5f3759dfで最初の推定3.43756%、1回0.175228%、2回4.66 × 10⁻⁴%。筆者がニュートン法の途中を倍精度で計算し、各回の終わりだけ単精度に丸めて数え直すと、1回0.175229%、2回4.660 × 10⁻⁴%で、3つの定数とも表とほぼ一致した。ロモンの2回目の値がこの記事の表より少し小さいのは、途中の丸め方の違いで説明がつく。 ↩
C11の委員会草案N1570の6.5節7項（オブジェクトの値を読むときに使える型）と4節2項（決まりを破ったときは未定義の動作）。共用体を通した読み替えは、6.5.2.3節の脚注95が「type punning」として説明している。C++20の作業草案N4861では26.5.3節のbit_cast。 ↩