#1167154 (C++17) No.3589 Make Ends Meet (Hard)

• 問題: No.3589 Make Ends Meet (Hard)
• ユーザー: marc2825
• 提出日時: 2026-05-29 16:57:33
• 言語: C++17 (gcc 15.2.0 + boost 1.90.0)
• 結果: AC
• 実行時間: 756 ms / 2,000 ms
• コード長: 6,537 bytes
• コンパイル時間: 1,511 ms
• コンパイル使用メモリ: 229,464 KB
• 実行使用メモリ: 9,216 KB
• 最終ジャッジ日時: 2026-07-10 20:55:21
• 合計ジャッジ時間: 8,039 ms
• ジャッジサーバーID （参考情報）: judge3_0 / judge1_1

ファイルパターン	結果
sample	AC * 3
other	AC * 47

ソースコード

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

static const int MOD = 998244353;
static const int G = 3;

ll modpow(ll a, ll e) {
    ll r = 1;
    while (e > 0) {
        if (e & 1) r = r * a % MOD;
        a = a * a % MOD;
        e >>= 1;
    }
    return r;
}

void ntt(vector<int> &a, bool inv) {
    int n = (int)a.size();

    for (int i = 1, j = 0; i < n; i++) {
        int bit = n >> 1;
        while (j & bit) {
            j ^= bit;
            bit >>= 1;
        }
        j ^= bit;
        if (i < j) swap(a[i], a[j]);
    }

    for (int len = 2; len <= n; len <<= 1) {
        int wlen = (int)modpow(G, (MOD - 1) / len);
        if (inv) wlen = (int)modpow(wlen, MOD - 2);

        for (int i = 0; i < n; i += len) {
            ll w = 1;
            int half = len >> 1;

            for (int j = 0; j < half; j++) {
                int u = a[i + j];
                int v = (int)(w * a[i + j + half] % MOD);

                int x = u + v;
                if (x >= MOD) x -= MOD;
                a[i + j] = x;

                x = u - v;
                if (x < 0) x += MOD;
                a[i + j + half] = x;

                w = w * wlen % MOD;
            }
        }
    }

    if (inv) {
        int inv_n = (int)modpow(n, MOD - 2);
        for (int &x : a) x = (ll)x * inv_n % MOD;
    }
}

struct Comb {
    vector<int> fact, ifact;

    void init(int n) {
        fact.assign(n + 1, 1);
        ifact.assign(n + 1, 1);
        for (int i = 1; i <= n; i++) fact[i] = (ll)fact[i - 1] * i % MOD;
        ifact[n] = (int)modpow(fact[n], MOD - 2);
        for (int i = n; i >= 1; i--) ifact[i - 1] = (ll)ifact[i] * i % MOD;
    }

    int C(int n, int r) const {
        if (r < 0 || r > n) return 0;
        return (ll)fact[n] * ifact[r] % MOD * ifact[n - r] % MOD;
    }
};

static inline void trim(vector<int> &a) {
    while (!a.empty() && a.back() == 0) a.pop_back();
}

vector<int> convolution_self(vector<int> a, vector<int> b) {
    if (a.empty() || b.empty()) return {};

    int need = (int)a.size() + (int)b.size() - 1;

    // 小さいものは素直な畳み込みの方が速い
    if ((ll)a.size() * b.size() <= 4096) {
        vector<int> c(need, 0);
        for (int i = 0; i < (int)a.size(); i++) {
            if (a[i] == 0) continue;
            for (int j = 0; j < (int)b.size(); j++) {
                if (b[j] == 0) continue;
                c[i + j] = (c[i + j] + (ll)a[i] * b[j]) % MOD;
            }
        }
        trim(c);
        return c;
    }

    int n = 1;
    while (n < need) n <<= 1;

    a.resize(n);
    b.resize(n);

    ntt(a, false);
    ntt(b, false);

    for (int i = 0; i < n; i++) {
        a[i] = (ll)a[i] * b[i] % MOD;
    }

    ntt(a, true);

    a.resize(need);
    trim(a);
    return a;
}

vector<int> multiply_poly(const vector<int> &a, const vector<int> &b, int E) {
    if (a.empty() || b.empty()) return {};

    vector<int> c = convolution_self(a, b);

    if ((int)c.size() > E + 1) {
        c.resize(E + 1);
    }

    trim(c);
    return c;
}

vector<int> shift_poly(const vector<int> &p, int shift, int E) {
    if (p.empty() || shift > E) return {};

    int need = min(E + 1, shift + (int)p.size());
    vector<int> res(need, 0);

    for (int i = 0; i < (int)p.size() && i + shift <= E; i++) {
        res[i + shift] = p[i];
    }

    trim(res);
    return res;
}

void add_scaled(vector<int> &dst, const vector<int> &src, int scale, int E) {
    if (src.empty() || scale == 0) return;

    if ((int)dst.size() < (int)src.size()) {
        dst.resize(src.size(), 0);
    }

    for (int i = 0; i < (int)src.size(); i++) {
        if (src[i] == 0) continue;
        dst[i] = (dst[i] + (ll)src[i] * scale) % MOD;
    }

    trim(dst);
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    int N, M, K;
    cin >> N >> M >> K;

    int E = N * (N - 1) / 2;
    int R = E - M;
    int S = N - 2;

    Comb comb;
    comb.init(E);

    // base_poly[s] = q^s - 1
    vector<vector<int>> base_poly(S + 2);

    for (int s = 0; s <= S + 1; s++) {
        vector<int> p(s + 1, 0);
        p[0] = MOD - 1;
        p[s] += 1;
        if (p[s] >= MOD) p[s] -= MOD;
        trim(p);
        base_poly[s] = p;
    }

    // pow_poly[s][t] = (q^s - 1)^t
    vector<vector<vector<int>>> pow_poly(S + 2, vector<vector<int>>(S + 1));

    for (int s = 0; s <= S + 1; s++) {
        pow_poly[s][0] = vector<int>{1};
        for (int t = 1; t <= S; t++) {
            pow_poly[s][t] = multiply_poly(pow_poly[s][t - 1], base_poly[s], E);
        }
    }

    // kernel[s][t] = (q^s - 1)^t * q^{t(t-1)/2}
    vector<vector<vector<int>>> kernel(S + 2, vector<vector<int>>(S + 1));

    for (int s = 0; s <= S + 1; s++) {
        for (int t = 0; t <= S; t++) {
            int shift = t * (t - 1) / 2;
            kernel[s][t] = shift_poly(pow_poly[s][t], shift, E);
        }
    }

    // dp[r][s] is a polynomial in q.
    vector<vector<vector<int>>> dp(S + 1, vector<vector<int>>(S + 2));
    dp[S][1] = vector<int>{1};

    // Build layers 1,2,...,K-1 without reaching vertex N.
    for (int step = 0; step < K - 1; step++) {
        vector<vector<vector<int>>> ndp(S + 1, vector<vector<int>>(S + 2));

        for (int r = 0; r <= S; r++) {
            for (int s = 0; s <= S + 1; s++) {
                if (dp[r][s].empty()) continue;

                for (int t = 0; t <= r; t++) {
                    int ways = comb.C(r, t);

                    vector<int> prod = multiply_poly(dp[r][s], kernel[s][t], E);
                    add_scaled(ndp[r - t][t], prod, ways, E);
                }
            }
        }

        dp.swap(ndp);
    }

    // H(q): generating polynomial for dist(1,N)=K in q=1+x.
    vector<int> H(E + 1, 0);

    for (int r = 0; r <= S; r++) {
        for (int s = 0; s <= S + 1; s++) {
            const vector<int> &p = dp[r][s];
            if (p.empty()) continue;

            int base = s * r + r + r * (r - 1) / 2;

            // final factor: (q^s - 1) q^base
            vector<int> final_factor = shift_poly(base_poly[s], base, E);
            vector<int> prod = multiply_poly(p, final_factor, E);

            add_scaled(H, prod, 1, E);
        }
    }

    // We need [x^R] H(1+x).
    // If H(q)=sum_c h_c q^c, then [x^R]H(1+x)=sum_c h_c*C(c,R).
    ll ans = 0;

    for (int c = R; c <= E; c++) {
        if (c < (int)H.size() && H[c] != 0) {
            ans += (ll)H[c] * comb.C(c, R) % MOD;
            ans %= MOD;
        }
    }

    cout << ans << '\n';
    return 0;
}