前回

前回は、2300個の処理を2300bitベクタへまとめ、カウンタもbit-sliced化しました。

さらに、最大700万まで数えられれば十分なので、カウンタ幅を32bitから23bitへ削減。

最大ケースの予測時間は、

約8時間50分
    ↓ 
約33分44秒

まで縮まりました。

ここまでは主に、データの持ち方をIcarus向けに変える高速化です。

今回はもう一歩進めて、実際に行うカウンタ更新そのものの回数を減らす方法を考えてみます。

※ 今回のコードや実装資料はGithubに上げてあります。文末のリンクからご参照ください。
※ 本記事のテストや測定はすべて筆者のローカル環境(Win 11, i7-12700H, MEM 64GB)で行っています。

高速化⑤: カウンタ加算の回数を減らす - 2-rise fusion

現在は、2300レーンが新しい入力を受けるたびにカウンタの更新をしています。

しかし、一つのレーンに連続して入ってくる入力をよく見てみると、2回連続で立ち上がりが検出されることはないことに気づきます。

ある入力で立ち上がりが検出されるのは、0/1列が1になったときです。

次の立ち上がりを見るためには、どんなに急いでも、0に下げてからまた1に上げる、という2回の入力が必要です。

だったら、2回分の入力をまとめて1回のカウンタ更新にできるんじゃない？というアイデアが生まれます。

どうやって実装する？

奇数番めの入力ではカウンタ更新を行わず、出力されたrise信号をそのままdelayレジスタに保存しておきます。

次の入力で受け取ったrise信号とdelay信号を各bitごとに見比べると、(0,0), (0,1), (1,0)のどれかになります。

(1,1)は存在しません。

ですから

pair_event = rise OR delay

で2入力ぶんの合計カウントアップ値を一つの2300bit信号にまとめることができます。

これで、入力ごとに23回のbit演算から2入力ごとに23+1回のbit演算まで計算回数が下がるはずです。

結果発表：2-rise fusionは効いた？

では結果です。

ここまで積み重ねてきた高速化①〜④に、今回の2-rise fusionを追加して比較します。

まず K=2300, N=10,000。

高速化④まで        2.819秒
2-rise fusion追加後 1.522秒

約1.85倍速くなりました。

2入力ぶんの立ち上がりを1回のカウンタ更新にまとめた効果が、かなり素直に出ています。

N=100,000でも、

高速化④まで         28.911秒
2-rise fusion追加後 14.861秒

まで短縮しました。

では最大ケース N=7,000,000 はどうなりそうでしょうか。

今回も N=100,000 の実測値から単純に線形外挿すると、

高速化④まで         約33分44秒
2-rise fusion追加後 約17分20秒

さらに半分近くまで縮まりました。

ここまでの最大ケース予測を並べると、

最初の2300レーン版  約8時間50分
高速化①＋②          約3時間23分
高速化③             約47分56秒
高速化④             約33分44秒
高速化⑤             約17分20秒

です。

もちろん今回も、それまで採用した高速化は全部そのままです。

解法そのものは、まだ愚直なまま。

それでも17分台。

Icarus君、だいぶ空を飛ぶようになってきました😁

高速化⑥: 16回の更新を1回にする - Carry Save Adder

2-rise fusionで、カウンタ更新は2入力に1回に減りました。

もっとカウンタ更新を減らす方法はないでしょうか。

考えてみましょう。

16回ぶんの更新イベントをまとめる

今回は、pair_eventを16回ぶん貯めてからカウンタへ反映する方法を検討してみます。

うまくいけばかなりの高速化が実現するかもしれません。

今回16回を選んだのは、5bitで表せて扱いやすいかな、と思ったからです。

32回や64回をまとめる構成も考えられますが、今回は比較していません。

まずは16回ぶんで試してみます。

まずはアイデアを考える

これまではイベント1回ごとに「1の位に1を足すか足さないか」を表す2300bitのpair_event出力を作り、23bit加算器に渡していました。

しかし、前回記事のbit slice加算器の構造を考えると、+1にこだわる必然はないような気がします。

例えば2の位に+1をするかしないか、4の位に+1をするかしないか、というbitsetを作れるとします。

2の位の計算時には1の位からの繰上りの加算をしてから、その結果にさらに新しく作る2の位のbitsetを加算、という順序で計算ができるはずです。

ということは、加算したい数字を2進数で持てれば、+3でも+7でも好きな数字の加算がちょっとした計算量の追加でできそうです。

具体的な持ち方を考える

例えば一番最初のレーンだけを見たときに、16回のpair_eventが、

1 0 1 0 1 0 0 0 0 0 0 0 0 1 0 1

だったとします。

1は5個あります。つまりこのレーンは+5をバイナリ表現できればいいのです。

早い話が、先頭レーンの加算値は +5、2進数なら 5'b00101 だと示せればOKです。

ですから、

1の位の2300bitset  [1, ...]
2の位の2300bitset  [0, ...]
4の位の2300bitset  [1, ...]
8の位の2300bitset  [0, ...]
16の位の2300bitset [0, ...]

のように表現できればこの目論見は成功しそうです。

どうやって実現する？

必要になる変換は、16本のパラレルな0/1入力を、その1の個数を表すバイナリ値へ変換する処理です。

これを3:2 Compressor（3:2圧縮器）と呼ばれる回路を使って実現してみます。

なんだかマッチョな名前ですね。モーターの振動が伝わってきそうです😁

3:2 Compressorってなんだ？

入力がパラレル3本、出力がバイナリ2本の加算器です。

3本の入力端子のうち、1になっている本数を2進数2桁の出力として表現します。

例えば入力が(1,1,1)であれば出力は3の2進数表現、つまり(1,1)を出力する仕組みです。

名前はごついですが、役割はかなりシンプルです。

具体的にどうするんだ？

まず仕組みを追いやすくするため、2300レーンのうち先頭レーンだけを抜き出して考えます。

実際の回路では、このあと説明する処理を2300bitのbitsetに対してまとめて行います。

先ほどの例では、先頭レーンの16回イベントぶんの入力が C0 = 1010100000000101です。

1の個数は最大16回ですから5bitで表現できますね。

ここでは、同じ重みのbitを一時的に置いておく列として、C0_2, C0_4, C0_8, C0_16も用意しておきます。

まずは、C0の下位3bitをCompressor入力にコピーします。

そしてC0を3bit右シフトします。

C0   = 1010100000000101
      ↓
comp = 101
C0   = 1010100000000

次にCompressor出力をC0とC0_2の位に戻します。

Compressor入力は101ですから、出力は2の2進数表現、10です。

出力の1の位をC0に、2の位をC0_2に戻します。

戻し方は C0, C0_2をそれぞれ1bit左シフトしてからORで再下位桁に加算です。

comp_out = 10
      ↓
C0   = 10101000000000
C0_2 = 1

さて、今の処理で何が変わったでしょうか。

答えは「C0～C0_16が持つイベント結果ビット数合計が16から15に減った」です。

入力用に3bitを取り除いてから、出力の2bitを戻していますから、一度の処理で1bitずつ合計bit数が減っていきます。

ということは、16bitパラレルのデータを5bitバイナリに圧縮するには、今の処理を11回程度繰り返せばいいことになります。

具体的には、C0をどんどん処理して、最後の1bitまで短縮したら、次はC0_2を減らしながら、C0_4とC0_2のビット列に出力を追加、終わったら、さらにC0_4、C0_8を減らす、という処理です。

ここまでは1レーンだけを16bitの列として表現しました。

実際のVerilogでは、16本のpair_eventはそれぞれ2300bitです。

そこで同じ重みの2300bitベクタを3本ずつCompressorへ入れ、2300レーンをまとめて3:2圧縮します。

やっていることは同じですが、Icarusから見ると2300レーンを一度にbit演算できる形になります。

これがCarry Save Adderです

この加算の仕組みが、Carry Save Adder、略してCSAです。

前に出てきたハーフアダーでは、カウンタへ+1するために1の位からcarryを上へ送っていました。

CSAではcarryをすぐ上の桁へ足し込まず、

「2倍の重みの値が1本できた」として、とりあえず保存しておく

のがポイントです。

だから Carry Save Adderです。

最後に23bitカウンタへ1回だけ足す

ここまで来れば、前回作ったbit-sliced counterへ反映できます。

前はriseが0か1だけだったので、1の位にだけ値を足してから、繰上りを23bitにわたって伝搬させていました。

つまり、16回のpair_eventに対して、

pair_event
    ↓
23bit counter update
 
 ×16回

でした。CSAを使った場合、

16 pair_event
      ↓
CSA tree
      ↓
1 / 2 / 4 / 8 / 16の加算bitset作成
      ↓
23bit counterへまとめて反映

です。

16回のカウンタ更新を行うより、

16回ぶんをCSAでまとめる処理 + 最後の1回のカウンタ更新

の方が十分軽ければ、高速化できるはずです。

これも解法を変えたわけではありません。

16回ぶんの立ち上がりを一つずつ足す代わりに、先にまとめてから一度に足す。

またIcarusへの見せ方を変えただけです。

結果発表：Carry Save Adderは効いた？

では結果です。

ここまで積み重ねてきた高速化①〜⑤に、今回のCSAを追加して比較します。

まず K=2300, N=10,000。

2-rise fusionまで  1.522秒
CSA追加後          0.627秒

約2.43倍速くなりました。

3入力を2出力へ地道に圧縮する内職仕事、ちゃんと効いてます😁

N=100,000でも、

2-rise fusionまで  14.861秒
CSA追加後           5.347秒

まで短縮しました。

こちらは約2.78倍です。

16回の23bitカウンタ更新を、11回程度の3:2圧縮と最後の1回の更新へ置き換える作戦は、Icarusではかなり有効だったようです。

ここまで速くなったので、N=1,000,000も実測してみました。

N=1,000,000  56.096秒

無事に完走です。

では最大ケース N=7,000,000 も行ってみましょう。

結果は、

300秒 TIMEOUT

でした。

5秒どころか、まだ5分では終わりません😂

最後の十分大きい完走点 N=1,000,000 から単純に線形外挿すると、

56.095892秒 × 7
= 392.671244秒
≈ 6分33秒

です。

300秒を超える予測なので、今回のTIMEOUTとも整合します。

ここまでの最大ケース予測を並べると、

最初の2300レーン版  約8時間50分
高速化①＋②          約3時間23分
高速化③             約47分56秒
高速化④             約33分44秒
高速化⑤             約17分20秒
高速化⑥              約6分33秒

まで来ました。

約8時間50分
    ↓
約6分33秒

Icarus君、ものすごく頑張った。

でもAtCoderの制限時間には、まだ遠い。

Fly! Icarus, Fly High!

……今のところ太陽に近づきすぎる心配はいらないぞ。

もう自分では思いつかない

約8時間50分から約6分33秒まで来ました。

でも、AtCoderの5秒にはまだ遠い。

そして、自分の頭では次の加速装置が思いつきません。

できそうなのは、せいぜい同時処理数のチューニングぐらいです。

ゴメン、Icarus。もうつかれたよ。

ちょっと立ち止まる

さて、ここまで来てふと我に返ります。

初めは「仮想理想FPGAならどう実装するか」を考えていたはずです。

なのに、いつの間にか Icarus Verilogが処理しやすいように回路構成を変えています。

シミュレータの都合で、FPGA側の回路構成を改造しているわけです。何をやっているんでしょうね😁

しかし、ここまで来てIcarusを見放すのもどうかなと思います。

そこで、最後の一押し、まだ改善の方法があるかどうか、ちょっと公式解説を見てみましょう。

公式解説には別の道がある

公式解説ページを見ると、いくつか高速な解法が紹介されています。

その中でまず気になったのが公式解説で紹介されているNTTです。

今までの愚直実装改善法、もうちょっとカッコよく言えばdirect法は、kごとにAを最初から最後まで見て、0→1の立ち上がりを数えていました。

最大ケースなら、700万個のAを最大2300個のkについて調べます。

これに対してNTT解法では、まずAを一度だけ先頭から走査して、

a → b

という隣接する値の組が、それぞれ何回現れたかをhistogramへ集約します。

例えば、

4 → 2 がx回出現
2 → 4 がy回出現
4 → 1 がz回出現
...

という情報を先に作ってしまうわけです。

立ち上がり判定に必要なのは、Aの何番目に現れたかではなく、

ひとつ前の値が a
今回の値が b

という組み合わせです。

ですからhistogramを作った後は、700万個のAをkごとに何度も見直す必要がありません。

あとは各shift k について、

V[(a+k) mod K] = 0
V[(b+k) mod K] = 1

になる遷移回数をhistogramから求めれば、0→1の回数を計算できます。

なお、先頭要素に対応する分は最後に別途加算します。

ここでkを変える操作は、Vを循環的にずらして積和を取る形になります。

この大量の積和計算をconvolutionへ変換し、NTTを使ってまとめて計算するのが今回試す解法です。

こちらは、

同じ仕事を速くするのではなく、仕事量そのものを減らす

方向ですね。

なお、NTT/FFTそのものについては、筆者も人様に数学的な解説ができるほど理解がこなれていません😁

ここでは仕組みの詳細には深入りしません。

NTT/FFTに興味がある方は、優れた解説記事がたくさん公開されていますのでそちらをご参照ください。

さて、これならIcarus君も迷宮を脱出できるかもしれません。

やってみます。

IcarusでNTTを試してみる

公式解説の方法1を参考に、full実装へ進む前に、MOD=998244353、NTT_N=8192 のNTT core、histogram、row convolutionなど、速度予測に必要な処理をSystemVerilogで段階的に実装しました。

NTTの往復やconvolution、K=2300でのrow convolutionなど、今回の速度予測に必要な部分の正しさは確認できました。

では速度測定です。

このNTT版はかなり時間がかかりそうだったため、処理をいくつかに分けて実測し、その結果から最大ケースの時間を見積もりました。

まず、最大ケース N=7,000,000, K=2300 のAを実際に生成し、隣接pairのhistogramを作るところまでは実測しています。

A生成 + histogram
約64.5秒

次に、K=2300で必要になる8192-point NTTを含む1行ぶんのconvolutionを測定しました。

1行、10行、100行と処理行数を増やして測ったところ、処理時間はほぼ行数に比例し、

1 rowあたり 約0.704秒

となりました。

K=2300では、このrow convolutionを2段、合計4600行ぶん実行します。

さらに、固定側のNTTと途中のtranspose処理を加えると、

A生成 + histogram        約  64.5秒
4600 row convolution    約3239.5秒
固定側NTTなど            約   0.7秒
transpose               約   4.6秒
--------------------------------
合計                     約3309秒
                        ≈ 55分

という予測になりました。

予測だけで300秒を大幅に超え、direct最終版の約6分33秒よりもかなり遅いことが分かったので、最大ケースのfull runは行いませんでした。

……。

Icarus君、墜落しました。

せっかく数学で仕事量を減らしたのに、Icarusではむしろ遅くなっています。

NTTの実装が悪いの？ Icarusが苦手なの？

ここで少し気になることがあります。

Icarusはmod乗算や大量のforループなどを回すのは苦手？

それとも

NTTの実装方法を何か間違えてる？

切り分けてみましょう。

SystemVerilog版で使ったNTTの処理を、なるべく同じ形のままnative C++へ移植します。

数学もデータ構造も同じ。

変えるのは、

Icarus vvp
↓
MSVC C++ /O2（速度最適化あり）

だけです。

同じNTTをC++で動かしてみた

最大ケースの結果は、

プロセス全体  約1.854秒
C++処理本体   約1.822秒

でした。

悪くないです。ACがとれそうな処理速度です。

Icarus版の約55分から一気に2秒弱。

どうやらNTTの良さを殺すような実装になっていたわけではなく、

Icarusが、このような計算の形をあまり得意としていなかった。

同じ算法でも、実行する計算機が変わると景色がまるで違います。

ところで、もう一つ気になる解法がある

公式解説ページには、HBitさんのAVX-512を使ったSIMD解法解説も載っています。

ざっと眺めると、

複数のkをbitでまとめる
↓
隣接するriseをまとめる
↓
CSAで数える

……。

これ、ずいぶん見覚えがあります。

我々がIcarusを速くするためにやってきたことと、かなり同じ方向です。

そして提出版のAC時間は30ms台。

無茶苦茶速い。

Icarus vvpではAVX命令は使えません。

でも、もうすでにC++実装に手を出しちゃっています。

手元マシンのCPUはi7-12700Hなので、AVX-512は使えません。

しかしAVX2（1レジスタ = 256bit）なら使えます。

だったら次は、

Icarus最終版をC++へ移植して、AVX2で同じ並列処理をやったらどうなる？

を試してみましょう。

この探検が終わるまで、Icarus君には迷宮の中でちょっと待機してもらいます。

次回

次回は、Icarusで約400秒まで鍛えたdirect解法を、そのままCPU SIMDの世界へ持っていきます。

まずはVerilog版の移植コピー。

そこから公式ページのSIMD解法を参考に、もう少し加速させてみます。

前回：
仮想理想FPGAでAtCoderを解く(2) - ABC471G② - Fly! Icarus, Fly High!

次回：
仮想理想FPGAでAtCoderを解く(4) - ABC471G④ - direct解法をAVX2で実装してみる（前編）

コード全文と実装資料：(コード実装やドキュメント作成は生成AIの助けを借りています)
第3回コード全文と実装資料