[洛谷 5828] 边双连通图计数（生成函数 + 扩展拉格朗日反演） | 错题本

题目

[洛谷 5828] 边双连通图计数

边双联通图：若一个无向图删去任意一条边，这个无向图仍联通，则这个无向图是边双连通图。

分析

同 点双连通图计数，先得到 有根无向连通图 的 EGF $R(x)$，设 有根边双连通图 的 EGF 为 $G(x)$。相对点双而言，边双比较简单，因为根最多属于一个 极大边双 中（否则，根所在的两个边双显然可以合成一个边双），因此只需考虑根 $r$ 所在的边双，然后对于一个以 $r$ 为根的 有根无向连通图，它就是将一些 有根无向连通图 分别通过一条割边挂在了 $r$ 所在的边双上，因此枚举根所在的 极大边双 的大小，由于每个挂上去的 有根无向连通图 可以选择挂在哪个点，于是可以得到 $$\begin{array}{rl}R(x)& =\sum _{i=1}^n{g}_i\frac{x^i}{i!}\exp (iR(x))\\ & =\sum _{i=1}^n{g}_i\frac{x^i}{i!}{\exp }^{i}R(x)\\ & =\sum _{i=1}^n{g}_i\frac{(x\exp R(x){)}^i}{i!}\\ & =G(x\exp R(x))\end{array}$$ 设 $H(x)=x\exp R(x)$，则 $$G(H(x))=R(x)$$ 代入 扩展拉格朗日反演 中可得 $$\begin{array}{rl}[x^n]G(x)& =\frac{1}{n}[x^{n-1}]R^{\prime }(x){\left(\frac{x}{H(x)}\right)}^n\\ & =\frac{1}{n}[x^{n-1}]R^{\prime }(x){\left(\frac{1}{\exp R(x)}\right)}^n\\ & =\frac{1}{n}[x^{n-1}]R^{\prime }(x)\exp (-nR(x))\end{array}$$ 直接算即可，注意我们算的是 EGF 并且是有根的，所以最后要除以 $n$ 乘 $n!$。

代码

变量名有出入。

#include <bits/stdc++.h>

#define RG register

typedef long long LL;

int Read() {
    int x = 0; bool f = false; char c = getchar();
    while (c < '0' || c > '9')
        f |= c == '-', c = getchar();
    while (c >= '0' && c <= '9')
        x = (x * 10) + (c ^ 48), c = getchar();
    return f ? -x : x;
}

template <const int _MOD> struct ModNumber { // 为了效率 (事实上还是不高) 省去了一些实用性
    int x;
    inline ModNumber() { x = 0; }
    inline ModNumber(const int &y) { x = y; }
    // 需保证 y 的范围! (如果在这 % _MOD 会 T, 因为代码中大量调用该构造函数)
    // (当然, 开 O2 可以起飞)
    inline int Int() { return x; }
    inline ModNumber Pow(LL y) const {
        RG int ret = 1, tmp = x;
        while (y) {
            if (y & 1) ret = ((LL)ret * tmp) % _MOD;
            y >>= 1; tmp = ((LL)tmp * tmp) % _MOD;
        }
        return ModNumber(ret);
    }
    inline bool operator == (const ModNumber &y) const { return x == y.x; }
    inline bool operator != (const ModNumber &y) const { return x != y.x; }
    inline bool operator < (const ModNumber &y) const { return x < y.x; }
    inline bool operator > (const ModNumber &y) const { return x > y.x; }
    inline bool operator <= (const ModNumber &y) const { return x <= y.x; }
    inline bool operator >= (const ModNumber &y) const { return x >= y.x; }
    inline ModNumber operator + (const ModNumber &y) const { return (x + y.x >= _MOD) ? (x + y.x - _MOD) : (x + y.x); }
    inline ModNumber operator - (const ModNumber &y) const { return (x - y.x < 0) ? (x - y.x + _MOD) : (x - y.x); }
    inline ModNumber operator * (const ModNumber &y) const { return ModNumber((LL)x * y.x % _MOD); }
    inline ModNumber operator / (const ModNumber &y) const { return *this * y.Pow(_MOD - 2); }
    inline ModNumber operator ^ (const LL &y) const { return Pow(y); }
    inline void operator += (const ModNumber &y) { *this = *this + y; }
    inline void operator *= (const ModNumber &y) { *this = *this * y; }
    inline void operator -= (const ModNumber &y) { *this = *this - y; }
    inline void operator /= (const ModNumber &y) { *this = *this / y; }
    inline void operator ^= (const LL &y) const { *this = *this ^ y; }
};

const int MAXN = 100000 * 4; // 所有 MAXN 都要开 4 倍!
const int MOD = 998244353;

typedef ModNumber<MOD> Int;

const Int __G = 3, One = 1, Two = 2, InvTwo = One / Two;

namespace Polynomial {
    // 好氧, 好氧, 好氧! (无氧可能原地去世)
    // 所有 n: 多项式的项数 (即次数 + 1)
    // 数组从 0 开始存
    int Rev[MAXN + 5];
    Int G0[2][MAXN + 5];

    void GetG0(const int &n) { // 使用前必须先初始化 G0
        for (RG int i = 2; i <= n; i <<= 1) {
            G0[0][i] = __G ^ ((MOD - 1) / i);
            G0[1][i] = G0[0][i] ^ (MOD - 2);
        }
    }

    inline void GetRev(const int n) { // Rev 在函数内初始化
        for (RG int i = 0; i < n; i++)
            Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) * (n >> 1));
    }

    inline int ToPow(const int &n) { // lim 在函数内初始化
        RG int ret = 1;
        while (ret < n)
            ret <<= 1;
        return ret;
    }

    void PrintPoly(Int *A, const int &n) {
        for (RG int i = 0; i < n; i++)
            printf("%d ", A[i].x);
        puts("");
    }

    void ReadPoly(Int *A, const int &n) {
        for (RG int i = 0; i < n; i++)
            A[i].x = Read();
    }

    void NTT(Int *A, const int &n, const int &opt) {
        for (RG int i = 0; i < n; i++)
            if (i < Rev[i])
                std::swap(A[i], A[Rev[i]]);
        for (RG int mid = 1; mid < n; mid <<= 1) {
            const int k = mid << 1;
            const Int g0 = G0[opt][k];
            for (RG int i = 0; i < n; i += k) {
                Int g = 1;
                for (RG int j = 0; j < mid; j++, g *= g0) {
                    Int tmp1 = A[i + j], tmp2 = A[i + j + mid] * g;
                    A[i + j] = tmp1 + tmp2, A[i + j + mid] = tmp1 - tmp2;
                }
            }
        }
        if (opt == 1) {
            const Int inv = One / n;
            for (RG int i = 0; i < n; i++)
                A[i] *= inv;
        }
    }

    Int A0[MAXN + 5], B0[MAXN + 5];

    void Multiply(const Int *A, const Int *B, Int *P, const int &n, const int &m) { // P = A * B
        int lim = ToPow(n + m - 1);
        GetRev(lim);
        for (int i = 0; i < lim; i++)
            A0[i] = B0[i] = 0;
        for (RG int i = 0; i < n; i++)
            A0[i] = A[i];
        for (RG int i = 0; i < m; i++)
            B0[i] = B[i];
        NTT(A0, lim, 0), NTT(B0, lim, 0);
        for (RG int i = 0; i < lim; i++)
            P[i] = A0[i] * B0[i];
        NTT(P, lim, 1);
    }

    Int A1[MAXN + 5], B1[MAXN + 5], C1[MAXN + 5];

    void Inverse(const Int *A, Int *B, const int &n) { // B = 1 / A, A 不变
        if (n == 1) {
            B[0] =  A[0] ^ (MOD - 2);
            return;
        }
        Inverse(A, B, (n + 1) >> 1);
        const int lim = ToPow(n + n - 1);
        GetRev(lim);
        for (RG int i = 0; i < n; i++)
            A1[i] = A[i], B1[i] = B[i];
        for (RG int i = n; i < lim; i++)
            A1[i] = B1[i] = 0;
        NTT(A1, lim, 0), NTT(B1, lim, 0);
        for (RG int i = 0; i < lim; i++)
            C1[i] = A1[i] * B1[i] * B1[i];
        NTT(C1, lim, 1);
        for (RG int i = 0; i < n; i++)
            B[i] = Two * B[i] - C1[i];
        for (RG int i = n; i < lim; i++)
            B[i] = 0;
    }

    Int C2[MAXN + 5], InvB[MAXN + 5], B2[MAXN + 5];

    void Sqrt(const Int *A, Int *B, const int &n) { // B = √A, A 不变 (默认开根的多项式常数项为 1)
        if (n == 1) {
            B[0] = 1;
            return;
        }
        Sqrt(A, B, (n + 1) >> 1);
        Inverse(B, InvB, n);
        Multiply(B, B, B2, n, n);
        for (RG int i = 0; i < n; i++)
            InvB[i] *= InvTwo;
        for (RG int i = 0; i < n; i++)
            B2[i] += A[i];
        Multiply(B2, InvB, C2, n, n);
        for (RG int i = 0; i < n; i++)
            B[i] = C2[i];
    }

    Int AR[MAXN + 5], BR[MAXN + 5], InvBR[MAXN + 5], A2[MAXN + 5], CR[MAXN + 5], C3[MAXN + 5];

    void Divide(const Int *A, const Int *B, Int *C, Int *R, const int &n, const int &m) { // C = A / B, R = A % B, A 不变, B 不变
        for (RG int i = 0; i < n; i++)
            AR[i] = A[n - i - 1], A2[i] = A[i];
        for (RG int i = 0; i < n - m + 1; i++)
            BR[i] = B[m - i - 1];
        Inverse(BR, InvBR, n - m + 1);
        Multiply(AR, InvBR, CR, n, n - m + 1);
        for (int i = 0; i < n - m + 1; i++)
            C[i] = CR[n - m - i];
        Multiply(B, C, C3, m, n - m + 1);
        for (RG int i = 0; i < m - 1; i++)
            R[i] = A[i] - C3[i];
    }

    void Derivative(const Int *A, Int *B, const int &n) { // B = A', A 不变
        for (RG int i = 0; i < n - 1; i++)
            B[i] = A[i + 1] * (i + 1);
        B[n - 1] = 0;
    }

    void Integral(const Int *A, Int *B, const int &n) { // B = ∫A, A 不变
        for (RG int i = 1; i < n; i++)
            B[i] = A[i - 1] / i;
        B[0] = 0;
    }

    Int AD[MAXN + 5], InvA[MAXN + 5], C4[MAXN + 5];

    void Ln(const Int *A, Int *B, const int &n) { // B = ln A, A 不变
        Derivative(A, AD, n);
        Inverse(A, InvA, n);
        Multiply(AD, InvA, C4, n, n);
        Integral(C4, B, n);
    }

    Int LnB[MAXN + 5], B3[MAXN + 5], B4[MAXN + 5];

    void Exp(const Int *A, Int *B, const int &n) { // B = e^A, A 不变
        if (n == 1) {
            B[0] = 1;
            return;
        }
        Exp(A, B, (n + 1) >> 1);
        Ln(B, LnB, n);
        const int lim = ToPow(n + n - 1);
        for (int i = 0; i < n; i++)
            B3[i] = A[i] - LnB[i], B4[i] = B[i];
        B3[0] += One;
        Multiply(B3, B4, B, n, n);
        for (int i = n; i < lim; i++)
            B[i] = 0;
    }

    Int KA[MAXN + 5], LnA[MAXN + 5];

    void Pow(const Int *A, Int *B, const int &n, const Int &k) { // B = A^k, A 不变
        Polynomial::Ln(A, LnA, n);
        for (int i = 0; i < n; i++)
            KA[i] = LnA[i] * k;
        Exp(KA, B, n);
    }

    Int A3[MAXN + 5], AP[MAXN + 5], C5[MAXN + 5];

    Int CompoundInverse(const Int *A, const int &k, const Int *C = NULL) { // A(B(x)) = x 或 A(B(x)) = C(x), 返回 [x^k]B
        for (int i = 0; i < k; i++)
            A3[i] = A[i + 1];
        Pow(A3, AP, k, k);
        Inverse(AP, A3, k);
        if (C == NULL)
            return A3[k - 1] / k;
        Derivative(A3, AP, k);
        Multiply(AP, C, C5, k, k);
        return C5[k - 1] / k;
    }
}

int N;
Int Fac[MAXN + 5], Inv[MAXN + 5];
Int G[MAXN + 5], D[MAXN + 5], D1[MAXN + 5], H[MAXN + 5],  HE[MAXN + 5], H1[MAXN + 5];

int main() {
    Polynomial::GetG0(MAXN);
    N = MAXN / 4;
    Fac[0] = 1;
    for (int i = 1; i <= N; i++)
        Fac[i] = Fac[i - 1] * i;
    Inv[N] = One / Fac[N];
    for (int i = N - 1; i >= 0; i--)
        Inv[i] = Inv[i + 1] * (i + 1);
    int T = 5;
    while (T--) {
        if ((N = Read()) == 1) {
            puts("1");
            continue;
        }
        G[0] = 1;
        for (int i = 1; i <= N; i++)
            G[i] = (Two ^ ((LL)i * (i - 1) / 2)) * Inv[i];
        Polynomial::Ln(G, D, N + 1);
        for (int i = 0; i <= N; i++) {
            D[i] = D[i] * i;
            H[i] = D[i] * (MOD - N);
        }
        Polynomial::Derivative(D, D1, N + 1);
        Polynomial::Exp(H, HE, N + 1);
        Polynomial::Multiply(HE, D1, H1, N + 1, N);
        Int Ans = H1[N - 1] / N * Fac[N - 1];
        printf("%d\n", Ans.x);
    }
    return 0;
}