なんとか解けた。若干エスパー気味に解いた。

問題へのリンク

問題概要

頂点数 $N$ 、辺数 $M$ の無向単純グラフが与えられる。

各 $k = 0, 1, \dots, N$ に対して、この誘導部分グラフ (頂点集合はそのまま、辺集合は部分集合) であって、次数が奇数の頂点が $k$ 個であるようなものの個数を 999244353 で割ったあまりを求めよ。

制約

$N, M \le 5000$

考えたこと

まずは色々なグラフで試してみることにした。そうしていくうちに

連結成分ごとに独立して考えて、得られた解を多項式としての掛け算をしたものが答えになる (それはそう)
頂点数 $N$ のパスグラフに対しては $({}_{N}{\rm C}_{0}, 0, {}_{N}{\rm C}_{2}, 0, {}_{N}{\rm C}_{4}, \dots)$ が答え

という風になってることがわかってきた。さらにパスグラフでなくても、木であれば答えが一定になるっぽいことがわかって来た。そしてさらに色んなケースを試すことで、

連結であれば、辺数が 1 増えると、各 $k$ に対する答えが全体的に 2 倍ずつになる

ということが見えて来た。つまり、連結グラフであれば、頂点数と辺数のみに依存して答えが求まるのだ。これを素直に実装すると通った。

背景

グラフの接続行列を用いた ${\rm F}_{2}$ 線形代数っぽい雰囲気の問題では、「全域木に帰着して考える」という方法が有効打になることが多いみたいだ。

参考資料：「接続行列を係数にもつ線形方程式」の競プロでの応用　(opt さん)

今回の問題も、連結グラフに対して

連結グラフに辺を 1 本追加すると、全体の答えが 2 倍になる
木についての答えが $({}_{N}{\rm C}_{0}, 0, {}_{N}{\rm C}_{2}, 0, {}_{N}{\rm C}_{4}, \dots)$ となる

ということを解決していけば OK。

前者 (全域木へと帰着)

こっちは比較的わかりやすい。追加する辺を $e$ としたとき、「次数が奇数となる頂点集合」を変えないように、

辺 $e$ を含まない誘導部分グラフ
辺 $e$ を含む誘導部分グラフ

との間に一対一対応を作ることができるのだ。具体的には、 $e$ を含むサイクルをとることができる (全域木の基本サイクル) ので、そのサイクル上の各辺について「部分グラフに選ぶ」「部分グラフに選ばない」を反転していけば OK。

よって、連結グラフに辺 $e$ を追加するとき、全体の答えは 2 倍になることがわかった。

後者 (木の場合)

まず奇数次数の頂点はかならず偶数個になることに注意する (有名問題)。そして、次のことが成立する！！！

$k$ が偶数のとき、 $N$ 個の頂点から $k$ 個の頂点を選べば、それらの頂点の次数を奇数にして、それ以外の頂点の次数を偶数にするようなものは、ただ一つ存在する。

具体的にはこれは、葉から順に決まって行くイメージなのだ。各頂点の偶奇を決めてあげると、葉に接続している各辺に対して、その辺を部分グラフとして採用すべきかどうかが決まって行く。そして葉から内側への辺に対しても Greedy に採用すべきかどうかが決まって行くのだ。

よって、 $({}_{N}{\rm C}_{0}, 0, {}_{N}{\rm C}_{2}, 0, {}_{N}{\rm C}_{4}, \dots)$ が答えになる。

補足

maspy さんの話が面白かった

これの Ker は、
・辺集合であって、任意の点の次数が偶数なもの全体
・cycle の disjoint sum として書けるもの全体といってもよい

これはサイクル空間と呼ばれて、全域森をとると基底が手に入ります。（森に含まれない辺をとるごとに、森と合わせてサイクルがひとつとれて、それらが基底） https://t.co/O1LXOdtyja
— maspy (@maspy_stars) 2021年3月22日

https://twitter.com/sdfi13jfx/status/1373969475465142272?s=20

コード

各連結成分に対する解は明示的に得られるので、それらを多項式として掛け算していけば OK。

特に、「二分木のような計算順序」 (この記事を参照) を用いて FFT を活用していくことで、計算量は $O(N (\log N)^{2})$ となる。

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

// modint
template<int MOD> struct Fp {
    long long val;
    constexpr Fp(long long v = 0) noexcept : val(v % MOD) {
        if (val < 0) val += MOD;
    }
    constexpr int getmod() const { return MOD; }
    constexpr Fp operator - () const noexcept {
        return val ? MOD - val : 0;
    }
    constexpr Fp operator + (const Fp& r) const noexcept { return Fp(*this) += r; }
    constexpr Fp operator - (const Fp& r) const noexcept { return Fp(*this) -= r; }
    constexpr Fp operator * (const Fp& r) const noexcept { return Fp(*this) *= r; }
    constexpr Fp operator / (const Fp& r) const noexcept { return Fp(*this) /= r; }
    constexpr Fp& operator += (const Fp& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Fp& operator -= (const Fp& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Fp& operator *= (const Fp& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Fp& operator /= (const Fp& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator < (const Fp& r) const noexcept {
        return this->val < r.val;
    }
    constexpr bool operator == (const Fp& r) const noexcept {
        return this->val == r.val;
    }
    constexpr bool operator != (const Fp& r) const noexcept {
        return this->val != r.val;
    }
    friend constexpr istream& operator >> (istream& is, Fp<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Fp<MOD>& x) noexcept {
        return os << x.val;
    }
    friend constexpr Fp<MOD> modpow(const Fp<MOD>& r, long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return modpow(modinv(r), -n);
        auto t = modpow(r, n / 2);
        t = t * t;
        if (n & 1) t = t * r;
        return t;
    }
    friend constexpr Fp<MOD> modinv(const Fp<MOD>& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Fp<MOD>(u);
    }
};

namespace NTT {
    long long modpow(long long a, long long n, int mod) {
        long long res = 1;
        while (n > 0) {
            if (n & 1) res = res * a % mod;
            a = a * a % mod;
            n >>= 1;
        }
        return res;
    }

    long long modinv(long long a, int mod) {
        long long b = mod, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        u %= mod;
        if (u < 0) u += mod;
        return u;
    }

    int calc_primitive_root(int mod) {
        if (mod == 2) return 1;
        if (mod == 167772161) return 3;
        if (mod == 469762049) return 3;
        if (mod == 754974721) return 11;
        if (mod == 998244353) return 3;
        int divs[20] = {};
        divs[0] = 2;
        int cnt = 1;
        long long x = (mod - 1) / 2;
        while (x % 2 == 0) x /= 2;
        for (long long i = 3; i * i <= x; i += 2) {
            if (x % i == 0) {
                divs[cnt++] = i;
                while (x % i == 0) x /= i;
            }
        }
        if (x > 1) divs[cnt++] = x;
        for (int g = 2;; g++) {
            bool ok = true;
            for (int i = 0; i < cnt; i++) {
                if (modpow(g, (mod - 1) / divs[i], mod) == 1) {
                    ok = false;
                    break;
                }
            }
            if (ok) return g;
        }
    }

    int get_fft_size(int N, int M) {
        int size_a = 1, size_b = 1;
        while (size_a < N) size_a <<= 1;
        while (size_b < M) size_b <<= 1;
        return max(size_a, size_b) << 1;
    }

    // number-theoretic transform
    template<class mint> void trans(vector<mint>& v, bool inv = false) {
        if (v.empty()) return;
        int N = (int)v.size();
        int MOD = v[0].getmod();
        int PR = calc_primitive_root(MOD);
        static bool first = true;
        static vector<long long> vbw(30), vibw(30);
        if (first) {
            first = false;
            for (int k = 0; k < 30; ++k) {
                vbw[k] = modpow(PR, (MOD - 1) >> (k + 1), MOD);
                vibw[k] = modinv(vbw[k], MOD);
            }
        }
        for (int i = 0, j = 1; j < N - 1; j++) {
            for (int k = N >> 1; k > (i ^= k); k >>= 1);
            if (i > j) swap(v[i], v[j]);
        }
        for (int k = 0, t = 2; t <= N; ++k, t <<= 1) {
            long long bw = vbw[k];
            if (inv) bw = vibw[k];
            for (int i = 0; i < N; i += t) {
                mint w = 1;
                for (int j = 0; j < t/2; ++j) {
                    int j1 = i + j, j2 = i + j + t/2;
                    mint c1 = v[j1], c2 = v[j2] * w;
                    v[j1] = c1 + c2;
                    v[j2] = c1 - c2;
                    w *= bw;
                }
            }
        }
        if (inv) {
            long long invN = modinv(N, MOD);
            for (int i = 0; i < N; ++i) v[i] = v[i] * invN;
        }
    }

    // for garner
    static constexpr int MOD0 = 754974721;
    static constexpr int MOD1 = 167772161;
    static constexpr int MOD2 = 469762049;
    using mint0 = Fp<MOD0>;
    using mint1 = Fp<MOD1>;
    using mint2 = Fp<MOD2>;
    static const mint1 imod0 = 95869806; // modinv(MOD0, MOD1);
    static const mint2 imod1 = 104391568; // modinv(MOD1, MOD2);
    static const mint2 imod01 = 187290749; // imod1 / MOD0;

    // small case (T = mint, long long)
    template<class T> vector<T> naive_mul 
    (const vector<T>& A, const vector<T>& B) {
        if (A.empty() || B.empty()) return {};
        int N = (int)A.size(), M = (int)B.size();
        vector<T> res(N + M - 1);
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < M; ++j)
                res[i + j] += A[i] * B[j];
        return res;
    }

    // mint
    template<class mint> vector<mint> mul
    (const vector<mint>& A, const vector<mint>& B) {
        if (A.empty() || B.empty()) return {};
        int N = (int)A.size(), M = (int)B.size();
        if (min(N, M) < 30) return naive_mul(A, B);
        int MOD = A[0].getmod();
        int size_fft = get_fft_size(N, M);
        if (MOD == 998244353) {
            vector<mint> a(size_fft), b(size_fft), c(size_fft);
            for (int i = 0; i < N; ++i) a[i] = A[i];
            for (int i = 0; i < M; ++i) b[i] = B[i];
            trans(a), trans(b);
            vector<mint> res(size_fft);
            for (int i = 0; i < size_fft; ++i) res[i] = a[i] * b[i];
            trans(res, true);
            res.resize(N + M - 1);
            return res;
        }
        vector<mint0> a0(size_fft, 0), b0(size_fft, 0), c0(size_fft, 0);
        vector<mint1> a1(size_fft, 0), b1(size_fft, 0), c1(size_fft, 0);
        vector<mint2> a2(size_fft, 0), b2(size_fft, 0), c2(size_fft, 0);
        for (int i = 0; i < N; ++i)
            a0[i] = A[i].val, a1[i] = A[i].val, a2[i] = A[i].val;
        for (int i = 0; i < M; ++i)
            b0[i] = B[i].val, b1[i] = B[i].val, b2[i] = B[i].val;
        trans(a0), trans(a1), trans(a2), trans(b0), trans(b1), trans(b2);
        for (int i = 0; i < size_fft; ++i) {
            c0[i] = a0[i] * b0[i];
            c1[i] = a1[i] * b1[i];
            c2[i] = a2[i] * b2[i];
        }
        trans(c0, true), trans(c1, true), trans(c2, true);
        static const mint mod0 = MOD0, mod01 = mod0 * MOD1;
        vector<mint> res(N + M - 1);
        for (int i = 0; i < N + M - 1; ++i) {
            int y0 = c0[i].val;
            int y1 = (imod0 * (c1[i] - y0)).val;
            int y2 = (imod01 * (c2[i] - y0) - imod1 * y1).val;
            res[i] = mod01 * y2 + mod0 * y1 + y0;
        }
        return res;
    }

    // long long
    vector<long long> mul_ll
    (const vector<long long>& A, const vector<long long>& B) {
        if (A.empty() || B.empty()) return {};
        int N = (int)A.size(), M = (int)B.size();
        if (min(N, M) < 30) return naive_mul(A, B);
        int size_fft = get_fft_size(N, M);
        vector<mint0> a0(size_fft, 0), b0(size_fft, 0), c0(size_fft, 0);
        vector<mint1> a1(size_fft, 0), b1(size_fft, 0), c1(size_fft, 0);
        vector<mint2> a2(size_fft, 0), b2(size_fft, 0), c2(size_fft, 0);
        for (int i = 0; i < N; ++i)
            a0[i] = A[i], a1[i] = A[i], a2[i] = A[i];
        for (int i = 0; i < M; ++i)
            b0[i] = B[i], b1[i] = B[i], b2[i] = B[i];
        trans(a0), trans(a1), trans(a2), trans(b0), trans(b1), trans(b2);
        for (int i = 0; i < size_fft; ++i) {
            c0[i] = a0[i] * b0[i];
            c1[i] = a1[i] * b1[i];
            c2[i] = a2[i] * b2[i];
        }
        trans(c0, true), trans(c1, true), trans(c2, true);
        static const long long mod0 = MOD0, mod01 = mod0 * MOD1;
        vector<long long> res(N + M - 1);
        for (int i = 0; i < N + M - 1; ++i) {
            int y0 = c0[i].val;
            int y1 = (imod0 * (c1[i] - y0)).val;
            int y2 = (imod01 * (c2[i] - y0) - imod1 * y1).val;
            res[i] = mod01 * y2 + mod0 * y1 + y0;
        }
        return res;
    }
};

// Binomial Coefficient
template<class T> struct BiCoef {
    vector<T> fact_, inv_, finv_;
    constexpr BiCoef() {}
    constexpr BiCoef(int n) noexcept : fact_(n, 1), inv_(n, 1), finv_(n, 1) {
        init(n);
    }
    constexpr void init(int n) noexcept {
        fact_.assign(n, 1), inv_.assign(n, 1), finv_.assign(n, 1);
        int MOD = fact_[0].getmod();
        for(int i = 2; i < n; i++){
            fact_[i] = fact_[i-1] * i;
            inv_[i] = -inv_[MOD%i] * (MOD/i);
            finv_[i] = finv_[i-1] * inv_[i];
        }
    }
    constexpr T com(int n, int k) const noexcept {
        if (n < k || n < 0 || k < 0) return 0;
        return fact_[n] * finv_[k] * finv_[n-k];
    }
    constexpr T fact(int n) const noexcept {
        if (n < 0) return 0;
        return fact_[n];
    }
    constexpr T inv(int n) const noexcept {
        if (n < 0) return 0;
        return inv_[n];
    }
    constexpr T finv(int n) const noexcept {
        if (n < 0) return 0;
        return finv_[n];
    }
};

const int MOD = 998244353;
using mint = Fp<MOD>;

using Graph = vector<vector<int>>;
int V, E;
void dfs(const Graph &G, int v, vector<bool> &seen) {
    seen[v] = true;
    ++V;
    for (auto e: G[v]) {
        ++E;
        if (seen[e]) continue;
        dfs(G, e, seen);
    }
}

int main() {
    int N, M;
    cin >> N >> M;
    BiCoef<mint> bc(N+1);
    Graph G(N);
    for (int i = 0; i < M; ++i) {
        int a, b;
        cin >> a >> b;
        --a, --b;
        G[a].push_back(b);
        G[b].push_back(a);
    }

    // 各連結成分ごとに求める
    priority_queue<pair<int,vector<mint>>, vector<pair<int,vector<mint>>>, greater<pair<int,vector<mint>>>> que;
    vector<bool> seen(N, false);
    for (int v = 0; v < N; ++v) {
        if (seen[v]) continue;
        V = 0, E = 0;
        dfs(G, v, seen);
        E /= 2;
        mint two = 1;
        for (int i = 0; i < E-V+1; ++i) two *= 2;
        vector<mint> f(V+1, 0);
        for (int i = 0; i <= V; ++i) {
            if (i % 2 == 1) continue;
            f[i] = bc.com(V, i) * two;
        }
        que.push({f.size(), f});
    }

    // 二分木のような計算順序で FFT
    while (que.size() >= 2) {
        auto f = que.top().second; que.pop();
        auto g = que.top().second; que.pop();
        auto h = NTT::mul(f, g);
        que.push({h.size(), h});
    }
    auto res = que.top().second;
    for (int i = 0; i <= N; ++i) cout << res[i] << endl;
}

けんちょんの競プロ精進記録

競プロの精進記録や小ネタを書いていきます

AtCoder ARC 115 D - Odd Degree (黄色, 600 点)

問題概要

制約

考えたこと

背景

前者 (全域木へと帰着)

後者 (木の場合)

補足

コード