BZOJ 3456: 城市规划 [多项式求逆元 DP]

题意：

求出n个点的简单(无重边无自环)无向连通图数目.
方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.

n<=130000

---

DP求方案

g(n) n个点所有图的方案数 显然2^C(n,2)=2^n(n-1)

f(n) n个点连通图的方案数 

然后枚举第一个点所在连通块的点数 

g(n)=∑_i=1..n-1{C(n-1,i-1)*f(i)*g(n-i)}

代入g(n) 两边同除(n-1)!消掉那个组合数上面那块，就变成了卷积的形式

我不写了直接看Miskcoo的公式啦 http://blog.miskcoo.com/2015/05/bzoj-3456

然后C(x)=A(x)*B(x)

A(x)=C(x)*B(x)^-1

放在mod (x^>n) 意义下求逆元就行了 因为需要的是a[n]

---

多项式求逆元

去看Miskcoo的教程吧 http://blog.miskcoo.com/2015/05/polynomial-inverse

简单的思路就是知道A(x) mod (x[n/2]) 下的逆元求mod (xn) 下的逆元

方法就是两个同余的式子写出来一减，两边平方再同乘A(x) 再移项

说一点关于意义的理解吧：

A(x)=Q(x)B(x)+R(x) degR<degB

A(x)Ξ0 (mod xⁿ) 就是说A(x)的0..n-1项系数都是0

A(x)B(x)Ξ1 (mod xⁿ) 它们每一项都有xⁿ,否则不可能余数只有1;所以也有x^n/2;

注意：

1.最后要乘(n-1)! 不要乘(n-1)

2.多项式求逆元每次长度都不确定，不能先预处理二进制反转

 

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=3e5+5;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
int P=1004535809,MOD=P;
ll Pow(ll a,ll b,ll MOD){
    ll ans=1;
    for(;b;b>>=1,a=a*a%MOD)
        if(b&1) ans=ans*a%MOD;
    return ans;
}
struct NTT{
    int n,rev[N];
    ll g;
    void ini(int m){
        n=1;
        while(n<m) n<<=1;
        /*
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            rev[i]=t;
        }
        */
        g=3;
    }
    void transform(int *a,int flag,int n){
        int k=0;
        while((1<<k)<n) k++;
        for(int i=0;i<n;i++){
            int t=0;
            for(int j=0;j<k;j++) if(i&(1<<j)) t|=(1<<(k-j-1));
            if(t<i) swap(a[i],a[t]);
        }
 
        for(int l=2;l<=n;l<<=1){
            int m=l>>1;
            ll wn=Pow(g,flag==1?(P-1)/l:P-1-(P-1)/l,P);
            for(int *p=a;p!=a+n;p+=l){
                ll w=1;
                for(int k=0;k<m;k++){
                    ll t=w*p[k+m]%P;
                    p[k+m]=(p[k]-t+P)%P;
                    p[k]=(p[k]+t)%P;
                    w=w*wn%P;
                }
            }
        }
        if(flag==-1){
            ll inv=Pow(n,P-2,P);
            for(int i=0;i<n;i++) a[i]=a[i]*inv%P;
        }
    }
    int c[N];
    void test(int *a,int n){for(int i=0;i<n;i++) printf("%d ",a[i]);puts("");}
    void polyInv(int deg,int *a,int *b){
        if(deg==1) b[0]=Pow(a[0],P-2,P);
        else{
            polyInv((deg+1)>>1,a,b);
            int n=1;
            while(n< deg<<1) n<<=1;
            copy(a,a+deg,c);
            fill(c+deg,c+n,0);
            transform(c,1,n);
            transform(b,1,n);
            for(int i=0;i<n;i++)
                b[i]=(ll)b[i]*(2-(ll)b[i]*c[i]%P+P)%P;
            transform(b,-1,n);
            fill(b+deg,b+n,0);
        }
    }
}fft;
int n,inv[N],invFac[N],poc[N],A[N],B[N],C[N];
void getInv(int n){
    inv[1]=invFac[0]=1;
    for(int i=1;i<=n;i++){
        if(i!=1) inv[i]=-(ll)P/i*inv[P%i]%P;
        if(inv[i]<0) inv[i]+=P;
        invFac[i]=(ll)invFac[i-1]*inv[i]%P;
    }
}
int main(){
    //freopen("in","r",stdin);
    n=read();
    fft.ini(n);
    getInv(n);
 
    poc[0]=poc[1]=1;
    for(int i=2;i<=n;i++) poc[i]=Pow(2,(ll)i*(i-1)>>1%(P-1),P);
    for(int i=0;i<=n;i++) B[i]=(ll)poc[i]*invFac[i]%P;
    for(int i=1;i<=n;i++) C[i]=(ll)poc[i]*invFac[i-1]%P;//printf("CC %d
",C[i]);
    fft.polyInv(fft.n,B,A);
 
    fft.n<<=1;
    fft.transform(A,1,fft.n);
    fft.transform(C,1,fft.n);
    for(int i=0;i<fft.n;i++) A[i]=(ll)A[i]*C[i]%P;//,printf("ABC %d %d %d
",i,A[i],C[i]);
 
    fft.transform(A,-1,fft.n);
    printf("%lld",(ll)A[n]*Pow(invFac[n-1],P-2,P)%P);
}