1 つ入れて、6 つを測った — Go の SIMD 実験から Mere に入れたもの

9 月に、Go のチームが SIMD の実験についての記事を書いた。Go 1.26 は amd64 向けの SIMD を GOEXPERIMENT=simd の下で出していて、その API は CPU の形に沿ったものだった。Go 1.27 はそれを、 arm64 の NEON と WebAssembly も覆う、CPU に依存しない simd パッケージにした。まだ実験機能で、 どう考えてそうしたかが丁寧に書かれている。

Mere は私が自分のために作っている小さな言語だ。インタプリタと、コンパイルする 4 つのバックエンド —— C、LLVM IR、WebAssembly、RISC-V の機械語 —— を持っていて、SIMD もすでにあった。普通のループを clang がベクトル化できるようにするコンパイラのパスと、128bit の型 3 つ（f64x2・f32x4・u8x16）だ （数字は前の記事にある）。だから Go の記事が投げかけた問いは 「Mere に SIMD は要るか」ではなく、「Go のどの判断を取り入れれば Mere の SIMD が良くなるか」だった。

Go の設計を 1 項目ずつ見て、決める前にそれぞれ測った。入ったのは 1 つだけ。この記事は、その 1 つと、 入らなかった 6 つと、確かめる途中で出てきたものの話だ。数字はすべて arm64（Apple clang 21、node 24）で測った。

先に表を

Go の機能	Mere	決め手
MulAdd（融合積和）	入れた：fma・f64x2_fma	レーンで書いた頂点変換が 140 ms → 100 ms
ループの端数を扱う LoadVPart / StorePart	入れない	使う唯一のプログラムでの上限が 0.16 s → 0.15 s
レーン幅ごとのマスク型と IfElse	今は入れない	select を要るものが無く、余分な命令は型を増やさずに消せる
「必ず動く」ためのエミュレーション	原則にはしない	厳密に再現できるものだけソフトウェアで計算する。今のところ Wasm の fma だけ
型から幅を消す	入れない	Mere のターゲットはどれも 128bit 幅
feature variants（1.28 で予定）	入れない	実行時に CPU を見分ける問題で、事前にビルドする Mere には起きない
総和系（Go 1.27 には無い）	すでにある	reduce_add・any_true・first_true

最初の推しは間違っていた

記事を読んだ時点では、端数の読み込みを 1 番に推すつもりだった。Go の LoadFloat32sPart(x[i:]) は、 ループの最後の短いベクトルを読み込み、足りない分をゼロで埋め、何要素読んだかを返す。だからループの 最後までレーンのまま処理できる。Mere のバイト単位のプログラムは、最後の数バイトのためにループを スカラで二重に書いているはずで、それが消える —— と思っていた。

作る前に確かめたら、その前提はファイルを 1 つ開いただけで崩れた。UTF-8 の検証器は、端数をすでに ユーザーコードでゼロ埋めして 1 ブロックとして処理していた。ゼロは ASCII なので、埋めても答えは変わらない。 端数と 16 バイト未満の行をスカラで処理しているのは、mgrep（Mere で書いた grep）の literal prefilter だけだった。

そこで、それで上限を測った。81 MB のソースコードの行（170 万行、平均 49 バイト）で、スカラの経路を 通るバイトは全体の 14%。

mgrep の変種	先頭バイトがまれ	先頭バイトが多い
今のまま（端数はスカラ）	0.16 s	0.31 s
端数をユーザーコードでゼロ埋め	0.26 s	0.38 s
端数を丸ごと省く —— 答えは間違う、上限	0.15 s	0.27 s

組み込みの端数読み込みで届くのは、よくて最後の行までだ。しかもその行は端数での照合まで省いているので、 実際より良く見えている。1 番に推すつもりだった機能は、それを使える唯一のプログラムにとって 6% ほどの価値だった。

この測定で、mgrep が本当に損をしている場所も分かった。Mere の SIMD の読み込みは bytes を読み、 行は str なので、行ごとに一度コピーしている。普通の長さの行では、そのコピーが取れるはずの効果の 4 割ほどを食っていた。これは Mere 自身の問題で、Go の設計が答えてくれる話ではないので、別の名前を付けて残してある。

マスクとエミュレーションを手短に

Mere では、2 つの u8x16 を比べると、一致したレーンが 0xFF の u8x16 が返る。NEON と WebAssembly の流儀だ。 RISC-V では、比較はマスクレジスタにビットのマスクを作るので、それをバイト列に戻すのに 2 命令かかる。 WebAssembly では、最初に一致したレーンを訊くと、本物のマスク型があった場合より 2 命令多くかかる。 一方で UTF-8 の検証器は、比較の結果をデータとして使っている。1 と AND を取って継続バイトを数えるのだ。 マスクを別の型にすると、そこに変換が要る。そして select を要るプログラムはまだ無い。 よく使う形をバックエンドでパターンとして拾えば、型を増やさずに余分な命令は消せる。 これは select が必要な JSON の走査器か base64 のコーデックが現れるまで待つ。

Go は、ハードウェアに無い SIMD の演算をすべてエミュレーションするので、simd を使う Go のプログラムは どこでも動く。Mere は、下ろせない演算をそのバックエンドで名前を挙げて拒否する。たとえば RISC-V の バックエンドは f64x2 を拒否する。それに Mere には、コンパイルされたすべてのバックエンドを突き合わせる 基準としてインタプリタがある。Go の GODEBUG=simd=0 がテストの道具として与えるものだ。 原則としてのエミュレーションは Mere にはあまり効かない。だが演算ごとの判断としてのエミュレーションは 大事になった。その演算が、Mere が取り入れた唯一の機能だった。

融合積和

fma a b c は a * b + c を丸め 1 回で計算する。書き下した a * b + c は、積のあとと和のあとで 2 回丸める。 2 つの答えは最後のビットで違うことがあり、積と c がほとんど打ち消し合うときにいちばん違う。 fma 0.1 10.0 (-1.0) は 5.55e-17 で、0.1 * 10.0 - 1.0 はちょうど 0.0 になる。

Mere はこれを意図して持たずにいた。C バックエンドは #pragma STDC FP_CONTRACT OFF を出力している。 arm64 の clang は -O2 で、a * b + c を勝手に fma 命令にまとめるからだ。それが起きたとき、内積が インタプリタと最適化したバイナリで違うビットを出した。Mere には「最適化で答えを変えてはいけない」という 規則がある。まとめるのを止めればその規則は守れるが、積和の連鎖を速くする唯一の命令も手放すことになった。

Go の記事は、別の問いを立てるきっかけになった。融合した答えを、プログラムが名前で頼めるようにしたらどうか。 それなら規則は保たれる。a * b + c はどこでも 2 回丸め、fma はどこでも 1 回丸める。

1 行も書く前に、どれだけの価値があるかを測った。頂点変換（20,000 頂点を 4x4 行列に通すのを 5,000 回）について、 Mere が出力した C の pragma を ON にしたのだ。

カーネル	pragma OFF（出荷時）	pragma ON	融合命令の数
スカラ	0.20 s	0.15 s	0 → 348
f64x2 のレーン	0.14 s	0.10 s	0 → 72

どちらも 1.4 倍。答えは変わる —— 変わるはずだ —— が、融合したスカラ版と融合したレーン版は、 互いに同じビットを出した。fma は決定的なので、融合した答えは丸めの偶然ではなく、きちんと定義された答えになる。 組み込みを書いたあと、f64x2_fma で書き直したレーン版は、C バックエンドで 140 ms に対して 100 ms で走る。 clang はレーンの fma を 2 つずつ fmla.2d 1 つにまとめる。

fma 命令の無いバックエンド

5 つのうち 4 つは素直だった。インタプリタは OCaml の Float.fma、C は fma(3)、LLVM は llvm.fma を呼ぶ。 RISC-V のバックエンドは倍精度の演算器を持たないので、名前で拒否する。本当の判断は WebAssembly だった。 fma 命令が無く、JavaScript のホストにも Math.fma が無い。

選択肢は 2 つあった。Mere がいつもするようにそこで拒否するか、Go がするようにソフトウェアで計算するか。 ソフトウェアを選んだ理由は 1 つだ。fma は IEEE-754 が正しい丸めを要求する演算の 1 つで、 +・-・*・/・sqrt と同じく、正しい答えがただ 1 つあり、それは再現できる。 拒否しても、誰も問題から逃れられない。fma を使いたいプログラムは WebAssembly 版だけ a * b + c と書き、 そこで違う答えを得ることになる。それこそ pragma が防ごうとしていたものだ。

最初のソフト fma は、Mere で 15 行ほどだった。Boldo と Melquiond による技法を使っている。積を 2 つの double として 正確に求め、足し、途中で一度だけ「奇数への丸め」をする。4 つのバックエンドすべてで、3 万件でハードウェアと一致した。 けれど、それが成り立つのは途中で何もオーバーフローもアンダーフローもしない間だけで、確かめたのは 2^-100 から 2^100 くらいの範囲だけだった。しかも WebAssembly では呼ぶたびに確保をして、200 万回ほどで メモリが尽きた。代わりに出荷したのは、WebAssembly のランタイムの中に整数で書いたもので、musl の fma.c と同じ構成だ。

1. 3 つの仮数を 54bit に正規化する。
2. 108bit の積を 32bit ずつに割って作る。WebAssembly には、64bit 同士を掛けて 128bit を返す命令が無いからだ。
3. c を積に揃える。はみ出したビットは sticky ビットにまとめる。
4. 足すか引く。
5. 丸めはちょうど 1 回。上位 63bit を double に変換するところが丸めになる。結果が非正規化数になるときは、 先に非正規化数の最後の桁で丸めておき、最後の 2 のべき乗倍でもう一度丸めないようにする。

確かめる方が、書くより長くかかった。

• まず C で、ハードウェアと 2 億 5,000 万件。特殊値、非正規化数、指数の上下端、打ち消し合い、タイから作った入力で、食い違いは 0。
• 分岐を 1 つずつわざと壊した。10 通りのうち 8 通りで食い違いが出た。残りの 2 通りは、どうやっても答えを変えられなかった。 理由も分かった。53bit 同士の積の真ん中に 64bit 以上のゼロが並び、その両側に 1 が立っている必要があるが、そういう積は存在しない。 なのでこの 2 つは「確かめた」とは言わず、「観測できない」理由を文章で残した。
• 次に WAT そのもの。出力されたモジュールから切り出し、node の上でハードウェアの答えとさらに 5,000 万件を突き合わせた。 食い違いは 0。node での突き合わせ自体も、WAT を壊して落ちることを見て確かめた。

このうちの 1 つには、危うくだまされるところだった。はみ出した c の sticky ビットを落とす壊し方は、 ランダムな入力 2,000 万件では通ってしまった。ランダムな入力は、そのビットだけが丸めを決める唯一の場合であるタイに、 ほぼ当たらない。わざと作ったタイの系列を足すと —— ちょうど半分で丸まる奇数の 54bit の正確な積と、 その遥か下にあって向きを決める小さな c —— 同じ壊し方が 2,000 万件中 10,272 件で落ちた。 タイは常設のゲートに入れてある。コンパイラの WebAssembly のランタイムからそのビットを外すと、ゲートは系列の名前を挙げて落ちる。 Wasm differs from fma(3) in: ties。

WebAssembly でのコストは 1 回あたり約 12 ns。1 回の繰り返しで fma を 8 段入れ子にすると 130 ms で、 融合しない式は 30 ms だった。繰り返しのたびに float を確保するループの中では、この差は雑音に埋もれる。

ゲートが見つけた、fma ではないもの

常設のゲートは、6 系列の入力をインタプリタ、2 つの最適化レベルの C、LLVM、WebAssembly に通し、 それぞれを C ライブラリの fma(3) とビット単位で比べる。このゲートは、初めて緑になる前に、 fma とは無関係な LLVM バックエンドのバグを 2 つ見つけていた。

x != x が NaN を「NaN ではない」と言った。IEEE-754 は NaN のペイロードを決めていないので、ゲートは NaN の結果を 1 つの値にまとめて比べる。NaN かどうかは r != r で見る。LLVM の足は、ハードウェアと 1 系列だけ食い違った。 NaN が出る唯一の系列だ。バックエンドは float の != を順序付きの比較にしていて、片側が NaN だと常に false になっていた。 インタプリタと C と WebAssembly は、どれも true を返す。float の端のケースを見る parity のテストは、NaN について ==・<・> を訊いていたが、!= は一度も訊いていなかった。

組み込みの名前を局所的に束縛すると、それがプログラム全体に漏れた。fma の parity のテストは、fma という名前の 局所関数がそのスコープの中だけのものであることを確かめる。LLVM では、それがファイル中のすべての fma を乗っ取った。 トップレベルのコードが、最後に出力した関数の「外に持ち上げた内側の関数」の表を残したまま出力されていたのだ。 同じバグで、どこかの関数が自分の atan2 を定義すると、atan2 1.0 1.0 が 2.0 になった。 最初に書いた回帰テストは、バグが残ったまま通った。漏れるのは最後の関数の表だけだったからだ。 今のテストは、上書きする関数を最後に置いている。

!= のバグの隣を見ていて、3 つ目も出てきた。f_min と f_max に NaN を渡すと、インタプリタと WebAssembly は NaN を返し、 C と LLVM は引数の順番によってもう一方の値を返した。これは既知の失敗の一覧に入れて、少しあとで直した。 直し方は、OCaml 自身の Float.min を手順どおりにすべてのバックエンドに写すことだ。最初の版は NaN を作るのに a + b を使っていて、ビットを比べるゲートが捕まえた。signaling NaN を渡すと OCaml はそのまま返すが、足し算は quiet NaN にする。

始める前の自分に言うこと

• 機能を作る前に、上限を測る。端数の読み込みは、5 つのバックエンドに組み込みを 1 つ足して 6% だった。端数を丸ごと省いた変種 —— 答えが間違うのも承知で —— が、半日で決着をつけた。
• 規則は「fma を使わない」ではなく「最適化で答えを変えない」だった。pragma がなぜあるのかを読んだことで、 1.4 倍が「禁止」から「使える」に変わった。
• 厳密に再現できるものはエミュレーションし、残りは拒否する。これは全体の方針ではなく演算ごとの判断で、 最初にそれを通ったのが fma だった。
• ランダムなテストはタイに届かない。わざと作る。
• ゲートの収穫は、しばしば対象の隣にある。fma は初日から正しかった。赤くなったのは、ゲートが結果をまとめるのに使った NaN の判定と、 parity のテストがスコープを確かめるのに使った局所束縛だった。

fma には、テストと上の測定のほかに、まだ実際の利用者がいない。いちばん自然なのは小さな 3D エンジンの線形代数で、 そこはすでに、正しく丸める演算についてすべてのバックエンドが同じビットを出すことを確かめている。 判断と測定をすべて書いた詳しい記録は、Mere のリポジトリの docs/simd-go-comparison.md にある。