BZOJ3456 城市规划（多项式求逆）

　　设f[i]为连通图的数量，g[i]为不连通图的数量，显然有f[i]=2^i*(i-1)/2-g[i]，g[i]通过枚举1所在连通块大小转移，有g[i]=Σf[j]*C(i-1,j-1)·2^{(i-j)*(i-j-1)/2}，也即f[i]=2^i*(i-1)/2-(i-1)!·Σf[j]·2^{(i-j)*(i-j-1)/2}/(j-1)!/(i-j)!。显然是一个卷积形式，可以分治NTT。

　　进一步将式子化的更优美一点。设h[i]=2^i*(i-1)/2，有f[i]=h[i]-(i-1)!·Σf[j]·h[i-j]/(j-1)!/(i-j)!。阶乘项也可以弄掉，设F[i]=f[i]/(i-1)!，H[i]=h[i]/(i-1)!，G[i]=h[i]/i!，则有F[i]=H[i]-ΣF[j]·G[i-j] (0<j<i)，而G[0]=1，不妨设F[0]=H[0]=0，则H[i]=ΣF[j]·G[i-j] (0<=j<=i)，标准的卷积。设F(x),G(x),H(x)为其各自的生成函数，则H(x)=F(x)·G(x)。H(x)和G(x)我们很容易就能算出来，所以求一发逆就有F(x)了。

　　预定本年度最弱智bug：定义了全局变量inv3预处理3的逆元，在主程序里给其赋值，写了int inv3=……。然后就盯着两份一模一样的代码调了一年。

　　upd：原来这个除了阶乘的东西就叫指数型生成函数？

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
#define N 540000
#define P 1004535809
char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;}
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int n,f[N],g[N],h[N],fac[N],inv[N],r[N],tmp[N],inv3;
int ksm(int a,int k)
{
    int s=1;
    for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P;
    return s;
}
void DFT(int *a,int n,int g)
{
    for (int i=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
    for (int i=2;i<=n;i<<=1)
    {
        int wn=ksm(g,(P-1)/i);
        for (int j=0;j<n;j+=i)
        {
            int w=1;
            for (int k=j;k<j+(i>>1);k++,w=1ll*w*wn%P)
            {
                int x=a[k],y=1ll*w*a[k+(i>>1)]%P;
                a[k]=(x+y)%P,a[k+(i>>1)]=(x-y+P)%P;
            }
        }
    }
}
void mul(int *a,int *b,int n,int op)
{
    for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1);
    DFT(a,n,3),DFT(b,n,3);
    if (op==0) for (int i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%P;
    else for (int i=0;i<n;i++) a[i]=1ll*a[i]*(P+2-1ll*a[i]*b[i]%P)%P;
    DFT(a,n,inv3);
    int u=ksm(n,P-2);
    for (int i=0;i<n;i++) a[i]=1ll*a[i]*u%P;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj3456.in","r",stdin);
    freopen("bzoj3456.out","w",stdout);
    const char LL[]="%I64d
";
#else
    const char LL[]="%lld
";
#endif
    n=read();inv3=ksm(3,P-2);
    fac[0]=fac[1]=1;for (int i=2;i<=n;i++) fac[i]=1ll*fac[i-1]*i%P;
    inv[0]=inv[1]=1;for (int i=2;i<=n;i++) inv[i]=P-1ll*(P/i)*inv[P%i]%P;
    for (int i=2;i<=n;i++) inv[i]=1ll*inv[i-1]*inv[i]%P;
    for (int i=0;i<=n;i++) g[i]=1ll*ksm(2,(1ll*i*(i-1)>>1)%(P-1))*inv[i]%P,h[i]=1ll*g[i]*i%P;
    int t=1;while (t<=(n<<1)) t<<=1;
    f[0]=1;
    for (int i=2;i<=t;i<<=1)
    {
        for (int j=0;j<i;j++) tmp[j]=g[j];
        mul(f,tmp,i<<1,1);
        for (int j=i;j<(i<<1);j++) f[j]=0;
    }
    for (int i=n+1;i<t;i++) f[i]=0;
    mul(f,h,t,0);
    cout<<1ll*f[n]*fac[n-1]%P;
    return 0;
}