20世紀の遺跡を訪れる

前世紀の遺物、といってもQuake IIIのソースコードの一部。単位ベクトルを計算するのに必要な $1 / \sqrt{x}$ の計算量が馬鹿にならない。なんとかならないかと先人が工夫した結果がこれだ。

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck? 
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

なんのこっちゃ判らないので、説明してくれてる動画をどうぞ。

www.youtube.com

この動画見てなるほどって言える人は情報系の学校行った人なんかだと思われます。ソフォモアまでの数学と計算機科学と英語の単位が必須です。これだけでは何なので、普通の人でも分かった気になれそうな説明をしてみます。暇ならどうぞ。

なにはともあれ、IEEE754

突然ですが、浮動小数点を計算機で扱うときの規格を理解するのが結果的に早道かと。浮動小数点つうのは、

$\displaystyle 1.38064852 \times 10^{-23}$

みたいな表記法で、小数の表現に便利。で、計算機内には 4バイト = 32ビットを使って、この小数が格納されており、そのための規則が IEEE754 つうわけ。詳しくはWikipediaでも見てもらうとして、大事なのは、MSBから、正負の符号 1ビット、指数部 8ビット、残りLSBまで仮数部 23ビットの、合計32ビットということだ。

符号 (Sign)は、0がプラス、1がマイナスを示す。
指数部 (Exponent)は 8ビットなので、[0 .. 255] の整数を表現できるが、規格では 127 を引いて、[-127 .. 128] に対応させる。
仮数部 (Mantissa)は 23ビットで、これを全て小数点以下に当てると 0 ≦ M ＜ 1 を表現するが、規格ではこれに 1 を足したものを仮数としている。つまり、1 ≦ M ＜ 2 となる。

例えば、1.75E-1 を例に取ると、1.75 から1を引いて 0.75、これを二進数にすると1100 0000...。仮数の-1 に127を足して126、二進数で0111 1110 となり、ビットは、0 0111 1110 1100 0000 0000 0000 0000 000 となる。

ひとつめのトリック

なぜ長々と浮動小数点の二進表記まで説明したのか。このトリックの説明にどうしても必要だから。先ずはさっき例に出した 1.75E-1 の内部表現を見てほしい。

0 0111 1110 1100 0000 0000 0000 0000 000

やっぱりなんのことだが分からないよね。それで良いんだ。計算機内部に保存されてるデータそのものに意味はなくて、意味を付けているのは使用している人間の方なんだ。だからこの01の並びをIEEE754に従って浮動小数と解釈するのも、整数と解釈するのもこちらの勝手なのだ。先程の浮動小数点の指数部を $E$ 、仮数部を $M$ として、

$\displaystyle x \equiv ( 1 + \frac {M}{2^{23}}) \times 2 ^{E - 127}$

が表現される小数だ。逆に、上の内部表現を無理やり整数として読み取るとこんな感じになる。

$\displaystyle M + E \times 2^{23}$

そうなんだ勝手なのよ。この２のべき乗は桁を左右にずらす意味だよ。んじゃ、さっきの小数を $x$ として、 $\log x$ を考える。ここで言う $\log$ は底が2なので、例えば、 $\log (8) = 3$ となる。

$\displaystyle \log ( ( 1 + \frac {M}{2^{23}}) \times 2 ^{E - 127} ) = \log ( 1 + \frac {M}{2^{23}}) + E - 127$

$\log ( 1 + \frac {M}{2^{23}})$ は、 $0 \leq \frac {M}{2 ^ {23}} \lt 1$ であることから、 $\log ( 1 + \frac {M}{2^{23}}) \approx \frac {M}{2 ^ {23}} + \delta$ と近似できる*1。 $\delta$ は補正定数、後で説明する。よって、

$\begin{eqnarray} \log ( 1 + \frac {M}{2^{23}}) + E - 127 & \approx & \frac {M}{2 ^ {23}} + \delta + E - 127 \\ & = & \frac {1}{2 ^ {23}} (M + E \times 2^{23}) + \delta -127 \\ \end{eqnarray}$

となる。で、 $M + E \times 2^{23}$ の部分に注目すると、これは、計算機の内部表現を無理やり整数として解釈したものと同じだ。やっやこしいが、こういうことだ。浮動小数の内部表現を無理やり整数として解釈したもの、それを $2^ {23}$ で割って、 $\delta$ 足して、127引くと、 $\log (x)$ が求められる。衝撃だよね、ダマサれた感あるね。引用した説明動画のコメントにsloppy math とblack magic の中間だ、と表現したものがあったけど、そのとおりだね。

マジックナンバー

平方根の逆数じゃなくって $\log$ のはなしばかりでしたが、いよいよ本題に入ります。 $\Gamma \equiv \frac {1}{\sqrt {x}}$ として、 $\log (\Gamma)$ を考えます。

$\begin{eqnarray} \log (\Gamma) & = & \log (\frac {1}{\sqrt {x}}) \\ & = & \log (x ^ {- \frac {1}{2}}) \\ & = & - \frac {1}{2} \log (x) \end{eqnarray}$

$\Gamma$ と $x$ の仮数部、指数部をそれぞれ、 $M_\Gamma$ , $E_\Gamma$ , $M_x$ , $E_x$ として、上の等式を書き換える。 $\log (x) = \frac {1}{2 ^ {23}} (M + E \times 2^{23}) + \delta -127$ であったので、

$\begin{eqnarray} \frac {1}{2 ^ {23}} (M_\Gamma + E_\Gamma \times 2^{23}) + \delta -127 & = & - \frac {1}{2} ( \frac {1}{2 ^ {23}} (M_x + E_x \times 2^{23}) + \delta -127 ) \\ \frac {1}{2 ^ {23}} (M_\Gamma + E_\Gamma \times 2^{23}) + (\delta -127) & = & - \frac {1}{2} (\delta -127) - \frac {1}{2} ( \frac {1}{2 ^ {23}} (M_x + E_x \times 2^{23} ) \\ M_\Gamma + E_\Gamma \times 2^{23} & = & - \frac {3}{2} 2 ^ {23} (\delta -127) - \frac {1}{2} (M_x + E_x \times 2^{23}) \end{eqnarray}$

$\begin{eqnarray} なんらかの整数 & = & マジックナンバー - \frac{1}{2} \times 小数xを無理に整数で解釈したもの \end{eqnarray}$

となります。ここがソースコードの

	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?

の部分です。二分の一倍はビットシフトで書き換えてますね*2。

小数の内部表現を無理に整数で解釈するというのが、この部分。

i  = * ( long * ) &y;                       // evil floating point bit level hacking

その逆、整数の内部表現を無理に小数で解釈するというのが、これ。

	y  = * ( float * ) &i;

ポインタのキャストのデリファレンス？で良いのか？よくわからんが確かにevil だ。現代風に、コンパイラが文句言わないようにunion とかで置き換えて理解するのが良いかも。

このマジックナンバー、ソースでは 0x5f3759df と記されています。仮に $\delta = 0$ とすると、0x5F400000 となりますから、何らかの味付けがされているようです。逆算すると $\delta = 0.045...$ となります。作者は数学的なアプローチではなくトライアンドエラーで決定したとも伝えられています。また、このマジックナンバーについての学術論文もあるようです。

仕上げも念入りに

狐にダマサれたまま $\frac {1}{\sqrt{x}}$ の近似値が出てしまいました。この近似値はなかなかイイ線いってて、こんな感じです。

f:id:rikipoco:20210820181949j:plain

赤線が近似値、黒線が計算値です。全体的に下ぶれしてます。これを補おうというのが $\delta$ だった訳です。しかし、ここで手を抜いてるようでは一流とは言えません。さらに仕上げの一発をぶっ込んできます。ニュートン法です。関数の解を繰り返しの数値演算で求めるものです。ソースのこの部分がそうです。

	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration

$f(y) \equiv \frac {1}{y^2} - x$ とすると、 $f^\prime (y) = - \frac {2}{y^3}$ となるので、

$\begin{eqnarray} y_{n + 1} & = & y_n - \frac {f (y)}{f^\prime (y)} \\ & = & y_n - \frac {1 / y_n^2 - x}{- 2 / y_n^3} \\ & = & y_n + \frac {y_n - xy_n^3}{2} \\ & = & \frac {3}{2} y_n - \frac {xy_n^3}{2} \\ & = & y_n \times (\frac {3}{2} - \frac {x}{2} y_n^2) \end{eqnarray}$

となる。ソースコードでは二度目のニュートン法はコメントアウトされてます。一回で充分なんす。

時は流れて

対数と浮動小数点の内部表現をつかったトリックは、数学者でもプログラマでもなく天才の閃きとしか言えない。さらに、ニュートン法でキッチリ落とし前をつける下りは見事、着地満点。しかし、1995年に発表されたIntelのx86の命令セットのSSEによってこのコードは墓場行きとなりました。そうなんです、本当に前世紀の遺跡だったのです。Quake IIIの戦場が未来の遺跡に見えてきましたか？見えませんか、そうですか。

*1:テイラー展開すれば確認できる

*2:ふだんつかってる十進法で、桁を左に、つまり小数点を右にひと桁ずらすと数は十倍に、その逆なら十分の一になるように、二進法にも同様の演算が有ります。二進法では、桁を左に、つまり小数点を右にひと桁ずらすと数は二倍に、その逆なら二分の一になります。