bzoj 3456: 城市规划 NTT+多项式求逆

题意

求出n个点的简单(无重边无自环)无向连通图数目。
你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.
n<=130000

分析

首先要知道递推式

$$f[i]=2^{c_i^2}-\sum_{j=1}^{i-1}f[j]*C_{i-1}^{j-1}*2^{C_{i-j}^2}$$
如果不优化直接算的话显然是不能够接受的。考虑化简式子。
将组合数展开，两边同除 $(i-1)!$ $$\frac{f[i]}{(i-1)!}=\frac{2^{C_i^2}}{(i-1)!}-\sum_{j=1}^{i-1}\frac{f[j]}{(j-1)!}*\frac{2^{C_{i-j}^2}}{(i-j)!}$$ $$\frac{2^{C_i^2}}{(i-1)!}=\sum_{j=1}^i\frac{f[j]}{(j-1)!}\frac{2^{C_{i-j}^2}}{(i-j)!}$$
容易发现这是一个卷积的形式。
设 $A[i]=\frac{f[i]}{(i-1)!},B[i]=\frac{2^{C_i^2}}{i!},C[i]=\frac{2^{C_i^2}}{(i-1)!}$
那么有 $C=A*B=>A=C*B^{-1}$
上多项式求逆即可。

多项式求逆：
现要求$A*G=1(\mod x^m)$
已求出$B$满足$A*B=1(\mod x^{\frac{m}{2}})$
因为$A*G=1(\mod x^{\frac{m}{2}})$
所以$(G-B)=0(\mod x^{\frac{m}{2}})$
两边平方得$G^2+B^2-2G*B=0(\mod x^m)$
$G^2=2G*B-B^2(\mod x^m)$
两边同乘A得$G=2B-AB^2$

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long LL;

const int N=300005;
const int MOD=1004535809;

int n,b[N],c[N],ny_b[N],jc[N],ny[N],rev[N],tmp[N];

int ksm(int x,int y)
{
    int ans=1;
    while (y)
    {
        if (y&1) ans=(LL)ans*x%MOD;
        x=(LL)x*x%MOD;y>>=1;
    }
    return ans;
}

void fft(int *a,int n,int f)
{
    int lg=log(n)/log(2)+0.1;
    for (int i=0;i<n;i++)
    {
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
        if (i<rev[i]) swap(a[i],a[rev[i]]);
    }
    for (int i=1;i<n;i*=2)
    {
        int wn=ksm(3,f==1?(MOD-1)/i/2:MOD-1-(MOD-1)/i/2);
        for (int j=0;j<n;j+=i*2)
        {
            int w=1;
            for (int k=0;k<i;k++)
            {
                int u=a[j+k],v=(LL)a[j+k+i]*w%MOD;
                a[j+k]=(u+v)%MOD;a[j+k+i]=(u-v+MOD)%MOD;
                w=(LL)w*wn%MOD;
            }
        }
    }
    if (f==-1)
    {
        int w=ksm(n,MOD-2);
        for (int i=0;i<n;i++) a[i]=(LL)a[i]*w%MOD;
    }
}

void get_ny(int *a,int *b,int n)
{
    if (n==1)
    {
        b[0]=ksm(a[0],MOD-2);
        return;
    }
    get_ny(a,b,n/2);
    memcpy(tmp,a,sizeof(a[0])*n);
    memset(tmp+n,0,sizeof(tmp[0])*n);
    fft(tmp,n<<1,1);fft(b,n<<1,1);
    for (int i=0;i<n*2;i++) tmp[i]=(LL)b[i]*(2-(LL)tmp[i]*b[i]%MOD+MOD)%MOD;
    fft(tmp,n<<1,-1);
    for (int i=0;i<n;i++) b[i]=tmp[i];
    memset(b+n,0,sizeof(b[0])*n);
}

int main()
{
    scanf("%d",&n);
    ny[0]=jc[0]=1;
    for (int i=1;i<=n;i++) jc[i]=(LL)jc[i-1]*i%MOD,ny[i]=ksm(jc[i],MOD-2);
    for (int i=0;i<=n;i++) b[i]=(LL)ksm(2,(LL)i*(i-1)/2%(MOD-1))*ny[i]%MOD;
    for (int i=1;i<=n;i++) c[i]=(LL)ksm(2,(LL)i*(i-1)/2%(MOD-1))*ny[i-1]%MOD;
    int m;
    for (m=1;m<=n;m*=2);
    get_ny(b,ny_b,m);
    fft(c,m<<1,1);fft(ny_b,m<<1,1);
    for (int i=0;i<m*2;i++) c[i]=(LL)c[i]*ny_b[i]%MOD;
    fft(c,m<<1,-1);
    printf("%d",(LL)c[n]*jc[n-1]%MOD);
    return 0;
}