Brainfuckでレイトレーサーを書く

C++ のシステムプログラミングコンテストに向けて、久しぶりに CMake を学び直すことにした。その間 Cargo に甘やかされすぎていたせいもある。そんな中、チュートリアルで面白い記述を見つけた。

多くの場合、正しいやり方は汎用プログラミング言語でツールを書いて問題を解決し、そのツールをビルド中にどう呼び出すかを CMake に教えることだ。コード生成、暗号署名ツール、さらにはレイトレーサーまで、CMake Language で書かれた例があるが、これは推奨されるやり方ではない。CMake Language Fundamentals

以前レイトレーサーを書いたことがあり、その後 GPU 向けに書き直したこともある。だからこの一文がやけに引っかかり、考え始めた。レイトレーサーを書くのに、これより適した言語があるだろうか？ 前回書き直したときは、ほとんど第一原理から出発せず、かなり複雑な API に依存したコードを書いていた。そこで、自分が知る中で最も単純な言語——BF を選んだ。単純な言語なら、当然ながら非常に単純なコードベースになるはずだ。実際、BF のコードベースは数行に収まることが多い。しかも Muller が README に残したあのコメントを見て、反例を挙げてみたくなった。コードは mTvare6/rayfuck にある。

はじめに

BF は極めて単純な言語で、8 つの命令と 1 つの「データ構造」しかない。片側に無限に伸びるセルのテープで、各セルには u8 を 1 つ格納できる。文字 > を見たら、テープ上のどこかのセルを指すデータポインタを右へ動かし、I/O is managed through なら左へ動かす。and。. The first stores the input byte where the data pointer points and the latter prints it out.The only primitives other than I/O which allow changing a value are increment and decrement at the data pointer, through +and-.The limitations should be obvious: there are no n > 1 registers as other machines tend to have, no instruction operating on more than one cell, and no instructions for addition or multiplication.The last ingredient required to make BF Turing-complete is its loop, written using [and]. When the token [is met, the runtime checks the cell at the data pointer: if it is zero, execution jumps past the matching], otherwise it enters the loop. At ], it returns to the matching [if the cell is non-zero and exits the loop otherwise.A quick exercise would be writing acat` プログラムを書いてみて、5 文字だけで書けるか試し、この環境がどんなものか掴んでほしい。

下準備

事前に関連資料を読んだ上で、結果や実装の詳細は一切調べず、できるだけ第一原理から書き下すことに決めた。範囲をできるだけ小さくしつつ、明らかにレイトレーサーだと分かるプログラムにするため、RIW の Metal の節がレンダリングするのとまったく同じ画像を出力させることにした。とはいえあの C コードはやはり複雑すぎて、C パーサを書くのは明らかに範囲外だ。誰も保守しない C パーサを書くような仕事は Anthropic に任せておけばいい。各 double（および bool など同種のデータ型）は複数のセルを組み合わせて表現し、ビットの半分を小数部、残り半分を整数部に充てることにした。つまり両者の間に固定の二進小数点を置くイメージだ。後にこれが Q 形式と呼ばれるものだと知った。安価な符号付き Q8.8 なら分解能は 1/256、範囲はおよそ [-128, 128) になる。だが明らかに足りない。地面として使う球体が平らに見えるには r=1000 が必要だからだ。そこで高価な符号付き Q16.16 形式を使った。分解能は 1/2^16、範囲は [-2^15, 2^15) で、これで足りる。コードは SSA に近い形へ変換することにした（そしてこれが LLM に任せる唯一の作業だと決めた）。再帰コードは反復に直し、関数内で定義する変数にはハンガリアン風の接頭辞を付けて、名前でアドレスを引くときの衝突を避ける。同様に、パースとコード生成を分離するのは必須だと判断した。こうすれば複雑さが 2 つのコード領域に分かれ、その間を「DSL」が IR として橋渡しする。この DSL には abs、add、and、call、copy、div、else、end、eq、func、ge、gt、if、int、le、lt、mul、neg、not、or、print2、print3、set、sqrt、sub、text、var、while といった単純な操作が含まれる。次に厄介だったのがいくつかのライブラリ呼び出しだ。使ったのは sqrt、rand、abs である。当初は次のような 2 状態の方式を考えていた。

A = (A-B) % 256
B = (B+1) % 256
or
B = (B+A+p) % 256 # for some prime p

だがこれらの変種は周期がおおむね悪い。結局、より単純な

A = (5*A + 1) % 256

を選んだ。256 個の値すべてを経てから初めて繰り返すことが保証されており、この用途（スーパーサンプリングによるアンチエイリアス）には十分だ。

sqrt明白な候補が 1 つある。Heron の公式だ（ちょうど同じ CMake チュートリアルで復習したばかりだった）。しかし除算が絡むため、結果がひどくなるのは明らかだ。繰り返し減算は生成されるコードが小さく（これは良い。インタプリタが動かす命令が減る）、それでも実行はかなり高くつく。 残る候補は Taylor 級数と「学校で習う方法」で、後者は長除法を使う。

sqrt(1 + x) = 1 + x/2 - x^2/8
y = 2^16*x
sqrt(2^16 + y) / 2^8 = (1 + y/2^17 - y^2/2^35 )
sqrt(y) / 2^8 = (1 + (y - 2^16)/2^17 - (y - 2^16)^2/2^35 )

Desmos にプロットしてみると、符号化値が 20k 未満、つまりおよそ0.305未満で当てはまりが悪い。しかもそこはかなり重要な区間だ。 こうして残るのは長除法だけで、方法は単純である。真の値をxとすると、表現値は次のようになる。

N = x * 2^16

sqrt(x)を表すには次が必要だ。

N' = sqrt(x) * 2^16
isqrt(N) = sqrt(x) * 2^8
isqrt(2^16 * N) = sqrt(x) * 2^16 = N'

ここでisqrtは妥当である。符号化結果が 1 ずれても、復号後の平方根の変化は1/2^16未満、およそ0.00001526だからだ。

最後に、変数マッピング全体と対応する BF アドレスを 1 つの辞書で管理する。

実装

理論部分が固まれば、残りの実装はどれも単純だ。重要なプリミティブは move と copy の 2 つである。

[a, 0]

Move は、初期データポインタのあるセルがゼロになるまで値を減らし続け、そのたびに別のセルを同じ量だけ増やす。

[ # start loop
 - # decrement
 >+ # move right and increment
 ++>+0, b->0, a + ... + a]

ここでaは毎回使い切られ、bはaがゼロになると減り、b個のaが生じる。

4 セルの値では、一方の数の各セルをもう一方の数の各セルと総当たりで組み合わせる。iとjのセルを掛け、その結果を 8 セルの結果のi+jに加える。

[a0, a1, a2, a3] * [b0, b1, b2, b3]

result[i+j] += a_i*b_j

2 つの入力はどちらも表現に 2^16 を持つため、結果をコピーする際に下位 2 ユニットは捨てられる。

N_1 = x_1 * 2^16
N_2 = x_2 * 2^16

( N_1 * N_2 ) / 2^16 = N = (x_1 * x_2) * 2^16

除算は少し込み入っているが、それでも筆算と同じ手順でできる。被除数の最上位ユニットから読み始め、各ステップで古い剰余を 1 ユニット分だけ次のユニットへ繰り下げる。

R = R * 10 + A[next] # school
R = R * 256 + A[next] # here

次に、その剰余から除数を繰り返し引き、1 回引くごとに結果のユニットを 1 増やす。剰余が負になったらそのステップを取り消し、次のユニットへ移る。

while R >= D:
 R -= D
 result += 1

除算はユニットごとに行い、後で統合するため、各ユニットで必要な減算は最大 255 回である。

乗算で右シフトが自己膨張を招くのと同じように、除算は左シフトによって情報を失う。そのため、除算の前に表現をいくらかシフトする必要がある。

N_1 = x_1 * 2^16
N_2 = x_2 * 2^16

(N_1 / N_2) * 2^16 = (x_1 / x_2) * 2^16 # bits already lost
(N_1 * 2^16) / N_2 = (x_1 / x_2) * 2^16

これらの演算では、まず符号を外し、入力の一方だけが負なら結果を負にする。

比較演算は同じ、より小さな操作を共有している。2 つのユニットを同時にデクリメントし、少なくとも一方がゼロになるまで続け、その後は差が現れるか、一時コピーが完全にゼロになるまで続ける。最上位ユニットに 2^7 を加える細かい処理もあり、これがないと負数はバイトとしてそのまま見たときに技術的にはより大きな値になってしまう。

00 ... 7f -> positive half
80 ... ff -> negative half

128 を加えてラップアラウンドした後

80 ... ff -> positive half
00 ... 7f -> negative half

ブール値のチェックでは各ユニットを読み、いずれかが非ゼロなら出力を 1 に設定する。これは C の真偽値の扱いに合わせたものだ。構築される表現では、最下位ビットが 1 なら真、すべてのビットが 0 なら偽となる。and、or、not といったブール演算は、このビット表現に基づいて動作する。

符号反転には -x = ~x + 1 のトリックを使う。各ユニット x を[255-x によって]反転し、最下位ユニットに 1 を加えて桁上げする。abs は最後のユニットの最上位ビットだけを調べ、そのビットが立っていればこの符号反転を実行する。

変数名を関数ごとに限定し、SSA 形式のコードを使うという以前の決定を踏まえると、関数の扱いはきわめて単純である。コード生成器は関数名の下に本体を記録し、call func を見るとそれをインライン展開して出力する。

ループがあれば、if はごく単純だ。

[- body ]

else では、別のフラグが 1 から始まり、最初のループ本体でクリアされる。

while も同様である。

condition
[ 
 body
 move to condition and calculate
]

成果物

このプログラムは最終的に 23MB になり、画像そのもの（約 0.9MB）より大きい。圧縮手段としてはまったくお粗末だ。

ざっと見積もると、1 分間に 100 回のレイ計算、つまり 1 分に 1 ピクセルである。画像は 400x225 なので、このまま最適化なしで計算すると、私のノート PC ではおよそ 62.5 日かかる。だが、自分が見ているのは空の部分だけで、光線が球の間で跳ね返るたびに予想時間ははるかに延びることに気づいた。そこで上の画像は C コードでレンダリングした近似結果である。執筆時点で生成された 1229 ピクセル（全体で 9 万）のうち、異なるものは 10 個だけで、しかもそのほとんどは 1 しか違わない。

まだ最適化の余地はある。地面をもっと粗く近似できるなら、精度を犠牲にしてアクセスするセル数を減らすのが第一候補だ。ランダムベクトルの正規化ステップも省略できるが、それでは散乱分布が変わり、ここで再現しようとしているコードではなくなってしまう。

更新

Reddit のスレッドで、fork/join プリミティブを使えばどう改善できるかと尋ねられた。面白い問題だったので飛び込み、使っている JIT インタプリタの最適化を進めた（以前の作業にも助けられた）。その結果、性能は大幅に向上した。実際にレンダリングすると、精度誤差のせいか、どこかゴッホの絵のようになった。

← すべての記事