Quake IIIの高速逆平方根(0x5f3759df)をApple M4とZen 3で実機ベンチ、std::sqrtより速くなる条件がプラットフォームで違った
目次
Xで Q_rsqrt の話が回ってきて、Geminiと一通り議論したあとに「現代のApple M4で動かしたら何が起きるんだろう」と気になって、ベンチを書いて回してみた。
予想は「ハードウェアに frsqrte がある時代に1999年のビットハックが速いわけがない」だった。回してみたら逆で、1.0f / sqrtf(x) より1.55倍速かった。ただ、アセンブリと数字を追っていくと、速い理由はもうビットハック由来じゃなくて、ループが自動でNEONベクトル化されるかどうかという別の話になっていた。自前で vrsqrte を関数化して呼ぶパスは、逆に5倍以上遅くなった。
ついでなので、ニュートン法0/1/2回の精度差と、Chris Lomontの最適化定数 0x5f375a86 も同じベンチに乗せた。さらに別のWindowsノート（Ryzen 7 5800HS / WSL2）でも同じコードを回したら、x86_64ではApple Siliconとは別のところで挙動が割れたので、それも合わせて追記した。最後にアルゴリズムの仕組みをIEEE 754のlog近似まで掘り直して終わる。
ベンチの作り
ベンチコードは記事末尾にフル掲載するが、骨は素直で、
static inline float q_rsqrt_with(uint32_t magic, float number, int newton_iters) {
float x2 = number * 0.5f;
float y = number;
uint32_t i;
memcpy(&i, &y, sizeof(i));
i = magic - (i >> 1);
memcpy(&y, &i, sizeof(y));
for (int n = 0; n < newton_iters; n++) {
y = y * (1.5f - (x2 * y * y));
}
return y;
}元コードはポインタキャストで型punningしていたが、現代のCではUBなので memcpy に書き換えてある。ループ最適化を阻害しないように static inline 化、memcpy は最近のclangならコンパイル時にビット転送に潰れる。
これに対して、
0x5f3759df（オリジナル定数）でニュートン0回 / 1回 / 2回0x5f375a86（Chris Lomont定数）でニュートン0回 / 1回1.0f / sqrtf(x)- 自前ラッパー: NEON
vrsqrte_f32または SSE_mm_rsqrt_ss+ ニュートン1回 / 2回
を、N = 2,000,000 のlog-uniform分布 ( 〜 ) なfloat配列に対して回す。スループットは反復の最小値 / N。精度は 1.0 / sqrt((double)x) を真値とした相対誤差を集計する。
Apple M4 (Mac mini) で回す
環境
| 項目 | 値 |
|---|---|
| 機種 | Mac mini (Mac16,10) |
| CPU | Apple M4 10コア (4 Performance + 6 Efficiency) |
| メモリ | 16 GB unified memory |
| OS | macOS 26.3 (build 25D125) |
| コンパイラ | Apple clang 17.0.0 |
| ビルド | clang -O2 -std=c11 bench.c -o bench -lm |
| 計測 | bench内ループ best-of-5 |
結果
| 実装 | ns/op (M4) | 最大相対誤差 | 平均相対誤差 |
|---|---|---|---|
| Q_rsqrt (0x5f3759df, Newton 0回) | 0.102 | 3.438e-02 | 2.294e-02 |
| Q_rsqrt (0x5f3759df, Newton 1回) | 0.134 | 1.752e-03 | 9.326e-04 |
| Q_rsqrt (0x5f3759df, Newton 2回) | 0.185 | 4.712e-06 | 1.812e-06 |
| Q_rsqrt (Lomont 0x5f375a86, Newton 0回) | 0.084 | 3.437e-02 | 2.295e-02 |
| Q_rsqrt (Lomont 0x5f375a86, Newton 1回) | 0.113 | 1.751e-03 | 9.332e-04 |
1.0f / sqrtf(x) | 0.208 | 8.940e-08 | 2.956e-08 |
NEON vrsqrte + Newton 1回（関数経由） | 0.755 | 1.614e-05 | 2.092e-06 |
NEON vrsqrte + Newton 2回（関数経由） | 0.925 | 1.304e-07 | 3.567e-08 |
-ffast-math を付けても本質的に同じ結果だった。M4の 1.0f / sqrtf(x) はデフォルトで既に十分最適化されていて、fast-math で書き換わる余地がほぼない。
数字を見て一番驚いたこと
Q_rsqrt（Newton 1回）が 0.134 ns/op で、1.0f / sqrtf(x) が 0.208 ns/op。M4のP-coreは概ね4GHz前後で動くので、0.208 ns/op はおおよそ0.8サイクル/要素、0.134 ns/op は0.5サイクル/要素。両方とも「1要素1サイクル未満」の領域に入っていて、これはもうスカラ実行ではなくSIMDで4要素同時並列に走らせている速度。
つまりループ全体が勝手にNEONに変換されている、ということ。手書き Q_rsqrt の中身は、memcpy 2回と整数算術、それにニュートン1回ぶんの3 FMAしかない。clangのループベクタライザはこれを int32x4_t / float32x4_t の幅で4要素同時に展開してくれる。
1.0f / sqrtf(x) の方もARMでは frsqrte + Newton refinement の組み合わせにベクトル化されるパスが存在しているが、sqrtf の標準実装の精度保証（max ulp）が Q_rsqrt より厳しい都合で、refinement のステップ数が多くなる。結果として「Q_rsqrtのほうが命令数が少なくて速い」という順位がそのまま出る。
ビットハック自体が速いのではなく、ビットハックが浅くて単純なせいでベクトル化されやすかった、というのが見ているものの正体。同じ理屈で、ニュートン2回の Q_rsqrt は0.185 ns/op まで落ちて、1/sqrtf との差は0.023 nsまで縮む。精度を上げて命令数を増やすほど、Q_rsqrt の優位は消えていく。
自前NEONラッパーが逆に遅い
予想外だったのがこれ。vrsqrte_f32 を直接呼ぶ自前関数を書くとこうなる。
static inline float neon_rsqrt_1(float x) {
float32x2_t v = vdup_n_f32(x);
float32x2_t e = vrsqrte_f32(v);
e = vmul_f32(vrsqrts_f32(vmul_f32(v, e), e), e);
return vget_lane_f32(e, 0);
}これが 0.755 ns/op、ニュートン2回版なら 0.925 ns/op。Q_rsqrt（Newton 1回）の5倍以上遅い。
static inline を付けているので関数呼び出しオーバーヘッド自体は消えるはずだが、ループ全体の最適化を見てみると、
vdup_n_f32(x)でスカラ→NEONレジスタへの転送が要素ごとに発生vget_lane_f32(e, 0)でNEON→スカラへ戻すのも要素ごとに発生- このスカラ↔NEONの行き来でループ単位のベクタライザが手を引いて、ループそのものはスカラのまま展開される
せっかくのハードウェア命令 vrsqrte を1要素ずつ呼ぶスカラループになってしまう。Q_rsqrt のループが int32x4_t 単位で4並列に走るのと正反対の挙動になっている。
NEONを使いたいときの方針はわかりやすくて、スカラ関数の中にNEONイントリンシックを書かず、ループの外側から float32x4_t 単位で書く。それが面倒なら、素直な式（1.0f / sqrtf(x)）を書いてコンパイラに任せる方が結果的に速い。
ニュートン法は何回入れるべきか
次に気になったのが、Quake IIIが「ニュートン1回」を選んだ判断は本当に正しかったのか、という点。これは数字でほぼ片が付く。
| ニュートン回数 | 最大相対誤差 | 平均相対誤差 | 速度 (ns/op, M4) |
|---|---|---|---|
| 0 | 3.44e-02 | 2.29e-02 | 0.102 |
| 1 | 1.75e-03 | 9.33e-04 | 0.134 |
| 2 | 4.71e-06 | 1.81e-06 | 0.185 |
ニュートン法は二次収束なので、1回ごとに有効桁数がほぼ倍になっている。0回時の誤差 ~3.4% から、1回で ~0.17%、2回で ~5e-6。3Dゲームのライティング・法線計算で0.17%誤差は実用上問題ない（人間の目では判別できない）が、3.4%だと面のシェーディングに目で見える歪みが出る。
Quake IIIが1回を選んだのは、当時のCPU予算的にも誤差曲線的にも妥当な判断だった。コメントの // 2回目は不要（と判断された） は単に「速くしたいから」ではなく、「1回入れれば精度が十分実用域に入る」という二次収束の知識に基づく判断だったと読める。
Lomont定数との差
Chris Lomontが2003年の論文 Fast Inverse Square Root で、ニュートン1回後の最大相対誤差を最小化する定数を解析的に再導出して、0x5f3759df よりわずかに良い 0x5f375a86 を提案した。
M4で実測した結果。
| 定数 | Newton 0回 max_rel | Newton 1回 max_rel |
|---|---|---|
| 0x5f3759df（オリジナル） | 3.438e-02 | 1.752e-03 |
| 0x5f375a86（Lomont） | 3.437e-02 | 1.751e-03 |
ほぼ区別がつかない。ニュートン1回のフィードバックが強くて、初期値の品質差はそこで吸収されている。差が見えるのは厳密にはニュートン0回時だけで、それでも有効桁3桁めまで一致している。
「神の数字」と呼ばれるのはオリジナルの 0x5f3759df のほうだが、ニュートン1回を入れる前提なら定数選びにそこまでこだわる必要はなかった、という話に近い。
x86_64 (Zen 3 / WSL2) でも回してみた
Apple Silicon側で挙動の感触がつかめたところで、別アーキでも同じコードを走らせたくなった。手元のWindowsノートのWSL2上で同じ bench.c をビルドして回す。NEON部分は #if defined(__aarch64__) で自動的に外れて、代わりに _mm_rsqrt_ss を使うx86の小細工を足した。
環境
| 項目 | 値 |
|---|---|
| CPU | AMD Ryzen 7 5800HS (Zen 3, 8C/16T, base 3.2GHz / boost 〜4.4GHz) |
| SIMD | AVX2, FMA, F16C, AES-NI（AVX-512非搭載） |
| ホストOS | Windows 11 Home (26200) |
| 実行環境 | WSL2 / Ubuntu 24.04.3 LTS (kernel 6.6.87.2-microsoft-standard-WSL2) |
| メモリ | WSL2に約7.6GB割当 |
| コンパイラ | clang 18.1.3 (Ubuntu) |
| ビルド | clang -O2 -std=gnu11 -mavx2 -mfma bench.c -o bench -lm |
| 計測 | taskset -c 0 ./bench、bench内ループ best-of-30 |
best-of-30 を3回ぶん取って、その最小値を載せる。WSL2はスケジューラのジッタがそれなりにあって、1回ぶんでは ±50% くらい平気で揺れた。CPUを1コアにピン留めして、十分多めに繰り返してようやく落ち着く程度。M4側の best-of-5 と直接比較する場合はこのぶんだけ条件差があることに注意する。
結果
| 実装 | ns/op (Zen 3 / WSL) | 最大相対誤差 | 平均相対誤差 |
|---|---|---|---|
| Q_rsqrt (0x5f3759df, Newton 0回) | 0.172 | 3.438e-02 | 2.295e-02 |
| Q_rsqrt (0x5f3759df, Newton 1回) | 0.173 | 1.752e-03 | 9.334e-04 |
| Q_rsqrt (0x5f3759df, Newton 2回) | 0.197 | 4.717e-06 | 1.815e-06 |
| Q_rsqrt (Lomont 0x5f375a86, Newton 0回) | 0.174 | 3.437e-02 | 2.295e-02 |
| Q_rsqrt (Lomont 0x5f375a86, Newton 1回) | 0.172 | 1.751e-03 | 9.340e-04 |
1.0f / sqrtf(x) (-O2 のみ) | 2.155 | 8.941e-08 | 2.957e-08 |
1.0f / sqrtf(x) (-O2 -ffast-math) | 0.169 | 1.957e-07 | 3.301e-08 |
SSE rsqrtss（Newton 0回, スカラ） | 0.291 | 2.586e-04 | 6.415e-05 |
SSE rsqrtss + Newton 1回（スカラ） | 0.663 | 2.040e-07 | 3.371e-08 |
Apple Siliconと挙動が割れたところ
M4で「ビットハックが速い謎」が「自動ベクトル化されやすかっただけ」だったのと同じ構図が、Zen 3だとさらに極端に出る。
Q_rsqrtが0.17 ns/op前後で揃うのは、ループがAVX2幅（floatで8要素並列）にきれいに展開されているから。Ryzen 7 5800HSのブースト時クロックを4.4GHz、AVX2を8要素並列とすると、理論ピークは1要素あたり 1/4.4GHz/8 ≒ 0.028 ns。0.17 ns/opは6サイクル/8要素 ≒ 0.75 サイクル/要素相当で、まあ妥当な数字。
問題は 1.0f / sqrtf(x) のほう。-O2 だけだと 2.155 ns/op で、Q_rsqrt の12倍以上遅い。M4の 1/sqrtf が 0.208 ns/op だった事実とは別世界。
理由は単純で、x86のclang 18は -O2 のままだと sqrtf を vrsqrtps（精度の低い近似命令）に置き換えてくれない。sqrtf はIEEE 754の正確丸めを保証する必要があって、vrsqrtps は約12bit精度しかないので、安全側に倒すと vsqrtss のスカラ呼び出しのままになる。実際にアセンブリを見ても、ループは1要素ずつ vsqrtss + vdivss を回している。
-ffast-math を足すと、1.0f / sqrtf(x) は 0.169 ns/op まで一気に縮んで、Q_rsqrt と同じ土俵に並ぶ。clangがこの式を vrsqrtps + ニュートン1回相当のシーケンスに置き換えていて、誤差も Q_rsqrt より2桁良い（max 1.96e-7）ところに収まる。
Apple Clang がM4で -O2 のままでも 1/sqrtf を frsqrte ベースに置き換えていたのと、x86 Clang が -O2 では置き換えない（-ffast-math を要求する）のがそのまま差になっている。同じclangでもターゲットによって浮動小数点に関する安全側のデフォルトが違うので、「M4で測ったときの結論をx86にそのまま持ってくると外す」というのが今回いちばん予想外だった部分。
SSE rsqrtss を直接呼ぶと逆に遅い
M4で vrsqrte_f32 をスカラ関数経由で呼ぶと遅くなる現象（ベクタライザが手を引く）が、そのままx86でも再現する。_mm_rsqrt_ss を関数化したものは 0.291 ns/op で、ニュートン1回足すと 0.663 ns/op まで落ちる。Q_rsqrt より3〜4倍遅い。
rsqrtss は 128bit XMM の下位1レーンだけを処理する命令なので、関数のシグネチャが float -> float のまま _mm_set_ss / _mm_cvtss_f32 で出し入れすると、コンパイラはループ全体を rsqrtps 幅で展開できなくなる。ARMで vget_lane_f32 が壁になっていたのと同じパターン。
x86でもSIMDイントリンシックを float -> float 関数に閉じ込めるとループ単位の最適化が消える。__m256 単位でループの外側から書くか、-ffast-math を許せる文脈なら 1.0f / sqrtf(x) で済ませるのが現実的。
2026年の Q_rsqrt が速い条件
x86_64 / Zen 3で測った結果をM4側の結論と合わせると、こうなる。
- M4 (Apple clang 17,
-O2):Q_rsqrt0.134 vs1/sqrtf0.208。Q_rsqrtが1.55倍速いが、両方とも自動ベクトル化済み - Zen 3 (clang 18,
-O2):Q_rsqrt0.173 vs1/sqrtf2.155。Q_rsqrtが12倍速い。1/sqrtfはvsqrtssでスカラのまま回されている - Zen 3 (clang 18,
-O2 -ffast-math):Q_rsqrt0.173 vs1/sqrtf0.169。ほぼ同等。誤差は1/sqrtfのほうが2桁良い
「2026年でも Q_rsqrt が速い」という見出しは、x86_64の -O2 -fno-fast-math という条件下なら確かに正しい。ただしその速さは「Q_rsqrt がアルゴリズム的に勝っている」のではなく、「Clangが安全側に倒した 1/sqrtf のほうが片手を縛られている」だけ、という読みになる。-ffast-math を許せる文脈では並ぶし、誤差はむしろ 1/sqrtf のほうが小さい。
rsqrtss / vrsqrte をスカラ関数に閉じ込めるとどちらの環境でも遅くなる、という点はプラットフォーム共通の落とし穴。
アルゴリズムを掘り直す
ベンチで挙動が腑に落ちたところで、改めて「なんでこれで動くのか」を数式で追う。ここからが本来 Q_rsqrt 解説の本体。
仕組みは大きく分けて、
- floatのビット表現を整数として読むと、近似的に になる
- は対数の世界で
- 対数空間での演算は、整数ビット表現での「右シフト1bit」+「定数からの引き算」で再現できる
- これでまずまずの近似値が得られるので、最後にニュートン法1回で詰める
IEEE 754 のビット表現が log₂ になる
IEEE 754単精度（float）の32ビット内訳は、
| sign 1bit | exponent 8bit | mantissa 23bit |正の数 を （、 は実際の指数）と書くと、フィールドに格納される値は、
- exponent フィールドの値
- mantissa フィールドの値
float のビット列をそのまま符号なし整数 として読むと、
ここで なので、 を の範囲で1次近似する。最小二乗フィット的に、
を使うと、
floatのビット列を整数として読むと、 にスケールと定数を足したものになる。Q_rsqrt のビットハックはここから出発する。
対数空間で 1/√x を計算する
ほしいのは 。両辺の を取ると、
これを、float のビット列を読んだ整数 の関係に持ち込みたい。さっきの近似 を と の両方に適用すると、
右辺の第2項を整数定数 として括り出す。
すると、
これがコードの1行、
i = 0x5f3759df - (i >> 1);ちょうどそのもの。(i >> 1) が 、0x5f3759df - ... が に対応していて、対数空間での割り算と定数オフセットがこの1行で同時に済む。
ニュートン法で詰める
ビットハックだけだと初期推定値は数%の誤差を持つので、最後に1回だけニュートン法を入れる。 は 、つまり の解として書ける。微分は、
ニュートン反復は、
コード末尾の y = y * (1.5f - x2 * y * y); がこれそのもの。x2 = x / 2 も冒頭で計算されている。
ニュートン法はこの形だと二次収束で、有効桁数がほぼ毎回倍になる。さっきのベンチ表の「0回 ~3% → 1回 ~0.17% → 2回 ~5e-6」がきれいに2倍になっている。
仕組みの全体像
ここまでの流れを図にしておく。
flowchart TD
A["float x の IEEE 754 ビット表現"] --> B["uint32_t として再解釈"]
B --> C["log2(x) のスケール済み近似"]
C --> D["右1bitシフト = log2(x)/2"]
D --> E["0x5f3759df から引く = -log2(x)/2 + 定数"]
E --> F["float として再解釈 ≒ 1/√x"]
F --> G["ニュートン法を1回適用"]
G --> H["最終出力 (最大相対誤差 0.175%)"]「対数 → 引き算 → 指数」を、floatの整数表現の読み替えだけで成立させているのが核。log2 や exp2 を陽に呼んでいないし、対数表もテーブルもいらない。ビットレイアウト自体が対数のスケールを内包しているのを利用している。
それで、現代では Q_rsqrt を書くべきか
書くな、というのが今回測った結果での結論。
Q_rsqrt が速かったのは事実だが、M4とZen 3の両方で、速さの本質は「ループ全体が自動ベクトル化されやすい単純な算術」だった。ビットハックそのものが現代CPUで有利ということではない。同じ理由で素直な式もコンパイラがベクトル化してくれるので、1.0f / sqrtf(x) を -O2 （x86なら -O2 -ffast-math）で書けば、精度・移植性・可読性のどれを取っても上回る。自前で vrsqrte / rsqrtss を float -> float 関数に閉じ込めると逆に遅くなる落とし穴がある点も、両プラットフォーム共通。
順序としては、まず 1.0f / sqrtf(x) を書いてベンチを取る。足りないと分かったらループの外側から float32x4_t / __m256 単位で手書きSIMDに降りる。それでも足りなければ、ループ全体のメモリレイアウトやデータ並列性のほうを見直す。Q_rsqrt はその選択肢のどこにも入らない。1999年の文脈では正解だったが、2026年では「読むコード」枠。
余談: // what the fuck? コメントの出どころ
Q_rsqrt の有名な // what the fuck? コメントは、Quake IIIのソースコードが2005年8月にGPLで公開されたときに表に出た。当のid Software内でもこの定数の由来は完全には伝わっていなかったらしく、John Carmackも後年「自分が書いたコードではない」とコメントしている。
その後の調査で、
- Greg Walsh（Ardent Computer）が80年代後半に類似のトリックを実装していた
- Cleve Moler（MATLAB の作者）が広めたバージョンがベース
- 似たアイデアはKahan/NgのSGI系コードや、IBMのbit-twiddling集に既出
など複数の起源が絡んでいて、1人の発明品ではなく、80–90年代の浮動小数点ハック文化の集大成がたまたまQuake IIIのソース公開で一気に表に出た、というのが現在の理解。
ベンチコード全文
bench.c として動かしたコード（ARM側）。x86側は _mm_rsqrt_ss を使うブロックを #if defined(__x86_64__) で追加した版を使った（記事の流れ上、ARM版だけ載せておく）。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <time.h>
#if defined(__aarch64__)
#include <arm_neon.h>
#endif
static const uint32_t MAGIC_Q3 = 0x5f3759dfu;
static const uint32_t MAGIC_LOMONT = 0x5f375a86u;
static inline float q_rsqrt_with(uint32_t magic, float number, int newton_iters) {
float x2 = number * 0.5f;
float y = number;
uint32_t i;
memcpy(&i, &y, sizeof(i));
i = magic - (i >> 1);
memcpy(&y, &i, sizeof(y));
for (int n = 0; n < newton_iters; n++) {
y = y * (1.5f - (x2 * y * y));
}
return y;
}
static inline float q_rsqrt_quake_n0(float x) { return q_rsqrt_with(MAGIC_Q3, x, 0); }
static inline float q_rsqrt_quake(float x) { return q_rsqrt_with(MAGIC_Q3, x, 1); }
static inline float q_rsqrt_quake_n2(float x) { return q_rsqrt_with(MAGIC_Q3, x, 2); }
static inline float q_rsqrt_lomont_n0(float x){ return q_rsqrt_with(MAGIC_LOMONT, x, 0); }
static inline float q_rsqrt_lomont(float x) { return q_rsqrt_with(MAGIC_LOMONT, x, 1); }
static inline float std_rsqrt(float x) { return 1.0f / sqrtf(x); }
#if defined(__aarch64__)
static inline float neon_rsqrt_1(float x) {
float32x2_t v = vdup_n_f32(x);
float32x2_t e = vrsqrte_f32(v);
e = vmul_f32(vrsqrts_f32(vmul_f32(v, e), e), e);
return vget_lane_f32(e, 0);
}
static inline float neon_rsqrt_2(float x) {
float32x2_t v = vdup_n_f32(x);
float32x2_t e = vrsqrte_f32(v);
e = vmul_f32(vrsqrts_f32(vmul_f32(v, e), e), e);
e = vmul_f32(vrsqrts_f32(vmul_f32(v, e), e), e);
return vget_lane_f32(e, 0);
}
#endif
static uint64_t now_ns(void) {
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC, &ts);
return (uint64_t)ts.tv_sec * 1000000000ull + (uint64_t)ts.tv_nsec;
}
#define N 2000000
static float input[N];
static float output[N];
typedef float (*rsqrt_fn)(float);
static void bench(const char* name, rsqrt_fn fn, int repeats) {
for (int i = 0; i < N; i++) output[i] = fn(input[i]);
uint64_t best = UINT64_MAX;
for (int r = 0; r < repeats; r++) {
uint64_t t0 = now_ns();
for (int i = 0; i < N; i++) output[i] = fn(input[i]);
uint64_t dt = now_ns() - t0;
if (dt < best) best = dt;
}
double ns_per_op = (double)best / (double)N;
double max_rel = 0.0, sum_rel = 0.0;
for (int i = 0; i < N; i++) {
double ref = 1.0 / sqrt((double)input[i]);
double got = (double)fn(input[i]);
double rel = fabs(got - ref) / ref;
if (rel > max_rel) max_rel = rel;
sum_rel += rel;
}
printf("%-32s %8.3f ns/op max_rel=%.3e avg_rel=%.3e\n",
name, ns_per_op, max_rel, sum_rel / (double)N);
}
static volatile float sink;
int main(void) {
srand(42);
for (int i = 0; i < N; i++) {
double u = (double)rand() / (double)RAND_MAX;
double e = -2.0 + u * 6.0;
input[i] = (float)pow(10.0, e);
}
bench("Q_rsqrt (0x5f3759df, N=0)", q_rsqrt_quake_n0, 5);
bench("Q_rsqrt (0x5f3759df, N=1)", q_rsqrt_quake, 5);
bench("Q_rsqrt (0x5f3759df, N=2)", q_rsqrt_quake_n2, 5);
bench("Q_rsqrt (Lomont, N=0)", q_rsqrt_lomont_n0,5);
bench("Q_rsqrt (Lomont, N=1)", q_rsqrt_lomont, 5);
bench("1.0f / sqrtf(x)", std_rsqrt, 5);
#if defined(__aarch64__)
bench("NEON vrsqrte + 1 Newton", neon_rsqrt_1, 5);
bench("NEON vrsqrte + 2 Newton", neon_rsqrt_2, 5);
#endif
sink = output[N - 1];
(void)sink;
return 0;
}ARM側のビルドと実行は普通に。
clang -O2 -std=c11 bench.c -o bench -lm
./benchx86側はNEON部分の代わりに _mm_rsqrt_ss のラッパーを足し、-mavx2 -mfma を付けてWSL2のclang 18でビルド、taskset -c 0 ./bench で1コアにピン留めしてジッタを抑える。
M4で Q_rsqrt を動かしたら今でも 1/sqrtf より速かったのは事実だが、Zen 3で -O2 -fno-fast-math を踏むと差が12倍に開いて、-ffast-math を足すと逆にほぼ並ぶ、という形で「速い理由」がプラットフォームによって入れ替わるのが今回の収穫だった。同じコード、同じclangでも、ターゲットとフラグで結論が動くというのを数字で見られたのは良い実験。