【ルート高速計算の仕組み】n乗根ビットハックの導出プロセス

計算機で \(n\) 乗根（\(\sqrt{x}\) や \(\sqrt[3]{x}\)）を計算する際、ニュートン法などの反復計算を用いるのが一般的です。しかし、反復計算の速度は「初期値をいかに真値に近づけるか」に大きく依存します。

ここで登場するのが、特定の定数を用いたビット演算で初期値を一気に求める「魔法」のような手法です。実はこの魔法の正体は、「浮動小数点数のビット表現 \(I_x\) は、実は \(\log_2(x)\) の近似値である」という性質を利用したものです。今回は、そのマジックナンバーがどのような理屈で導出されるのか、その数学的背景を解説します。

マジックナンバーの仕組み

STEP 1: 浮動小数点数を数式化する

IEEE 754形式（\(f64\) C言語におけるdoubleと同等の浮動小数点数）の正の数 \(x\) は、符号部を除くと以下のように表せます。
\[x = (1 + f) \cdot 2^{E – B} { …[1]}\]

\(f\): 仮数部（\(0 \le f < 1\)）
\(E\): 指数部の整数値
\(B\): 指数バイアス（\(f64\) の場合は \(1023\)）

これをビット表現（整数として解釈したもの） \(I_x\) にすると、以下のようになります。
\[I_x = L \cdot E + L \cdot f = L(E + f){ …[2]}\]

\(L\): 仮数部のビット数による重み（\(2^{52}\)）

STEP 2: 対数との橋渡し（対数近似）

一番初めに提示した式[1]に対して、\(x\) の底を 2 とする対数を取ります。
\[\log_2(x) = \log_2(1 + f) + (E – B){ …[3]}\]

ここで、\(\log_2(1 + f)\) は \(0 \le f < 1\) の範囲において、直線 \(f\) で近似できます。ここに補正値 \(\sigma\) を加えると精度が上がります。
\[\log_2(1 + f) \approx f + \sigma \quad (\sigma \approx 0.04303){ …[4]}\]

y = log2(1+f )と直線 y = fのグラフ

これを代入すると：
\[\log_2(x) \approx f + \sigma + E – B = (E + f) + \sigma – B\]
STEP 1 の \(I_x = L(E + f)\) を変形して \((E + f) = \frac{I_x}{L}\) を代入すると、「ビット表現と対数の関係式」が完成します。
\[\log_2(x) \approx \frac{I_x}{L} – B + \sigma\]

STEP 3: n乗根の計算

求めたい値を \(y = x^{1/n}\) とします。両辺の対数を取ると：
\[\log_2(y) = \frac{1}{n} \log_2(x)\]
ここに STEP 2 で作った関係式を \(x\) と \(y\) それぞれに当てはめます。
\[\frac{I_y}{L} – B + \sigma \approx \frac{1}{n} \left( \frac{I_x}{L} – B + \sigma \right)\]

STEP 4: \(I_y\) について解く

あとは \(I_y\) について整理するだけです。

両辺に \(L\) を掛ける：\[I_y – L(B – \sigma) \approx \frac{1}{n} (I_x – L(B – \sigma))\]
整理する：\[I_y \approx \frac{1}{n} I_x – \frac{1}{n} L(B – \sigma) + L(B – \sigma)\]
\[I_y \approx \frac{1}{n} I_x + \left( 1 – \frac{1}{n} \right) L(B – \sigma)\]
最終形態：\[I_y \approx \frac{1}{n} I_x + \frac{n – 1}{n} L(B – \sigma)\]

これが、一般式となり、マジックナンバーの導出ができます。

逆根のマジックナンバー

平方根の超高速で安定した計算ができ、その便利さからCPUの命令レベルで組み込まれたアルゴリズム、1990年代の伝説的ゲーム『Quake III Arena』で使われた、「魔法の数字」によるハックの64ビット版も、今回と同じような導出で逆平方根を求めると計算できます。(導出は略しますがほぼ同じです)

\[I_y \approx -\frac{1}{2} I_x + \frac{3}{2} L(B – \sigma)\]

計算例：平方根（\(n=2\)）の場合

例えば \(f64\) で平方根（\(n=2\)）を求める場合、マジックナンバー \(M\) は以下のようになります。
\[M = \frac{2 – 1}{2} \cdot 2^{52} \cdot (1023 – 0.04303)\]
これを計算すると、よく知られた 0x1ff.. で始まる巨大な定数が導き出されます。(補正項の設定によって、0x1ff7a7dfe32a0600; や有名なものとして0x1ff7a7d51d1520f2などがあります)

逆平方根の計算は以下のようになり、以前の記事で紹介した電卓の計算に使用しました。

\(L = 2^{52}\) （\(f64\) の仮数部ビット数）
\(B = 1023\) （\(f64\) の指数バイアス）
\(\sigma \approx 0.0450465\) （最適化された補正値逆平方根用）

\[M = 1.5 \times 2^{52} \times (1023 – 0.0450465)\]
\[M \approx 1.5 \times 4,503,599,627,370,496 \times 1022.9549535\]
\[M \approx 6,910,469,410,797,115,305\]
16進数への変換:
\[M \approx \mathbf{0x5fe6eb50c7b537a9}\]

0x5fe6eb50c7b537a9は平方根の計算機能作成で使った数と一致しました。