以下の内容はhttps://nononmath.hatenablog.com/entry/2024/06/11/193707より取得しました。


上限付き単調増加列の数え上げを実装する

概要

これを実装します.なお,実装の都合上 $0 \leq B_{i} \lt A_{i}$ を満たすものを数え上げています. noshi91.hatenablog.com

実装

早速書いていきますが,思っているより長くなってしまったので部分問題に分けて書いていきます.

基本的な実装方針

まずはどのように情報を持つかの話です.元記事ではグリッドの線上を移動することを考えていましたが,それでは良い感じに書けなかったのでグリッドのマス部分を移動することを考えます. しかし,これでもうまくいきません(うまくいくかもしれませんが私は出来ませんでした).そこで,マスの辺に情報を乗せます.イメージはこんな感じです.

図 1

枠の外からの寄与は基本的に $ 0 $ ですが,左下のマスからだけは異なることに注意して下さい.

再帰部分を書く

次のような関数を書きます.

  • 下の辺の情報を受け取り,右の辺の情報を返す

まずは base-case として,幅が $ 1 $ のものを考えます.このときは,右側のどの辺にも,下から渡された値を入れてあげれば良いです. ここで注意ですが,基本的に左側からの寄与は $0$ になりますが,一番左にいるとき,すなわち l==0 のときだけ $ 1 $ が与えられます.
次に左下部分を計算します.下の辺の情報はすでに持っているのでそれを渡せば良いです.次は長方形領域の計算です.今は一旦ナイーブに寄与を二項係数 binom で与えてあげます.
これにより,上部分の情報が得られるのでその値を右上に渡してあげることで右上も計算できます.最後に右の辺の情報をマージすれば完成です.
多分以下のコードを読む方が理解しやすいと思います.今の位置が必要になるので,左端,右端,下端として l, r, d を与えています. また,各 edge の添え字は下,もしくは左が $ 0 $ となっています.

auto rec = [&](const auto &rec, int l, int r, int d, const vector<mint> &bottom_edge) -> vector<mint> {
    if (l + 1 == r) {
        // base-case
        return vector<mint>(A[l] - d, l == 0 ? mint::raw(1) : bottom_edge[0]);
    }
    int m = (l + r) / 2;
    int h = A[m] - d, w = r - m;
    // 左下の計算
    auto left_edge = rec(rec, l, m, d, vector<mint>(bottom_edge.begin(), bottom_edge.begin() + m - l));
    left_edge.resize(h);
    vector<mint> top_edge(w);
    // 左から上への寄与
    {
        for (int i = 0; i < w; i++) {
            for (int j = 0; j < h; j++) {
                top_edge[i] += left_edge[j] * binom(h - 1 - j + i, i);
            }
        }
    }
    // 下から上への寄与
    {
        for (int i = 0; i < w; i++) {
            for (int j = 0; j < w; j++) {
                top_edge[i] += bottom_edge[j + m - l] * binom(h - 1 + i - j, i - j);
            }
        }
    }
    vector<mint> right_edge(A[r - 1] - d);
    // 左から右への寄与
    {
        for (int i = 0; i < h; i++) {
            for (int j = 0; j < h; j++) {
                right_edge[i] += left_edge[j] * binom(w - 1 + i - j, i - j);
            }
        }
    }
    // 下から右への寄与
    {
        for (int i = 0; i < h; i++) {
            for (int j = 0; j < w; j++) {
                right_edge[i] += bottom_edge[j + m - l] * binom(w - 1 - j + i, i);
            }
        }
    }
    // 右上の計算
    auto upper_right = rec(rec, m, r, A[m], top_edge);
    int k = upper_right.size();
    // 右側の情報のマージ
    for (int i = 0; i < k; i++) {
        right_edge[i + h] += upper_right[i];
    }
    return right_edge;
};

畳み込み部分を書く

本質的に変わらないので「左から上への寄与」だけを計算します.

左側の辺の情報を,上から順に $L_{0},L_{1},\dots,L_{H - 1}$ ,上側の辺の情報を左から順に $T_{0},T_{1},\dots,T_{W - 1}$ とします(上の実装例とは添え字の向きが異なることに注意して下さい).このとき $$\begin{align} T_{i} &= \displaystyle\sum_{j = 0}^{H - 1}\displaystyle\binom{i + j}{i}L_{j} \\ & = \displaystyle\sum_{j = 0}^{H - 1}\dfrac{(i + j)!}{i!j!}L_{j} \\ & = \dfrac{1}{i!}\displaystyle\sum_{j = 0}^{H - 1} (i + j) ! \dfrac{L_{j}}{j!} \\ \end{align}$$ と計算すればこれは畳み込める形をしています.添え字に気をつけながら全体を実装すると次のようになります.

template <typename mint>
mint bounded_increasing_sequence(const vector<int> &A) {
    const int n = A.size();
    const int m = A[n - 1];
    // 階乗とその逆元の前計算
    vector<mint> fac(n + m + 1), finv(n + m + 1);
    {
        fac[0] = 1;
        for (int i = 1; i <= n + m; i++) fac[i] = i * fac[i - 1];
        finv[n + m] = fac[n + m].inv();
        for (int i = n + m; i >= 1; i--) finv[i - 1] = i * finv[i];
    }
    auto rec = [&](const auto &rec, int l, int r, int d, const vector<mint> &bottom_edge) -> vector<mint> {
        if (l + 1 == r) {
            return vector<mint>(A[l] - d, l == 0 ? mint::raw(1) : bottom_edge[0]);
        }
        int m = (l + r) / 2;
        int h = A[m] - d, w = r - m;
        // 左下の計算
        auto left_edge = rec(rec, l, m, d, vector<mint>(bottom_edge.begin(), bottom_edge.begin() + m - l));
        left_edge.resize(h);
        vector<mint> top_edge(w);
        // 左から上への寄与
        if (h) {
            vector<mint> f(h), g(h + w);
            for (int i = 0; i < h; i++) f[i] = left_edge[i] * finv[h - 1 - i];
            for (int i = 0; i < h + w; i++) g[i] = fac[i];
            f = convolution(f, g);
            for (int i = 0; i < w; i++) top_edge[i] += finv[i] * f[h - 1 + i];
        }
        // 下から上への寄与
        if (h) {
            vector<mint> f(w), g(w);
            for (int i = 0; i < w; i++) f[i] = bottom_edge[i + m - l];
            for (int i = 0; i < w; i++) g[i] = fac[h - 1 + i] * finv[i];
            f = convolution(f, g);
            for (int i = 0; i < w; i++) top_edge[i] += finv[h - 1] * f[i];
        } else {
            for (int i = 0; i < w; i++) top_edge[i] = bottom_edge[i + m - l];
        }
        vector<mint> right_edge(A[r - 1] - d);
        // 左から右への寄与
        if (h) {
            vector<mint> f(h), g(h + w);
            for (int i = 0; i < h; i++) f[i] = left_edge[i];
            for (int i = 0; i < h + w; i++) g[i] = fac[w - 1 + i] * finv[i];
            f = convolution(f, g);
            for (int i = 0; i < h; i++) right_edge[i] += finv[w - 1] * f[i];
        }
        // 下から右への寄与
        if (h) {
            vector<mint> f(w), g(h + w);
            for (int i = 0; i < w; i++) f[i] = bottom_edge[m - l + i] * finv[w - 1 - i];
            for (int i = 0; i < h + w; i++) g[i] = fac[i];
            f = convolution(f, g);
            for (int i = 0; i < h; i++) right_edge[i] += finv[i] * f[w - 1 + i];
        }
        vector<mint> upper_right = rec(rec, m, r, A[m], top_edge);
        int k = upper_right.size();
        // 右側の情報のマージ
        for (int i = 0; i < k; i++) right_edge[i + h] += upper_right[i];
        return right_edge;
    };
    vector<mint> right_edge = rec(rec, 0, n, 0, vector<mint>(n));
    mint res = 0;
    for (mint x : right_edge) res += x;
    return res;
}

終わりに

verify はどこですればいいんでしょうか(一応手元で愚直と合わせましたが).ちなみに convolutionmodintACL を用いた場合,手元での実行時間は $N=M=10^{5}$ のケースで 360 ~ 380 ms 程度でした.




以上の内容はhttps://nononmath.hatenablog.com/entry/2024/06/11/193707より取得しました。
このページはhttp://font.textar.tv/のウェブフォントを使用してます

不具合報告/要望等はこちらへお願いします。
モバイルやる夫Viewer Ver0.14