Rust 1.98 の algebraic_* でオーディオ DSP は速くなるのか検証してみた

はじめに

オーディオ DSP では、サンプルごとの浮動小数点演算を短い時間内に何度も繰り返します。通常の Rust コードが CPU の SIMD 性能をどこまで使えているのか分かりにくく、CPU 固有の intrinsics へ進むべきか迷うことがあります。

Rust 1.98 では、f32 と f64 へ algebraic_* 系のメソッドが追加されました。 浮動小数点演算の順序などに関する制約を緩め、コンパイラーがより積極的に最適化できる API です。標準 Rust のまま積和演算を速くできる可能性がある一方、通常演算と同じ数値になるとは限りません。

そこで、gain・mixing・除算・dot product・64 タップ FIR で挙動を比較しました。結果、gain と mixing は通常演算とほぼ同じでしたが、dot product は 2.751〜6.936 倍となり、FIR は 3.466〜3.513 倍となりました。一方で、高速化と引き換えに通常演算との差が生じました。

algebraic_* とは

Rust 1.98 で stable になったメソッドは次の 5 種類です。

• algebraic_add
• algebraic_sub
• algebraic_mul
• algebraic_div
• algebraic_rem

通常の浮動小数点加算では、a + b + c + d は ((a + b) + c) + d の順で評価されます。各演算で丸められるため、加算順序を変えると結果も変わることがあります。

algebraic_add をつなげた場合は、コンパイラーが (a + b) + (c + d) のように部分和を同時に求められます。公式ドキュメント では、演算の結合や並べ替え、除算と逆数乗算の変換、負のゼロを区別しない最適化などが許可されています。最適化の厳密な集合は規定されておらず、同じ入力でも結果が変わる可能性があります。 (ただし、未定義動作にはならない。)

検証環境

• MacBook Pro、Apple M4、10 コア、16 GB メモリー
• macOS 26.6.2、arm64
• rustc 1.98.1、LLVM 22.1.8
• Rust 2024 Edition、opt-level = 3
• codegen-units = 1、target-cpu=apple-m4
• 外部クレートなし

対象読者

• Rust で数値計算やオーディオ DSP を実装している方
• SIMD 化のために CPU 固有 intrinsics を書くか迷っている方
• algebraic_* の性能と、通常演算との結果の違いを知りたい方

参考資料

• Announcing Rust 1.98.0
• Primitive Type f32 の Algebraic operators

検証方法

性質の異なる 5 種類の処理を用意しました。

• gain: 各サンプルへ同じ係数を乗算する
• mixing: 2 信号の対応するサンプルを加算する
• 除算: 各サンプルを同じ値で除算する
• dot product: 2 配列の要素を乗算し、1 つの値へ加算する
• 64 タップ FIR: 各出力サンプルについて 64 回の積和演算を行う

通常演算版と algebraic_* 版は、演算以外のループ構造をそろえました。除算には、除数の逆数を一度求めて通常の乗算を行う実装も加えました。

64・128・256・512 frames を対象に、1 試行が 50 ms 以上になる反復数を求めた後、実行順を変えながら各条件を 15 回測りました。掲載値は 1 関数呼び出し当たりの中央値です。

数値差については、固定シードから生成した [-1, 1] の有限値を使い、通常演算・algebraic_*・f64 で同じ積和を計算した値を比較しました。dot product では各サイズ 10,000 件を、FIR では各サイズ 1,000 ブロックを実行しました。

実装

違いが最も分かりやすい dot product の本質部分は次のとおりです。

fn dot_normal(left: &[f32], right: &[f32]) -> f32 {
    left.iter()
        .zip(right)
        .fold(0.0, |sum, (&left_sample, &right_sample)| {
            sum + left_sample * right_sample
        })
}

fn dot_algebraic(left: &[f32], right: &[f32]) -> f32 {
    left.iter()
        .zip(right)
        .fold(0.0, |sum, (&left_sample, &right_sample)| {
            sum.algebraic_add(left_sample.algebraic_mul(right_sample))
        })
}

FIR も入力と係数を zip し、同じ形の fold で 64 要素を加算しました。入力範囲はタップループの前にスライス化し、内側へ境界検査を残さないようにしています。

結果

速度比は 通常演算 / algebraic_* です。1 より大きいほど algebraic_* が速いことを示します。

gain と mixing においては 0.980〜1.007 倍であり、性能差はほとんどありませんでした。どちらも要素間の依存関係がなく、通常演算版も algebraic_* 版も同種の NEON 命令になりました。通常演算のままでもベクトル化できる処理です。

dot product は要素数が増えるほど速度比が大きくなり、512 要素では通常演算版の 174.995 ns に対して algebraic_* 版は 25.231 ns でした。64 タップ FIR はフレーム数によらず約 3.5 倍で、512 frames では 5.024 µs から 1.430 µs になりました。

除算については、64 要素では差がありませんでしたが、512 要素では 1.361 倍になりました。アセンブリを見ると、通常演算版ではループ内でベクトル除算を実行し、algebraic_div 版では逆数を求めてからベクトル乗算を実行していました。手動の逆数乗算版では 512 要素で 1.398 倍となり、同じ傾向を示しました。

dot product のアセンブリには、次の違いがありました。通常演算版は 4 要素をまとめて乗算した後、各要素をスカラーレジスターへ取り出して順番に加算しています。

fmul.4s  v1, v1, v5
mov      s18, v1[1]
fadd     s0, s0, s1
fadd     s0, s0, s18

algebraic_* 版では、4 本のベクトル部分和を fmla.4s で並行して更新し、最後に水平加算していました。

fmla.4s  v0, v16, v4
fmla.4s  v1, v17, v5
fmla.4s  v2, v18, v6
fmla.4s  v3, v19, v7
fadd.4s  v0, v1, v0
fadd.4s  v1, v3, v2
fadd.4s  v0, v1, v0
faddp.4s v0, v0, v0
faddp.2s s0, v0

LLVM IR では、algebraic_* 側に reassoc nsz arcp contract が付きました。nnan と ninf は付いておらず、すべての fast-math フラグを有効にする実装ではありませんでした。

性能と引き換えに、通常演算との数値差が生じました。dot product の結果は次のとおりです。

「最大絶対誤差」は、入力を f64 へ拡張して計算した参照値との差です。今回の入力では algebraic_* 版の方が参照値に近くなりました。

64 タップ FIR の通常演算との最大絶対差は、全サイズで 5.215406e-8〜6.705523e-8 でした。除算 10,000 件では、通常除算と algebraic_div の差は最大 1 ULP、絶対値で 5.960464e-8 でした。

演算順序とゼロの符号による差も確認できました。

有限値だけの加算でも、再結合によって 1.0 と 4.0 に分かれました。中間結果がオーバーフローする加算では、通常演算と algebraic_* が Infinity と NaN に分かれました。負のゼロは数値比較では正のゼロと等しいものの、ビット列は異なりました。

考察

gain や mixing のように各要素を独立して処理する場合は、最初から algebraic_* へ置き換えず、通常ループが自動ベクトル化されたかを確認する方法で十分だと思います。今回のように同種の命令が生成されるなら、計算結果の違いを受け入れてまで置き換える必要はないでしょう。

一方、dot product や FIR では、一つ前の加算結果を次の加算に使うため、通常演算では加算を順番に進めます。順序を変えてよい処理なら、複数の部分和を並行して計算できるため、試す価値がありそうです。ただし、境界検査やループの形がベクトル化を妨げる場合もあります。コンパイラーが出力した命令と、代表的なブロックサイズでの処理時間を確認する必要があります。

共通の値による除算では、逆数乗算へ変換したことで、要素数が多い条件の処理時間が短くなりました。除数がループ中で変わらないことをコードで表せるなら、手動で逆数を一度求める実装も比較対象になります。どちらを選ぶ場合も、1 回ごとの除算と、先に丸めた逆数の乗算が異なる結果を返してよいか決める必要がありそうです。

リアルタイムオーディオ処理では、演算順序や丸め結果の違いを許容できる積和処理を候補とし、対象 CPU でコールバックの処理時間と出力信号を確認してから採用するのがよいでしょう。

まとめ

Rust 1.98 の algebraic_* を使い、オーディオ DSP を想定した 5 種類の処理を比較しました。通常演算でもベクトル化できた gain と mixing はほぼ同じ性能でした。一方、積和を集約する dot product は 2.751〜6.936 倍、64 タップ FIR は 3.466〜3.513 倍になりました。

algebraic_* は、すべての浮動小数点処理を速くする API ではありません。演算順序や丸め結果の違いを許容できる積和処理へ絞り、代表的な入力での結果・コンパイラーが出力した命令・処理時間・誤差を一緒に確認するのがよさそうです。