求$N$个点组成的有标号简单无向连通图个数，答案模 10^9+7 1998585857

由这个问题引入，顺便总结一下多项式求逆的内容。

首先可以设 $f(n)$ 表示有 $n$ 个点的有标号简单连通无向图的个数， $g(n)$ 表示有 $n$ 个点的有标号简单无向图的个数（也就是不要求连通）

显然
$$g(n) = 2 ^ {C_n^2}$$

又因为一个有标号简单无向图是由很多连通分量组成的，为了避免重复计数，我们枚举点 $1$ 所在的连通块大小（其余的点随便连边，因为 $1$ 号点所在连通块已经确定，其它怎么连都不会重复）

$$ g(n) = \sum_{i=1}^n C_{n-1}^{i-1}f(i)g(n-i) $$

把组合数写成阶乘的形式并移项：

$$ \frac{g(n)}{(n-1)!} = \sum_{i=1}^n \frac{f(i)}{(i-1)!} * \frac{g(n-i)}{(n-i)!} $$

这就是一个很显然的卷积形式了，利用母函数的知识，考虑这么三个多项式：

$$ A(x) = \sum_{i>=1} \frac{g(i)}{(i-1)!} * x^i $$

$$ B(x) = \sum_{i>=1} \frac{f(i)}{(i-1)!} * x^i $$

$$ C(x) = \sum_{i>=0} \frac{g(i)}{i!} * x^i $$

$$ A = B \otimes C$$

那么只要求出 $A$ 在模 $x^{n+1}$ 意义下的逆元，乘上 $C$ ，就可以得到答案

复杂度$O(NlogN)$

模$10^9+7$的话，现在直接用MTT做就好了，不再需要NTT+CRT合并这种鬼畜套路了，不过我对MTT理解不深，只会套板，多项式求逆写的比较丑，如果是模$1998585857=953*2^{21}+1$, 原根为$3$（关于原根，可以看ACdreamer的原根讲解,以及鄙队维护的gitbook有素数原根表）,直接用NTT做就好了，可以写的比较优雅。

然后总结一下多项式求逆。

多项式求逆元

基本概念

在介绍多项式的逆元之前，先说明一些概念：多项式的度、多项式的逆元、多项式的除法和取余

对于一个多项式$A(x)$，称其最高项的次数为这个多项式的度（degree），记作$degA$

对于多项式$A(x),B(x)$存在唯一的$Q(x)$满足$A(x) = Q(x)B(x) + R(x)$，其中$degR < degB$，我们称$Q(x)$为 $B(x)$ 除 $A(x)$ 的商，$R(x)$ 为 $B(x)$ 除 $A(x)$ 的余数，可以记做

$$A(x) \equiv R(x)(modB(x))$$

多项式的逆元

对于一个多项式$A(x)$，如果存在$B(x)$满足$degB \leq degA $并且

$$ A(x)B(x) \equiv 1\ (mod \ x^n) $$

那么称 $B(x)$ 为 $A(x)$ 在 $mod \ x^n$ 意义下的逆元（inverse element），记作$A^{-1}(x)$

多项式逆元求法

这里我们将使用一种倍增的思想来完成。

首先，当$n=1 $时， $B[0]=A[0]^{-1}$.

然后，假如我们已经求得在模$x^{\lceil t/2 \rceil}$意义下的逆元$B_0(x)$，现在要求在模$x^t$意义下的逆元。

那么
$$A(x) * B_0(x) \equiv 1 \ (mod \ x^{\lceil t/2 \rceil}) $$

$$A(x) * B(x) \equiv 1 \ (mod \ x^{\lceil t/2 \rceil}) $$

由于 A(x) 逆元 B0(x) 存在，故得：

$$B(x) - B_0(x) \equiv 0 \ (mod \ x^{\lceil t/2 \rceil}) $$

那么我们将等式各部分平方，展开后得：

$$B^2(x) - 2 B(x) * B_0(x) + B_0^2(x) \equiv 0 \ (mod \ x^t) $$

两边同乘 A(x)，利用逆元的性质移项便可得：

$$B(x) \equiv B_0(x) (2 - A(x) * B_0(x)) \ (mod \ x^t) $$

上式便可在$O(tlogt)$时间内求出了。

$T(n)=T(n^2)+O(nlogn)$

由主定理易得$T(n) = O(nlogn)$，不过也需要注意，这里其实有很大常数的……

顺便一提，由这个过程可以看出，一个多项式有没有逆元完全取决于其常数项是否有逆元

PS: 对于多项式开方，即求$求B(x)$，使得$B(x) * B(x)=A(x) \ (mod\ x^m)$，同样由倍增的方式，简单分析可得$B(x) = \dfrac{B_0^2(x)+A(x)}{2B_0(x)}$，需要求$B(x)$的逆，常数更大……有生之年系列

代码实现

为了便于理解，可以先来一个裸的实数版本的(实际上应用不多):

注意一下实现技巧……在单次运算时可以只用3次FFT……

\\ cp为复数结构体或者STL::complex
\\ 假设已经有了计算快速傅里叶变换的函数 void fft(cp x[], int n, bool flag);
cp tmp[maxn];
void polyInv(int a[], int b[], int deg) {
    if (deg == 1) { b[0] = cp(1.0) / a[0]; return; }
    polyInv(a, b, (deg + 1) >> 1);

    int p = 1; while(p < deg << 1) p <<= 1;
    copy(a, a + deg, tmp);
    fill(tmp + deg, tmp + p, cp(0.0) );
    fft(tmp, p, 1);
    fft(b, p, 1);
    for (int i = 0; i < p; i++) 
        b[i] *= cp(2.0) - tmp[i] * b[i];
    fft(b, p, -1);
    for (int i = 0; i < p; i++)
        b[i] /= p;
    fill(b + deg, b + o, cp(0.0) );
}

下面是别人的NTT版本代码……主要是我自己现在都是直接MTT没怎么用NTT……

而且NTT版本和上面FFT版本的区别很小，求个逆元就好了。

#include <cstdio>
#include <complex>
#include <cmath>
#include <algorithm>
#include <iostream>
using std::copy;
using std::fill;

const long long mod_v = 17ll * (1 << 27) + 1;
const int MaxN = 10010;
long long a[MaxN], b[MaxN], c[MaxN];
long long eps[MaxN], inv_eps[MaxN];
int tot;

long long power(long long x, long long p)
{
    long long v = 1;
    while(p)
    {
        if(p & 1) v = x * v % mod_v;
        x = x * x % mod_v;
        p >>= 1;
    }

    return v;
}

void init_eps(int n)
{
    tot = n;
    long long base = power(3, (mod_v - 1) / n);
    long long inv_base = power(base, mod_v - 2);
    eps[0] = 1, inv_eps[0] = 1;
    for(int i = 1; i < n; ++i)
    {
        eps[i] = eps[i - 1] * base % mod_v;
        inv_eps[i] = inv_eps[i - 1] * inv_base % mod_v;
    }
}

long long inc(long long x, long long d) 
{
    x += d; 
    return x >= mod_v ? x - mod_v : x; 
}

long long dec(long long x, long long d) 
{
    x -= d; 
    return x < 0 ? x + mod_v : x; 
}

void transform(int n, long long *x, long long *w)
{
    for(int i = 0, j = 0; i != n; ++i)
    {
        if(i > j) std::swap(x[i], x[j]);
        for(int l = n >> 1; (j ^= l) < l; l >>= 1);
    }

    for(int i = 2; i <= n; i <<= 1)
    {
        int m = i >> 1;
        for(int j = 0; j < n; j += i)
        {
            for(int k = 0; k != m; ++k)
            {
                long long z = x[j + m + k] * w[tot / i * k] % mod_v;
                x[j + m + k] = dec(x[j + k], z);
                x[j + k] = inc(x[j + k], z);
            }
        }
    }
}

void polynomial_inverse(int deg, long long* a, long long* b, long long* tmp)
{
    if(deg == 1)
    {
        b[0] = power(a[0], mod_v - 2);
    } else {
        polynomial_inverse((deg + 1) >> 1, a, b, tmp);

        int p = 1;
        while(p < deg << 1) p <<= 1;
        copy(a, a + deg, tmp);
        fill(tmp + deg, tmp + p, 0);
        transform(p, tmp, eps);
        transform(p, b, eps);
        for(int i = 0; i != p; ++i)
        {
            b[i] = (2 - tmp[i] * b[i] % mod_v) * b[i] % mod_v;
            if(b[i] < 0) b[i] += mod_v;
        }
        transform(p, b, inv_eps);
        long long inv = power(p, mod_v - 2);
        for(int i = 0; i != p; ++i)
            b[i] = b[i] * inv % mod_v;
        fill(b + deg, b + p, 0);

    }
}

int main()
{
    init_eps(2048);
    int n;
    std::cin >> n;
    for(int i = 0; i != n; ++i)
        std::cin >> a[i];
    polynomial_inverse(n, a, b, c);
    std::cout << "inverse: ";
    for(int i = 0; i != n; ++i)
        printf("%lld ", (b[i] + mod_v) % mod_v);
    std::cout << std::endl;
    return 0;
}

Update：把自己魔改的TLS的MTT板子放上来……

（其实也就是修正了一些错误加上了求逆）……

// TO be modified
const int maxn = 1e5 + 7;
const int maxLen = 18, maxm = 1 << maxLen | 1;
const ll maxv = 1e10 + 6; // 1e14, 1e15
const DB pi = acos(-1.0); // double is enough
ll mod = 313, nlim, sp, msk;
//

#define _ %mod
#define __ %=mod

ll qpow(ll x, ll p) {
    ll ret = 1;
    while (p) {
        if (p & 1) (ret *= x) __;
        (x *= x) __;
        p >>= 1;
    }
    return ret;
}

namespace FFT{
    struct cp {
        DB r, i;
        cp() {}
        cp(DB r, DB i) : r(r), i(i) {}
        cp operator + (cp const &t) const { return cp(r + t.r, i + t.i); }
        cp operator - (cp const &t) const { return cp(r - t.r, i - t.i); }
        cp operator * (cp const &t) const { return cp(r * t.r - i * t.i, r * t.i + i * t.r); }
        cp conj() const { return cp(r, -i); }
    } w[maxm], wInv[maxm];

    void init() {
        for(int i = 0, ilim = 1 << maxLen; i < ilim; ++i) {
            int j = i, k = ilim >> 1; // 2 pi / ilim
            for( ; !(j & 1) && !(k & 1); j >>= 1, k >>= 1);
            w[i] = cp(cos(pi / k * j), sin(pi / k * j));
            wInv[i] = w[i].conj();
        }
        nlim = std::min(maxv / (mod - 1) / (mod - 1), maxn - 1LL);
        for(sp = 1; 1 << (sp << 1) < mod; ++sp);
        msk = (1 << sp) - 1;
    }

    void FFT(int n, cp a[], int flag) {
        static int bitLen = 0, bitRev[maxm] = {};
        if(n != (1 << bitLen)) {
            for(bitLen = 0; 1 << bitLen < n; ++bitLen);
            for(int i = 1; i < n; ++i)
                bitRev[i] = (bitRev[i >> 1] >> 1) | ((i & 1) << (bitLen - 1));
        }
        for(int i = 0; i < n; ++i)
            if(i < bitRev[i])
                std::swap(a[i], a[bitRev[i]]);
        for(int i = 1, d = 1; d < n; ++i, d <<= 1)
            for(int j = 0; j < n; j += d << 1)
                for(int k = 0; k < d; ++k) {
                    cp &AL = a[j + k], &AH = a[j + k + d];
                    cp TP = w[k << (maxLen - i)] * AH;
                    AH = AL - TP, AL = AL + TP;
                }
        if(flag != -1)
            return;
        std::reverse(a + 1, a + n);
        for(int i = 0; i < n; ++i) {
            a[i].r /= n;
            a[i].i /= n;
        }
    }

    void polyMul(int a[], int aLen, int b[], int bLen, int c[]) { // c not in {a, b}
        static cp A[maxm], B[maxm], C[maxm], D[maxm];
        int len, cLen = aLen + bLen - 1; // optional: parameter
        for(len = 1; len < aLen + bLen - 1; len <<= 1);
        if(std::min(aLen, bLen) <= nlim) {
            for(int i = 0; i < len; ++i)
                A[i] = cp(i < aLen ? a[i] : 0, i < bLen ? b[i] : 0);
            FFT(len, A, 1);
            cp tr(0, -0.25);
            for(int i = 0, j; i < len; ++i)
                j = (len - i) & (len - 1), B[i] = (A[i] * A[i] - (A[j] * A[j]).conj()) * tr;
            FFT(len, B, -1);
            for(int i = 0; i < cLen; ++i) c[i] = (ll)(B[i].r + 0.5) % mod;
            return;
        } // if min(aLen, bLen) * mod <= maxv
        for(int i = 0; i < len; ++i) {
            A[i] = i < aLen ? cp(a[i] & msk, a[i] >> sp) : cp(0, 0);
            B[i] = i < bLen ? cp(b[i] & msk, b[i] >> sp) : cp(0, 0);
        }
        FFT(len, A, 1), FFT(len, B, 1);
        cp trL(0.5, 0), trH(0, -0.5), tr(0, 1);
        for(int i = 0, j; i < len; ++i) {
            j = (len - i) & (len - 1);
            cp AL = (A[i] + A[j].conj()) * trL;
            cp AH = (A[i] - A[j].conj()) * trH;
            cp BL = (B[i] + B[j].conj()) * trL;
            cp BH = (B[i] - B[j].conj()) * trH;
            C[i] = AL * (BL + BH * tr);
            D[i] = AH * (BL + BH * tr);
        }
        FFT(len, C, -1), FFT(len, D, -1);
        for(int i = 0; i < cLen; ++i) {
            int v11 = (ll)(C[i].r + 0.5) % mod, v12 = (ll)(C[i].i + 0.5) % mod;
            int v21 = (ll)(D[i].r + 0.5) % mod, v22 = (ll)(D[i].i + 0.5) % mod;
            c[i] = (((((ll)v22 << sp) + v12 + v21) << sp) + v11) % mod;
        }
    }

    int c[maxm], tmp[maxm];
    // y should clear to 0
    void polyInv(int x[], int y[], int deg) {
        if (deg == 1) {
            y[0] = qpow(x[0], mod - 2);
            return;
        }
        polyInv(x, y, (deg + 1) >> 1);

        copy(x, x + deg, tmp);
        int p = ((deg + 1) >> 1) + deg - 1;
        polyMul(y, (deg + 1) >> 1, tmp, deg, c);

        for (int i = 0; i < p; i += 1) c[i] = (- c[i] + mod) _;
        (c[0] += 2) __;

        polyMul(y, (deg + 1) >> 1, c, deg, tmp);
        copy(tmp, tmp + deg, y);
    }
};