- Common LispのSBCLには
sb-simdというSIMDインターフェースが組み込まれており、AVX命令を使うとsingle-floatを8個まとめて1命令で処理できます。 - typed版を基準にSSE（4要素）で1.4倍、AVX（8要素）で1.5倍程度の高速化と、メモリ帯域もボトルネックなため、演算幅を単純に倍にしても効果が頭打ちになりました。
- ただ、型宣言なしのコードに
(declare (type ...))を加えるだけで約20倍も速度が違うので、SIMDよりも、まずは型宣言の効果が大きかったです。 - SIMD pack操作には
SB-SIMD-SSE::F32.4+のような内部パッケージの関数を使う必要があり、do-all-symbolsで探さないと見つかりませんでした。
(defun saxpy/avx8 (a x y)
(declare (type single-float a)
(type (simple-array single-float (*)) x y))
(let* ((n (length x))
(step 8)
(main-end (* step (floor n step))))
(declare (type fixnum n step main-end))
(loop for i of-type fixnum from 0 below main-end by 8
do (setf (sb-simd-avx::f32.8-row-major-aref y i)
(sb-simd-avx::f32.8+
(sb-simd-avx::f32.8* a (sb-simd-avx::f32.8-row-major-aref x i))
(sb-simd-avx::f32.8-row-major-aref y i))))
(loop for i of-type fixnum from main-end below n
do (setf (aref y i) (+ (* a (aref x i)) (aref y i))))
y))
Code language: Lisp (lisp)1. SIMDとsb-simdパッケージ
SBCLには sb-simd というSIMDインターフェースが組み込まれています。
ELS 2022の論文 “Closing the Performance Gap Between Lisp and C”（Marco HeisigとHarald Köstler）で紹介されたもので、数値計算のベンチマークで「最適化されたCコードの最大94%の性能」が達成されたと報告されています1。
「SIMD」は、1つの命令で複数のデータを同時に処理するCPUの仕組みです。
Single Instruction Multiple Dataの略で、MMXやSSE、AVX、NEONなどの命令セットの総称です。
マルチメディア処理（音声・動画・画像）の需要が急増した2000年前後から、CPU命令セットの拡張として実装されてきました。
SIMD（概念）
└── x86実装
├── MMX (1997) 64bit 整数のみ
├── SSE (1999) 128bit 単精度float
├── SSE2 (2001) 128bit 倍精度floatも
├── SSE3/4 (2004〜) 128bit 命令追加
├── AVX (2011) 256bit
└── AVX-512 (2016) 512bit
└── ARM実装
└── NEON (2004〜) 128bit
└── その他
└── AltiVec（PowerPC）などたとえば、通常の加算は1要素ずつ処理しますが、AVX（Advanced Vector eXtension）を使うと single-float を8個まとめて1命令で加算できます。
1.1. SIMDが効く条件
効果が出る条件はかなり限定されています。
- 型が固定された数値の連続配列（
(simple-array single-float (*))） - 各要素の計算が互いに独立している
- 十分に大きいデータ量（数千〜数百万要素）
- 分岐が少ない
逆に、グラフ探索、ポインタをたどる構造、型が混在するデータ、小さいNには向きません。
1.2. SAXPY（Single-precision A·X Plus Y）
今回の比較対象として選んだのは「SAXPY（Single-precision A·X Plus Y）」です。
y[i] = a * x[i] + y[i]a がスカラー、x と y がベクトルです。
SAXPYは、昔から使われている行列・ベクトル計算の基本部品で、BLAS（Basic Linear Algebra Subprograms）の Level 1 ルーチンの一つです2。
「Single-precision」は、32ビット浮動小数点版の名前で、double-float 版は DAXPY と呼ばれます。
SAXPYは、機械学習の重み更新、音声波形の音量調整、画像の明るさ補正など、「同じ線形演算を配列全体に適用する」処理の代表です。
2. :sb-simdパッケージを読み込む
sb-simdを使うには、:sb-simdパッケージを読み込みます。
その後、:saxpy-benchという自分のパッケージ内にコードを書きました。
(eval-when (:compile-toplevel :load-toplevel :execute)
(require :sb-simd))
(defpackage :saxpy-bench
(:use :cl))
(in-package :saxpy-bench)Code language: Lisp (lisp)SLIMEで C-c C-k（ファイル全体のコンパイル）を使うとき、2つの落とし穴がありました。
1つ目は require の問題です。(require :sb-simd) をトップレベルに書いても、コンパイラがファイルを読み始める時点ではまだ実行されていません。
リーダーが sb-simd-avx::f32.8-row-major-aref という名前を読もうとした瞬間に「SB-SIMD-AVX パッケージが存在しない」というエラーになります。
解決策は eval-when を使うことです。
(eval-when (:compile-toplevel :load-toplevel :execute)
(require :sb-simd))Code language: Lisp (lisp):compile-toplevel を指定することで、コンパイル時のリードフェーズの前に sb-simd が読み込まれます3。
2つ目は delete-package の問題です。
REPLで試行錯誤しているときに「パッケージが既に存在します」という警告を避けるために書いていましたが、コンパイル時に走るとパッケージを定義した直後に消してしまい、その後の in-package が失敗します。
ファイルには書かないのが正解です。
また、C-c C-k 後にREPLのプロンプトはそのままです。SAXPY-DEMO> に移動するには、REPLで (in-package :saxpy-demo) を打ちます。
2.1. SB-SIMD-SSE::F32.4-ROW-MAJOR-AREF
sb-simd は元々独立リポジトリとして開発されていましたが、現在はSBCL本体に統合されています4。
たとえば、x86-64環境（Intel/AMD）であればSSEやAVXを使えますが、Apple Silicon（M1/M2/M3）では、AVX系の最適化はそのまま期待できません。
(require :sb-simd) が通るかどうかは環境次第なので、まず試してみるのが確実です。
pack版は SB-SIMD-SSE:: や SB-SIMD-AVX:: という内部パッケージに入っているため、do-all-symbols で探す必要があります。
(do-all-symbols (s)
(let ((name (symbol-name s))
(pkg (symbol-package s)))
(when (and pkg
(search "SB-SIMD" (package-name pkg))
(search "F32." name))
(format t "~&~A::~A" (package-name pkg) name))))Code language: Lisp (lisp)こうすると SB-SIMD-SSE::F32.4-ROW-MAJOR-AREF や SB-SIMD-AVX::F32.8-ROW-MAJOR-AREF が出てきます。
名前の形は 型.要素数 です。F32.4 は「single-float 4個をまとめたSIMD pack」、F32.8 は8個まとめた版です。
このような X.Y 形式の関数が、実際にCPUのSIMDレジスタを使う関数です。
(defun saxpy/sse4 (a x y)
(declare (type single-float a)
(type (simple-array single-float (*)) x y))
(let* ((n (length x))
(step 4)
(main-end (* step (floor n step))))
(declare (type fixnum n step main-end))
(loop for i of-type fixnum from 0 below main-end by 4
do (setf (sb-simd-sse::f32.4-row-major-aref y i)
(sb-simd-sse::f32.4+
(sb-simd-sse::f32.4* a (sb-simd-sse::f32.4-row-major-aref x i))
(sb-simd-sse::f32.4-row-major-aref y i))))
(loop for i of-type fixnum from main-end below n
do (setf (aref y i) (+ (* a (aref x i)) (aref y i))))
y))
Code language: PHP (php)ただし、SB-SIMD-SSE:: や SB-SIMD-AVX:: という内部パッケージにあります5。
公開APIとして安定して使える保証はなく、SBCL実験寄りの低水準コードとして扱うのが安全です。
2.2. 最終的なコード
(eval-when (:compile-toplevel :load-toplevel :execute)
(require :sb-simd))
(defpackage :saxpy-bench
(:use :cl))
(in-package :saxpy-bench)
(declaim (optimize (speed 3) (safety 0) (debug 0)))
;;; ----------------------------------------------------------------
;;; ユーティリティ
;;; ----------------------------------------------------------------
(defun make-f32 (n &optional (init 0.0f0))
(declare (type fixnum n) (type single-float init))
(make-array n :element-type 'single-float :initial-element init))
(defun fill-test-data (x y)
(declare (type (simple-array single-float (*)) x y))
(dotimes (i (length x))
(setf (aref x i) (coerce (/ (mod i 100) 100.0) 'single-float)
(aref y i) 1.0f0))
(values x y))
(defun copy-f32 (src)
(declare (type (simple-array single-float (*)) src))
(let ((dst (make-f32 (length src))))
(dotimes (i (length src) dst)
(setf (aref dst i) (aref src i)))))
(defun same-p (a b &optional (eps 1.0e-5))
(declare (type (simple-array single-float (*)) a b)
(type single-float eps))
(loop for i of-type fixnum below (length a)
always (< (abs (- (aref a i) (aref b i))) eps)))
;;; ----------------------------------------------------------------
;;; SAXPY バリアント
;;; ----------------------------------------------------------------
;;; 1. naive（型宣言なし）
(defun saxpy/untyped (a x y)
(dotimes (i (length x) y)
(setf (aref y i) (+ (* a (aref x i)) (aref y i)))))
;;; 2. naive（型宣言あり）
(defun saxpy/typed (a x y)
(declare (type single-float a)
(type (simple-array single-float (*)) x y))
(dotimes (i (length x) y)
(declare (type fixnum i))
(setf (aref y i) (+ (* a (aref x i)) (aref y i)))))
;;; 3. SSE 4要素版
(defun saxpy/sse4 (a x y)
(declare (type single-float a)
(type (simple-array single-float (*)) x y))
(let* ((n (length x))
(step 4)
(main-end (* step (floor n step))))
(declare (type fixnum n step main-end))
(loop for i of-type fixnum from 0 below main-end by 4
do (setf (sb-simd-sse::f32.4-row-major-aref y i)
(sb-simd-sse::f32.4+
(sb-simd-sse::f32.4* a (sb-simd-sse::f32.4-row-major-aref x i))
(sb-simd-sse::f32.4-row-major-aref y i))))
(loop for i of-type fixnum from main-end below n
do (setf (aref y i) (+ (* a (aref x i)) (aref y i))))
y))
;;; 4. AVX 8要素版
(defun saxpy/avx8 (a x y)
(declare (type single-float a)
(type (simple-array single-float (*)) x y))
(let* ((n (length x))
(step 8)
(main-end (* step (floor n step))))
(declare (type fixnum n step main-end))
(loop for i of-type fixnum from 0 below main-end by 8
do (setf (sb-simd-avx::f32.8-row-major-aref y i)
(sb-simd-avx::f32.8+
(sb-simd-avx::f32.8* a (sb-simd-avx::f32.8-row-major-aref x i))
(sb-simd-avx::f32.8-row-major-aref y i))))
(loop for i of-type fixnum from main-end below n
do (setf (aref y i) (+ (* a (aref x i)) (aref y i))))
y))
;;; ----------------------------------------------------------------
;;; 計測ヘルパー
;;; ----------------------------------------------------------------
(defun measure-ms (fn x y-init a runs)
"fn を runs 回実行し、最小実行時間をミリ秒で返す。"
(declare (type fixnum runs))
(let ((best most-positive-double-float))
(dotimes (_ runs best)
(let ((y (copy-f32 y-init)))
(let ((start (get-internal-real-time)))
(funcall fn a x y)
(let ((elapsed (* 1000.0d0
(/ (- (get-internal-real-time) start)
internal-time-units-per-second))))
(when (< elapsed best)
(setf best elapsed))))))))
;;; ----------------------------------------------------------------
;;; ベンチマーク実行
;;; ----------------------------------------------------------------
(defun print-disassembly ()
(dolist (fn '(saxpy/typed saxpy/sse4 saxpy/avx8))
(format t "~%;;; ============================================================~%")
(format t ";;; disassembly: ~A~%" fn)
(format t ";;; ============================================================~%")
(let ((asm (with-output-to-string (*standard-output*)
(disassemble (symbol-function fn)))))
(write-string asm)
(terpri))))
(defun run-bench (&key (n 10000000) (runs 5))
(let* ((x (make-f32 n))
(y0 (make-f32 n))
(a 0.01f0))
(fill-test-data x y0)
;; 結果の正しさチェック用（typed を基準）
(let ((ref (copy-f32 y0)))
(saxpy/typed a x ref)
(format t "~%SAXPY benchmark n=~:D runs=~D (best of)~%" n runs)
(format t "~60A ~8A ~6A ~6A~%"
"variant" "ms" "ratio" "ok?")
(format t "~60A ~8A ~6A ~6A~%"
(make-string 56 :initial-element #\-)
"--------" "------" "------")
(let (baseline)
(dolist (entry
(list (list "naive / untyped" #'saxpy/untyped)
(list "naive / typed (declaim)" #'saxpy/typed)
(list "SIMD / SSE 4-wide" #'saxpy/sse4)
(list "SIMD / AVX 8-wide" #'saxpy/avx8)))
(let* ((label (first entry))
(fn (second entry))
(ms (measure-ms fn x y0 a runs))
(ratio (if baseline (/ baseline ms) 1.0d0))
(ok (let ((y (copy-f32 y0)))
(funcall fn a x y)
(same-p ref y))))
(unless baseline (setf baseline ms))
(format t "~60A ~8,3F ~6,2Fx ~6A~%"
label ms ratio (if ok "OK" "FAIL"))))))))
#-swank(print-disassembly)
#-swank(run-bench)
Code language: Lisp (lisp)3. 実行結果（約1.5倍）
4つのバリアントを比較しました。
naive / typed をベースライン（1.00x）として、手元のMacBook Air（Intel 2018）と AtCoder のジャッジサーバー両方で計測しています。
SAXPY benchmark n=10,000,000 runs=5 (best of)
variant ms ratio ok?
-------------------------------------------------------- -------- ------ ------
naive / untyped 178.473 1.00x OK
naive / typed (declaim) 8.526 20.93x OK
SIMD / SSE 4-wide 5.993 29.78x OK
SIMD / AVX 8-wide 5.385 33.14x OKMacBook Air Intel 2018、n=10,000,000、runs=5 の最小値
SAXPY benchmark n=10,000,000 runs=5 (best of)
variant ms ratio ok?
-------------------------------------------------------- -------- ------ ------
naive / untyped 129.001 1.00x OK
naive / typed (declaim) 8.000 16.13x OK
SIMD / SSE 4-wide 6.000 21.50x OK
SIMD / AVX 8-wide 5.000 25.80x OK
AtCoder ジャッジサーバー、n=10,000,000、runs=5 の最小値
3.1. 読み取れること
- untyped は型宣言だけで20倍遅い。
- typed と SSE 4-wide の差は約1.4倍。
- SSE → AVX の差は1.1倍前後で頭打ち。
まずは、もっとも効率に影響していたのは、「型宣言」でした。
(declare (type single-float a)
(type (simple-array single-float (*)) x y))Code language: Lisp (lisp)ボクシング（single-float をヒープオブジェクトとして包む処理）とアンボクシングのコストが、ループ全体に乗ります。
数値計算で (declaim (type ...)) を書く意味がここに出ています。
ただ、型宣言を与えてもSBCLは自動ベクトル化せず、スカラー命令（MULSS / ADDSS）にとどまります。
そこで、sb-simd を明示的に呼ぶことで初めてpacked命令（MULPS / ADDPS）が出て、実測でも差がつきます。
しかし、SIMD幅を128bitから256bitに倍にしても、速度は1.1倍程度しか伸びませんでした。
SAXPYは y[i] = a * x[i] + y[i] という軽い演算に対してメモリアクセスが多い処理で、ボトルネックがCPUの演算ユニットではなくメモリ帯域にあります。
演算が速くなっても、メモリの読み書き速度が律速になるため効果が頭打ちになります。
両環境でバリアント間の比率がほぼ一致していることから、この構造はSBCLの最適化とSAXPYの性質によるもので、環境依存ではないと言えます。
3.2. アセンブリで確認する
disassemble で実際に出力される命令を確認しました。
saxpy/typed（スカラー）
L0: F30F104C4F01 MOVSS XMM1, [RDI+RCX*2+1]
F30F59CB MULSS XMM1, XMM3
F30F10544E01 MOVSS XMM2, [RSI+RCX*2+1]
F30F58D1 ADDSS XMM2, XMM1
F30F11544E01 MOVSS [RSI+RCX*2+1], XMM2Code language: CSS (css)MULSS / ADDSS はスカラー命令で、1要素ずつ処理します。(declaim (optimize (speed 3))) と型宣言を与えても、SBCLは自動ベクトル化しません。
saxpy/sse4（SSE packed）
L0: 0F105C4F01 MOVUPS XMM3, [RDI+RCX*2+1]
0F59CB MULPS XMM1, XMM3
410F105C4801 MOVUPS XMM3, [R8+RCX*2+1]
0F58CB ADDPS XMM1, XMM3
410F114C4801 MOVUPS [R8+RCX*2+1], XMM1Code language: CSS (css)MULPS / ADDPS は128bit packed命令で、single-float を4個まとめて処理します。MOVUPS（Unaligned）を使っているのはSBCLのヒープ配列がアライメント保証を持たないためで、適切な選択です。
saxpy/avx8（AVX packed）
L0: C4E27D18C1 VBROADCASTSS YMM0, XMM1
C5FC10544F01 VMOVUPS YMM2, [RDI+RCX*2+1]
C5FC59C2 VMULPS YMM0, YMM0, YMM2
C4C17C10544801 VMOVUPS YMM2, [R8+RCX*2+1]
C5FC58C2 VADDPS YMM0, YMM0, YMM2
C4C17C11444801 VMOVUPS [R8+RCX*2+1], YMM0Code language: CSS (css)VMULPS / VADDPS は256bit packed命令で、single-float を8個まとめて処理します。
ループ先頭の VBROADCASTSS はスカラー a を256bitレジスタの全8レーンに複製する命令で、SSE版の SHUFPS に相当します。
両環境でアセンブリ出力が一致しており、コードが意図通りの命令に落ちていることが確認できました。
4. 学んだこと
sb-simd を使いはじめて最初に引っかかりやすいのは、「:sb-simd パッケージの公開APIを使えばSIMDになる」という誤解です。
実際には：
do-external-symbolsで見えるSB-SIMD:F32+はスカラー演算（SIMD packではない）- 本当のSIMD pack操作は
SB-SIMD-SSE::F32.4+やSB-SIMD-AVX::F32.8+のようなX.Y形式の関数 - これらは
do-all-symbolsを使わないと見つからない内部パッケージにある X.YのYが要素数で、ループの刻みと一致させる必要があるC-c C-kで使うにはeval-whenが必要で、delete-packageは書かない
また、SIMDを使えば必ず速くなるわけではなく、処理のボトルネックがどこにあるかで効果が大きく変わります。
- 論文はErlangen-Nürnberg大学（FAU）の2名による発表で、2022年3月21〜22日にポルトガル・ポルトで開催された第15回European Lisp Symposium（ELS 2022）の予稿集に収録されています。DOI: 10.5281/zenodo.6335627。 – Closing the Performance Gap Between Lisp and C（Zenodo）
- BLASはLevel 1（ベクトル演算）、Level 2（行列ベクトル演算）、Level 3（行列行列演算）の3段階に分かれています。Level 1は1970年代初頭に定義され、当時のハードウェアではこれで十分なピーク性能を出せましたが、その後のマシンではLevel 2（1984〜86年）、Level 3（1987〜88年）が追加されました。SAXPYはLevel 1の代表的な演算です。 – Basic Linear Algebra Subprograms（Wikipedia）
eval-whenはANSI Common Lispの標準スペシャルオペレータです。:compile-toplevelはファイルコンパイラがトップレベルフォームを処理する際にコンパイル時に評価、:load-toplevelはFASLのロード時に評価、:executeはそれ以外のコンテキスト（REPLなど）で評価することを指定します。この3つを全て指定するのが、どのコンテキストでも確実に実行させる標準的な書き方です。 – EVAL-WHEN（Common Lisp HyperSpec）- 元々のGitHubリポジトリ（marcoheisig/sb-simd）はアーカイブ済みで、「コードはSBCLにマージされました」と記載されています。現在のコードは
https://github.com/sbcl/sbcl/tree/master/contrib/sb-simdにあります。SBCLのメーリングリストでは、Quicklisp経由での配布ではSBCL内部の変更で壊れやすいという理由から、contrib化が議論・決定されました。 – marcoheisig/sb-simd（GitHub） SB-SIMD-SSEやSB-SIMD-AVXは命令セットごとに分けられた内部パッケージです。SBCLメーリングリストでの議論では、これらはSBCLの内部変更に依存するため、外部公開APIとしての安定性は保証されていないとされています。 – Potential contrib: sb-simd（sbcl-devel）