9.3 SIMD化π計算

Gregory-Leibniz seriesと呼ばれる以下の数列を用いて円周率1を計算しよう。

\[ \frac{\pi}{4} = \frac{1}{1} - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \frac{1}{11} + \cdots \]

1回の除算を1回の乗算に置き換えると

\[ \frac{\pi}{4} = \frac{2}{1\times 3} + \frac{2}{5\times 7} + \frac{2}{9\times 11} + \cdots \]

あるいは、monadic ÷の使用による1文字削減を期待して、 \[ \frac{\pi}{8} = \frac{1}{1\times 3} + \frac{1}{5\times 7} + \frac{1}{9\times 11} + \cdots \]

となる。

なお、2行めの変換で項数は半分になる。以下ではこの式に基づいて項を数えることにする。

Iverson記法への書き換え

↕n0項目1項目2項目3項目一般化
各項の添字048124×↕n
分母×´1‿3×´5‿7×´9‿11×´13‿15{×´1‿3+𝕩}
項÷×´1‿3÷×´5‿7÷×´9‿11÷×´13‿15{÷×´1‿3+𝕩}

この項列を8×+´すればいいのだから、nを繰り返し上限として、以下のようになる。

pi ← 8×+´{÷×´1‿3+𝕩}¨↕n

ただしこのプログラム化には↕nが大変大きな配列となり実行不能になってしまうという空間計算量の問題がある。 そこでリスト上で計算を行い最後にFold '´'するという方法はやめて、（それまでの部分和を更新する）次状態計算関数を繰り返し実行するRepeat '⍟'を使った方法に切り替えることにする。

n ← 1000×1000×1000
pq ← ¯3‿¯1              # 行は引数として渡すのではなく、親環境中に保持
Term ← {𝕩+÷×´pq 4⊸+ ↩}  # 引数として部分和をとり更新された部分和を返す
pi ← 8×Term⍟n 0

→ Try it online⌨️

SIMD命令を使った高速化

完成したプログラムはnに関わらず実行できるがリストを使ってないためSIMD命令を使った高速化の余地がない。 そこで妥当な長さのリストを使って空間計算量の問題を起こすことなく高速化しよう。

n‿chunk ← ⟨1000×1000×1000,1000⟩
seed‿span ← ⟨3+4×↕chunk,4×chunk⟩
sum ← 0
Term ← {
  diff ← +´{÷×⟜(¯2⊸+)}¨𝕩+seed
  sum diff⊸+ ↩
  𝕩+span
}
Term⍟(n÷chunk)0
pi ← 8×sum

→ Try it online⌨️

BQNでは配列はimmutableなので計算結果は毎回割り当てられることになるが、その計算の基となるリストseedは再利用したいので親環境で定義している。

これでおよそ5倍の高速化になる。 しかし、実はこれはまだ十分に並列化できていない。5行目の明示的なEachをやめて以下のように変更する。

--- pi-simd0.bqn        2023-06-21 21:40:58
+++ pi-simd.bqn 2023-06-21 21:40:51
@@ -2,7 +2,7 @@
 seed‿span ← ⟨3+4×↕chunk,4×chunk⟩
 sum ← 0
 Term ← {
-  diff ← +´{÷×⟜(¯2⊸+)}¨𝕩+seed
+  diff ← +´÷×⟜(¯2⊸+)𝕩+seed
   sum diff⊸+ ↩
   𝕩+span
 }

すると更に8倍速くなる。下のプログラムから40倍以上の高速化を実現できたということでrayonを使ったrustプログラム2にわずか2倍程度しか遅くないという速度を得ることができた。 Foldはどうしても逐次処理になってしまうのでできるだけ避けるべきであるというのはわかるが、なぜだかrank多相を使ってEachも避ける方がよいという結論になってしまった。

$ time cbqn pi-simple.bqn
Error: Out of memory
pi-simple.bqn:2:
  •Show pi ← 8×+´{÷×´1‿3+𝕩}¨↕n
                            ^
cbqn pi-simple.bqn  0.01s user 0.00s system 87% cpu 0.015 total
$ time cbqn pi-repeat.bqn
3.1415926445762157
cbqn pi-repeat.bqn  92.21s user 0.06s system 99% cpu 1:32.43 total
$ time cbqn pi-simd0.bqn
3.1415926530824487
cbqn pi-simd0.bqn  17.73s user 0.01s system 99% cpu 17.750 total
$ time cbqn pi-simd.bqn
3.1415926530824487
cbqn pi-simd.bqn  2.58s user 0.01s system 99% cpu 2.587 total
2

当然--release付き

golf

上記のプログラムを短くしよう。 元々の実行可能なプログラムが4行だったので対比しやすいようそれに合わせて書き換えると以下のようになる。

n‿chunk ← ⟨1000×1000×1000,1000⟩
span‿seeds ← ⟨chunk×4,3+4×↕chunk⟩
Term ← {⟨𝕨++´÷×⟜(¯2⊸+)𝕩+seeds,𝕩+span⟩}     # 添字と部分和のペアに対する繰り返し実行
pi ← 8×⊑Term´⍟(n÷chunk)0‿0

→ Try it online⌨️

これでプレゼンテーションの1枚のシートに上下並べて2つのバージョンを載せることができるようになった。 APL（語族）を使うならこうでなければ。