はじめに

前回、ABC471 E問題をShrike-Liteへ実装しました。

配列を保存せず、

prefix_sum
pair_sum
square_sum

の3個を更新すれば計算できます。公式解説の裏ワザを使えば pair_sum も消せますが、Baselineでは、

まあ3個でも十分少ないし、裏ワザはまだいいでしょう。

と、そのまま実装しました。

ところがPNRすると、

Type=L : 577 / 140

全然入りません。

今回は、ここから始まります。

この記事の実装とその資料は以下のリンクで公開しています。

第38回コード全文と実装資料

ここで用語整理：SynthとPNR

AtCoder寄り読者の方には、ここから何度も出てくる Synth や PNR は馴染みが薄いかもしれません。

FPGAに回路を書き込む手順は、ざっくり言えば、

Verilog
    どの信号を、どのタイミングで、どう処理するかを書く

Synth（合成）
    Verilogの記述を、FPGA上で実現できる論理回路へ変換する

PNR（Place and Route）
    Synthでできた回路を、実際のFPGA内のどこへ配置し、
    どう配線するかを決める

という流れです。

この記事では、

Post-SynthesisのLUT数
    -> 合成した論理回路がどのくらい大きいか

PNRのType=L
    -> Shrike-Liteの実機に配置するときの必要面積はどれぐらいか

を見ることが多くなります。

Synthでは入りそうに見えたのに、PNRすると入らない、というのも今回の重要なポイントです。

まず普通に面積を削ってみる

まずは3レジスタ構成のまま、細かい見直しを行いました。

演算器や一時レジスタの共有
operand stagingの整理
状態やMUXの整理

といった変更です。

一見小さくなりそうな変更でも、MUXやdecodeが増えて逆にLUTが増えることもありました。

最終的にこの3レジスタ版では、

Type=L
Baseline      : 577 / 140
見直し後      : 571 / 140

となりました。

6 CLB減りました。

減ってはいますが、全然足りません。

そこで、ここで初めて前編で紹介した公式解説の裏ワザを使ってみます。

pair_sumをやめて2個だけにする

前編で紹介した、

2 * pair_sum
= prefix_sum^2 - square_sum

を使えば、pair_sum を保持する必要はありません。

ここで3レジスタから、

prefix_sum
square_sum

の2レジスタ構成へ変更します。

「レジスタを1個丸ごと消せるなら、今度こそ効くのでは？」

Post-Synthesis結果は、

LUT    : 2121
FF     : 948
CARRY4 : 121

でした。

さらにPNRすると、

Type=L : 561 / 140

です。

3レジスタ見直し版の571から、さらに10 CLB減りました。

Baselineから見ると、

577 -> 571 -> 561

です。

16 CLB減りました。

でも必要なのは140以下です。

全然景色が変わりません。

何が面積を食っているのか

そこで30bit・2個管理版を基準に、一部の計算結果を固定値へ置き換えて合成しました。

                         LUT   FF   CARRY4   LUT差
通常                    2121  948   121       -
powの結果を固定          1845  881   119    -276
係数とpowの結果を固定    1223  625    73    -898

固定すると周辺の制御やMUXも一緒に消え、ABC9の再マッピングも入るため、差分をそのまま各回路の面積とは読めません。

それでも、組合せ係数・modular pow・その周辺がかなり大きいことは分かります。

モジュラ問題、割り算が怖い

今回、C(n,k) を求めるためにmodular inverseが必要です。

Pythonなら、

inv = pow(x, MOD-2, MOD)

で終わります。

しかしFPGAでは、その pow が30bitの回路になります。

主要な集計値は2～3個しかないのに、実際には、

30bit modular multiplier
modular add/sub
combination計算
Fermat inverse
pow sequencer
operand MUX
state/control
SPI protocol

まで住んでいました。

住人2,3人のつもりで間取りを考えていたら、管理会社と工務店まで同居していたようなものです。

30bitはさすがに厳しい

そこで演算器ではなく、問題制約そのものを縮小してみます。

元のABC471Eは、

K <= N <= 200000
Ai <= 10^9
MOD = 998244353

です。

これを8bit幅で扱える、

K <= N <= 250
Ai <= 250
MOD = 251

へ縮めました。

言ってしまえば、

ABC471Eの1/4縮尺模型

です。

これならさすがに入るでしょう。

8bitでも入りません

普通に8bit化しただけの実装は、

803 LUT

程度でした。まだ大きい。

ベッドや冷蔵庫は小さくなりましたが、壁や廊下などの固定費が思った以上に残っています。

Codex、限界まで狭小住宅に詰め込んでみてくれ

ここからは、

RTLを少し直す
Icarus
合成
resourceを見る
また直す

の繰り返しです。

正直、もう飽きました。

合成までCodexにやってもらえばいいんじゃない？

ForgeFPGA Workshopは合成にYosysを使っています。Codexに調べてもらうと、ForgeFPGA付属YosysとWorkshop生成のsynth scriptを使い、CLIからPost-Synthesisを再現できるようです。

そこで、

Codex君、合成結果を見ながら削ってみて。PNRが見えてきたら呼んで

という方式へ変更しました。

この時点ではPNRはまだ人力です。

557 LUTまで来た

試行を繰り返した結果、W8版は、

557 LUT
163 FF
31 CARRY4

まで減りました。

Shrike-LiteのType=L容量は140なので、

ceil(557 / 4) = 140

理論上はちょうど入る可能性があります。

今度こそ。

PNRしてみる

ForgeFPGAの通常PNRフローを人力で回します。

Type=L : 165 / 140

入りません。

557 LUTから計算した140 CLBは、あくまで理論的な下限でした。実際のmapping・packing・配置結果では、Type=Lは165となりました。

じゃあ、どこまで問題を小さくすれば入るのか

ある規模の問題はShrike-Liteに実装できるか、を試してきましたが、先行きが見えません。

そこで方針を変えて「どのサイズまで縮めれば実装できるのか」を探索することにします。

各WIDTHでMODを固定し、N_MAX だけを変えてみます。

WIDTH=5  MOD=31
WIDTH=6  MOD=61
WIDTH=7  MOD=127
WIDTH=8  MOD=251
WIDTH=9  MOD=509

有効範囲は、

1 <= N <= N_MAX
1 <= K <= N
0 <= Ai < MOD

です。

各WIDTHで N_MAX=MOD-1 から始め、入らなければ段階的に下げます。

PNRまでCodex君にやってもらえないか

しかし候補ごとに人力PNRするのも面倒です。

Renesas社の公開マニュアルを見ても、PNRをCLIから実行する方法は見つけられませんでした。

Codex君、これもそっちで何とかならない？

乱暴な依頼にもかかわらず、CodexはForgeFPGA Workshopの動きを調べ、専用のコピーprojectから通常のGUI PNR flowを実行し、ログを回収する方法を見つけてきました。

W8の既知結果、

Post-Synthesis : 557 LUT / 163 FF / 31 CARRY4
PNR             : Type=L 165 / 140

も再現。

すごいなCodex君。

自動境界テスタ完成

これで、

候補生成
Icarus
Yosys synthesis
resource取得
必要ならForgeFPGA PNR
結果回収

まで完全自動化できました。

人間はボンヤリお茶を飲みながら待っているだけです。

よし。

RUN!

なお、このPNR自動化はRenesasの公開マニュアルに記載された操作方法ではなく、Workshopの通常flowをこちらで調べて利用した実験用の方法です。メーカー保証の方法ではないため、試す場合は自己責任でお願いします。

そして結果はガタガタだった

結果が出ました。代表的な観測点です。

W5 N=30  : 149 Type=L  FAIL
W5 N=29  : 109 Type=L  FIT / 54.174 MHz
W5 N=15  : 121 Type=L  FIT / 51.634 MHz

W6 N=60  : 138 Type=L  FIT / 48.757 MHz

W7 N=126 : 152 Type=L  FAIL
W7 N=63  : 139 Type=L  FIT / 33.118 MHz

W8 N=250 : 165 Type=L  FAIL
W8 N=63  : 141 Type=L  FAIL
W8 N=31  : 151 Type=L  FAIL
W8 N=15  : 141 Type=L  FAIL
W8 N=7   : 122 Type=L  FIT / 45.781 MHz

W9 N=508 : Post-Synthesis段階でSCREEN FAIL

何だこれは。

Nを減らしても素直に小さくならない

W8では、

N=63 -> 141 Type=L
N=31 -> 151 Type=L
N=15 -> 141 Type=L

N_MAXを減らしているのに、減る、増える、また減る。

W5では、

N=29 -> 109 Type=L
N=30 -> 149 Type=L

N_MAXが1違うだけで40 CLB違います。COUNT_WIDTHはどちらも5bitです。

Yosys/ABC9 mapping、carry-chain packing、MUX、placer geometryなどが絡んだ結果と考えられますが、今回の観測だけで原因を1つに決めることはできません。

logical LUTだけ見てもわからない

さらに、

W5 N=30 : 452 LUT -> 149 Type=L FAIL
W6 N=60 : 508 LUT -> 138 Type=L FIT

logical LUTが少ない方が、物理CLBは多い。

W8 N=250も557 LUTなのに165 Type=Lでした。

少なくとも今回のForge FPGAでは、

LUT数を4で割れば物理CLB数がだいたい分かる

とは考えない方がよさそうです。

W8で確認できたfit点はN=7

W8で確認できたPNR_COMPLETE点は、

N_MAX=7
Type=L=122/140
post-route=45.781 MHz

でした。

N=8～14を調べていないため、「W8の最大Nは7」とは言えません。

ただ、

8bit世界でNが一桁まで来たら、もう別の解法を考えた方がよくないか？

となり、探索はここで終了させることにました。

W=8, N=7なら、250以下の数字が7個、そこからk個選ぶ組み合わせの和の二乗の和、です。

この規模になると、そろばん名人の愚直解法が競合相手です。

k<=2なら負けそうです。

ここまで縮めるなら別解法でよい

N_MAX=7なら C(n,k) も小さな表で済み、Fermat逆元や汎用pow回路を残す意味も薄くなります。

問題を縮めればいつかは入るでしょう。

でも、縮めすぎれば元の一般解法そのものがオーバースペックです。

縮尺模型にも限度があります。

今回は失敗です

ABC471Eの原問題をShrike-Liteへ実装する、という目的には失敗しました。

入りませんでした。

今回の教訓を一言で書くなら、

モジュラ問題、とくに割り算と逆元が出てきたら、そっと目をそらす

でしょうか。

それでも、いい経験にはなりました。

副産物としてForgeFPGA WorkshopのCodex自動操作までできるようになりました。

おまけ：それでもVerilogでAtCoderには出してみたい

Shrike-Liteには入りませんでした。

しかしAtCoderが要求しているのは、Icarus vvp環境で正しい答えを出してください、です。

そこで次は、回路面積もクロックもクリティカルパスも完全に無視した、競プロコードを書いてみます。

module main;
    localparam signed [63:0] MOD = 64'd998244353;
    localparam integer STDIN = 32'h8000_0000;

    integer n;
    integer k;
    integer i;
    integer scan_result;

    reg signed [63:0] x;
    reg signed [63:0] prefix_sum;
    reg signed [63:0] pair_sum;
    reg signed [63:0] square_sum;
    reg signed [63:0] choose_single;
    reg signed [63:0] choose_pair;
    reg signed [63:0] inverse_n_minus_one;
    reg signed [63:0] pair_term;
    reg signed [63:0] answer;

    // base^exponent mod MODを二分累乗法で求める。
    // 組合せ係数では、Fermatの小定理による逆元計算に使用する。
    function automatic signed [63:0] mod_pow;
        input signed [63:0] base_input;
        input signed [63:0] exponent_input;
        reg signed [63:0] base;
        reg signed [63:0] exponent;
        reg signed [63:0] result;
        begin
            base = base_input % MOD;
            exponent = exponent_input;
            result = 64'd1;

            while (exponent > 0) begin
                if (exponent[0]) begin
                    result = (result * base) % MOD;
                end
                base = (base * base) % MOD;
                exponent = exponent >> 1;
            end

            mod_pow = result;
        end
    endfunction

    // 分子の積、分母の積、分母の逆元から
    // C(input_n - 1, input_k - 1)を求める。
    function automatic signed [63:0] combination_n_minus_one;
        input integer input_n;
        input integer input_k;
        integer r;
        integer j;
        reg signed [63:0] numerator;
        reg signed [63:0] denominator;
        begin
            r = input_k - 1;
            if (r > input_n - input_k) begin
                r = input_n - input_k;
            end

            numerator = 64'd1;
            denominator = 64'd1;
            for (j = 1; j <= r; j = j + 1) begin
                numerator = (numerator * (input_n - 1 - r + j)) % MOD;
                denominator = (denominator * j) % MOD;
            end

            combination_n_minus_one =
                (numerator * mod_pow(denominator, MOD - 2)) % MOD;
        end
    endfunction

    initial begin
        scan_result = $fscanf(STDIN, "%d %d", n, k);
        if (scan_result != 2) begin
            $finish(0);
        end

        prefix_sum = 64'd0;
        pair_sum = 64'd0;
        square_sum = 64'd0;

        // 前編で導出した3つの集計値を、そのまま更新する。
        for (i = 0; i < n; i = i + 1) begin
            scan_result = $fscanf(STDIN, "%d", x);
            x = x % MOD;

            pair_sum = (pair_sum + prefix_sum * x) % MOD;
            prefix_sum = (prefix_sum + x) % MOD;
            square_sum = (square_sum + x * x) % MOD;
        end

        choose_single = combination_n_minus_one(n, k);

        // C(N-2,K-2) = C(N-1,K-1) * (K-1) / (N-1)
        // K=1では2要素の積の寄与はない。
        if (k == 1) begin
            choose_pair = 64'd0;
        end else begin
            inverse_n_minus_one = mod_pow(n - 1, MOD - 2);
            choose_pair = (choose_single * (k - 1)) % MOD;
            choose_pair = (choose_pair * inverse_n_minus_one) % MOD;
        end

        // 最終式:
        // C(N-1,K-1) * square_sum + 2 * C(N-2,K-2) * pair_sum
        pair_term = (choose_pair * pair_sum) % MOD;
        answer = ((choose_single * square_sum) % MOD
                  + (2 * pair_term) % MOD) % MOD;

        $display("%0d", answer);
        $finish(0);
    end
endmodule

この実装は前編で紹介した3レジスタ構成そのままです。

AtCoderコードテスト環境でも動作を確認していますので、まず動かしてみたい方は、そのままコピペして試せます。

もっと遊びたい方は、2レジスタ化や % MOD の削減などを試して正式提出してみてください。

こちらで作った最適化参考版もGitHubへ置いておきます。

なお、正式提出前のコードなので、AC保証はありません。その点はご了承ください。

前回：
Shrike-LiteでAtCoder問題を解く(37) ABC471E（前編） - レジスタ数個なら入ると思った

次回：
Shrike-LiteでAtCoder問題を解く(39)：Interlude - ABC472の各問題をFPGAなどの視点で見てみる

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