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

• 問題: No.3589 Make Ends Meet (Hard)
• ユーザー: marc2825
• 提出日時: 2026-05-29 17:33:10
• 言語: C++17 (gcc 15.2.0 + boost 1.90.0)
• 結果: AC
• 実行時間: 221 ms / 2,000 ms
• コード長: 6,366 bytes
• コンパイル時間: 1,674 ms
• コンパイル使用メモリ: 239,736 KB
• 実行使用メモリ: 9,088 KB
• 最終ジャッジ日時: 2026-07-10 20:56:01
• 合計ジャッジ時間: 4,759 ms
• ジャッジサーバーID （参考情報）: judge3_0 / judge2_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 long long BRUTE_LIMIT = 2000000;

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;
}

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;
    }
};

long long comb_limited(int n, int r, long long limit) {
    if (r < 0 || r > n) return 0;
    r = min(r, n - r);

    long double val = 1.0;
    for (int i = 1; i <= r; i++) {
        val = val * (n - r + i) / i;
        if (val > limit) return limit + 1;
    }

    long long ans = 1;
    for (int i = 1; i <= r; i++) {
        ans = ans * (n - r + i) / i;
        if (ans > limit) return limit + 1;
    }

    return ans;
}

int brute_solve(int N, int M, int K) {
    vector<pair<int,int>> edges;
    for (int u = 0; u < N; u++) {
        for (int v = u + 1; v < N; v++) {
            edges.push_back({u, v});
        }
    }

    int E = (int)edges.size();

    vector<int> removed(E, 0);
    for (int i = 0; i < M; i++) removed[i] = 1;
    sort(removed.begin(), removed.end());

    int ans = 0;

    do {
        vector<vector<int>> g(N);

        for (int i = 0; i < E; i++) {
            if (removed[i]) continue;

            auto [u, v] = edges[i];
            g[u].push_back(v);
            g[v].push_back(u);
        }

        vector<int> dist(N, -1);
        queue<int> q;

        dist[0] = 0;
        q.push(0);

        while (!q.empty()) {
            int v = q.front();
            q.pop();

            for (int to : g[v]) {
                if (dist[to] == -1) {
                    dist[to] = dist[v] + 1;
                    q.push(to);
                }
            }
        }

        if (dist[N - 1] == K) {
            ans++;
            if (ans >= MOD) ans -= MOD;
        }

    } while (next_permutation(removed.begin(), removed.end()));

    return ans;
}

// p(q) * (q^s - 1)
// 多項式は常に長さ E+1 固定
vector<int> mul_qs_minus_1_full(const vector<int> &p, int s, int E) {
    vector<int> res(E + 1, 0);

    // p(q) * q^s
    if (s <= E) {
        for (int i = 0; i + s <= E; i++) {
            res[i + s] = p[i];
        }
    }

    // - p(q)
    for (int i = 0; i <= E; i++) {
        res[i] -= p[i];
        if (res[i] < 0) res[i] += MOD;
    }

    return res;
}

// dst += scale * q^shift * src
// 多項式は常に長さ E+1 固定
void add_shift_scaled_full(vector<int> &dst, const vector<int> &src, int shift, int scale, int E) {
    if (scale == 0 || shift > E) return;

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

int fast_solve(int N, int M, int K) {
    int E = N * (N - 1) / 2;
    int R = E - M;
    int S = N - 2;

    Comb comb;
    comb.init(E);

    // dp[r][s] is a polynomial in q of fixed length E+1.
    // r: number of ordinary vertices not reached yet
    // s: size of current BFS frontier
    vector<vector<vector<int>>> dp(
        S + 1,
        vector<vector<int>>(S + 2, vector<int>(E + 1, 0))
    );

    vector<vector<char>> active(S + 1, vector<char>(S + 2, 0));

    dp[S][1][0] = 1;
    active[S][1] = 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, vector<int>(E + 1, 0))
        );

        vector<vector<char>> nactive(S + 1, vector<char>(S + 2, 0));

        for (int r = 0; r <= S; r++) {
            for (int s = 0; s <= S + 1; s++) {
                if (!active[r][s]) continue;

                vector<int> cur = dp[r][s]; // cur * (q^s - 1)^t

                for (int t = 0; t <= r; t++) {
                    int shift = t * (t - 1) / 2;
                    int ways = comb.C(r, t);

                    add_shift_scaled_full(ndp[r - t][t], cur, shift, ways, E);
                    nactive[r - t][t] = 1;

                    if (t != r) {
                        cur = mul_qs_minus_1_full(cur, s, E);
                    }
                }
            }
        }

        dp.swap(ndp);
        active.swap(nactive);
    }

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

    // N must have at least one edge to the current frontier: (q^s - 1).
    // All edges among the remaining r ordinary vertices, between them and N,
    // and between them and the current frontier are free.
    for (int r = 0; r <= S; r++) {
        for (int s = 0; s <= S + 1; s++) {
            if (!active[r][s]) continue;

            const vector<int> &p = dp[r][s];

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

            // p(q) * (q^s - 1) * q^base
            for (int i = 0; i <= E; i++) {
                int val = p[i];
                if (val == 0) continue;

                // + p(q) * q^{base+s}
                if (i + base + s <= E) {
                    H[i + base + s] += val;
                    if (H[i + base + s] >= MOD) H[i + base + s] -= MOD;
                }

                // - p(q) * q^base
                if (i + base <= E) {
                    H[i + base] -= val;
                    if (H[i + base] < 0) H[i + base] += MOD;
                }
            }
        }
    }

    ll ans = 0;

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

    return (int)ans;
}

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

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

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

    long long brute_count = comb_limited(E, M, BRUTE_LIMIT);

    if (brute_count <= BRUTE_LIMIT) {
        cout << brute_solve(N, M, K) << '\n';
    } else {
        cout << fast_solve(N, M, K) << '\n';
    }

    return 0;
}