【bzoj3456】城市规划 【FFT/NTT】【多项式求逆】

题目大意:你需要求出n个点的简单(无重边无自环)无向连通图的数目模1004535809。 n130000 n ≤ 130000
题解:我们令 fn f n 表示n个点的简单无向图的数目, gn g n 表示n个点的简单无向连通图的数目。
易得 fn=2C2n f n = 2 C n 2 ,意思是任意两点之间的边都可以取或不取。
同时, fn=ni=1Ci1n1gifni f n = ∑ i = 1 n C n − 1 i − 1 g i f n − i ,意思是枚举点1所在连通块的大小i,从n-1个点里再选i-1个,连边组成一个联通块,剩下n-i个点随意连边。
我们变一下这个式子。
2C2n=ni=1Ci1n1gi2C2ni 2 C n 2 = ∑ i = 1 n C n − 1 i − 1 g i 2 C n − i 2
=> 2C2n=ni=1(n1)!(i1)!(ni)!gi2C2ni 2 C n 2 = ∑ i = 1 n ( n − 1 ) ! ( i − 1 ) ! ( n − i ) ! g i 2 C n − i 2
=> 2C2n=(n1)!ni=1(1(i1)!gi)(1(ni)!2C2ni) 2 C n 2 = ( n − 1 ) ! ∑ i = 1 n ( 1 ( i − 1 ) ! g i ) ( 1 ( n − i ) ! 2 C n − i 2 )
=> 2C2n(n1)!=ni=1(gi(i1)!)(2C2ni(ni)!) 2 C n 2 ( n − 1 ) ! = ∑ i = 1 n ( g i ( i − 1 ) ! ) ( 2 C n − i 2 ( n − i ) ! )
可以发现这是一个卷积的形式。
我们令
A(x)=ni=12C2i(i1)!xi A ( x ) = ∑ i = 1 n 2 C i 2 ( i − 1 ) ! x i
B(x)=ni=1gi(i1)!xi B ( x ) = ∑ i = 1 n g i ( i − 1 ) ! x i
C(x)=n1i=02C2ii!xi C ( x ) = ∑ i = 0 n − 1 2 C i 2 i ! x i
A(x)B(x)×C(x) A ( x ) ≡ B ( x ) × C ( x )
=> B(x)A(x)C1(x)(mod xn+1) B ( x ) ≡ A ( x ) ∗ C − 1 ( x ) ( m o d   x n + 1 )
因此我们可以通过多项式求逆和多项式乘法得到B。
多项式求逆的关键公式:
A(x)B(x)1(mod xn) A ( x ) B ( x ) ≡ 1 ( m o d   x n )
A(x)B(x)10(mod xn) A ( x ) B ( x ) − 1 ≡ 0 ( m o d   x n )
(A(x)B(x)1)20(mod x2n) ( A ( x ) B ( x ) − 1 ) 2 ≡ 0 ( m o d   x 2 n )
A(x)(2B(x)B(x)2A(x))1(mod x2n) A ( x ) ( 2 B ( x ) − B ( x ) 2 A ( x ) ) ≡ 1 ( m o d   x 2 n )
最后的答案为 B[n](n1)! B [ n ] ∗ ( n − 1 ) !
代码

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=270005;
const ll mod=1004535809;
int n,m,rev[N];
ll jc[N],a[N],b[N],c[N],d[N];
ll fastpow(ll a,ll x){
    a%=mod;
    ll res=1;
    while(x){
        if(x&1){
            res=res*a%mod;
        }
        x>>=1;
        a=a*a%mod;
    }
    return res;
}
ll getinv(ll a){
    return fastpow(a,mod-2);
}
void ntt(ll *a,int n,int dft){
    for(int i=0;i<n;i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)*(n>>1));
        if(i<rev[i]){
            swap(a[i],a[rev[i]]);
        }
    }
    for(int i=1;i<n;i<<=1){
        ll wn=fastpow(3,(mod-1)/i/2);
        if(dft==-1){
            wn=getinv(wn);
        }
        for(int j=0;j<n;j+=i<<1){
            ll w=1,x,y;
            for(int k=j;k<j+i;k++,w=w*wn%mod){
                x=a[k];
                y=w*a[k+i]%mod;
                a[k]=(x+y)%mod;
                a[k+i]=(x-y+mod)%mod;
            }
        }
    }
    if(dft==-1){
        ll inv=getinv(n);
        for(int i=0;i<n;i++){
            a[i]=a[i]*inv%mod;
        }
    }
}
int main(){
    scanf("%d",&m);
    jc[0]=1;
    for(int i=1;i<=m;i++){
        jc[i]=jc[i-1]*i%mod;
    }
    for(n=1;n<=m;n<<=1);
    a[0]=1;
    for(int i=1;i<m;i++){
        a[i]=fastpow(2,1LL*i*(i-1)/2)*getinv(jc[i]);
    }
    b[0]=getinv(a[0]);
    for(int k=2;k<=n;k<<=1){
        for(int i=0;i<k;i++){
            c[i]=a[i];
            if(i<(k>>1)){
                d[i]=b[i];
            }else{
                d[i]=0;
            }
        }
        ntt(c,k<<1,1);
        ntt(d,k<<1,1);
        for(int i=0;i<(k<<1);i++){
            c[i]=(2*d[i]%mod-d[i]*d[i]%mod*c[i]%mod+mod)%mod;
        }
        ntt(c,k<<1,-1);
        ntt(d,k<<1,-1);
        for(int i=0;i<k;i++){
            b[i]=c[i];
        }
    }
    a[0]=0;
    for(int i=1;i<=m;i++){
        a[i]=fastpow(2,1LL*i*(i-1)/2)*getinv(jc[i-1]);
    }
    ntt(a,n<<1,1);
    ntt(b,n<<1,1);
    for(int i=0;i<(n<<1);i++){
        a[i]=a[i]*b[i]%mod;
    }
    ntt(a,n<<1,-1);
    printf("%lld\n",a[m]*jc[m-1]%mod);
    return 0;
}
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值