[bzoj3456]城市规划

题目大意

给你n个点(存在顺序性),初始无边,你可以任意加边。求满足以下条件的连通图数量:无重边无自环。答案模 479221+1

DP

我们设f[i]表示i个点对应的答案。
正难则反,我们可以用总数量减去不合法数量。
总数量为 2C2i
如何不重复不遗漏统计不合法数量?
我们可以枚举i点所在连通块的数量,然后其他乱连,
也就是 f[i]=2C2ii1j=12C2ijCj1i1f[j]
注意到 C2i=i(i1)/2
我们预处理 v[i]=2i(i1)/2
那么 f[i]=v[i]i1j=1v[ij]Cj1i1f[j]
可以降一个log n的复杂度。
总复杂度n^2。

往FFT方面进军

我们设 ff[i]=i!,gg[i]=(i!)ff[i]
Cxy=ff[y]gg[x]gg[yx]
那么 f[i]=v[i]i1j=1v[ij]ff[i1]gg[j1]gg[ij]f[j]
我们设 u[i]=v[i]gg[i],g[i]=f[i]gg[i1],h[i]=v[i]gg[i1]
f[i]=v[i]ff[i1]i1j=1u[ij]g[j]
注意到 f[i]=f[i]u[0]=f[i]gg[i1]ff[i1]u[0]=u[0]g[i]
所以两边同时减去f[i]
v[i]ff[i1]ij=1u[ij]g[j]=0
两边同乘gg[i-1]并移项,得到标准卷积形式。
h=ug
h与u我们是可以处理出来的,要求的多项式为g
g=hu1
我们要做的是求多项式u的逆元,然后进行有取模运算的FFT。

带取模运算的FFT

这题的FFT需要进行取模,于是我们没必要再用复数运算,那该怎么办呢?
原来的w数组满足以下几个性质,我们只要满足这个性质,FFT的分治策略便依然正确:
1、w[0]=w[n]=1
2、w[i]!=w[j]当i!=j时
3、w[i]=w[i-1]*w[1]
解决方法如下:
设GG为mo的原根(本题GG=3),原根是什么呢?
假如a是m的原根,那么满足m-1是最小满足 axMod=1 的x值。
那好办了,根据定义, w[i]=GG(mo1)i/n
注意只有mo-1满足一定是n的倍数时才可以这样做。这题mo-1中含有2的21次幂,n一定是2的幂数且不大于21次幂,所以本题可以使用带取模运算的FFT。
接下来,我们需要解决如何求多项式的逆元。
代码如下:

void DFT(ll *a,ll sig){
    fo(i,0,len-1){
        ll p=0;
        for(ll j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
        tt[p]=a[i];
    }
    for(ll m=2;m<=len;m*=2){
        ll half=m/2,bei=len/m;
        fo(i,0,half-1){
            ll wi=(sig>0)?w[i*bei]:w[len-i*bei];
            for(ll j=i;j<len;j+=m){
                ll u=tt[j],v=tt[j+half]*wi%mo;
                tt[j]=(u+v)%mo;
                tt[j+half]=((u-v)%mo+mo)%mo;
            }
        }
    }
    if (sig==-1)
        fo(i,0,len-1) tt[i]=tt[i]*ni%mo;
    fo(i,0,len-1) a[i]=tt[i];
}
void FFT(ll *a,ll *b,ll *c,ll x){
    ll i;
    len=x*2;
    fo(i,0,len) e[i]=f[i]=0;
    fo(i,0,x) e[i]=a[i],f[i]=b[i];
    ni=quicksortmi(len,mo-2);
    w[1]=quicksortmi(GG,(mo-1)/len);
    w[0]=1;
    fo(i,2,len) w[i]=w[i-1]*w[1]%mo;
    ce=double(log(len)/log(2));
    DFT(e,1);DFT(f,1);
    fo(i,0,len-1) e[i]=e[i]*f[i]%mo;
    DFT(e,-1);
    fo(i,0,len-1) c[i]=e[i];
}

FFT(a,b,c,x)表示次数界限为x的两个多项式a与b相乘,得到次数界限为x*2的多项式c。

多项式求逆元

多项式求逆元同样用到了倍增思想。
假设A的逆元是B,则
AB1(mod xn)
如果我们已经求出在 A 在模xn2的逆元 B0
AB01(mod xn2)
因为 AB1(mod xn)
所以 AB1(mod xn2)
那么就有 BB00(mod xn2)
两边平方并移项
B2B0(2BB0)(mod xn)
两边同乘 A
BB0(2B0A)(mod xn)
代码如下

void getny(ll *a,ll x){
    if (x==1){
        ny[0]=1;
        return;
    }
    getny(a,x/2);
    FFT(ny,ny,b,x/2);
    FFT(b,a,b,x);
    fo(i,0,x-1) ny[i]=((2*ny[i]%mo-b[i])%mo+mo)%mo;
}

注意FFT时次数界限应和当前x有关。
这样本题就可以解决了。

参考程序

#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const ll maxn=130000+10,mo=1004535809,GG=3;
ll g[maxn],u[maxn],h[maxn],ff[maxn],gg[maxn],v[maxn],w[maxn*5];
ll i,j,k,l,t,n,m,len,wdc,ni;
double ce;
ll e[maxn*5],f[maxn*5],a[maxn*5],b[maxn*5],c[maxn*5],tt[maxn*5],ny[maxn*5];
ll quicksortmi(ll x,ll y){
    if (!y) return 1;
    ll t=quicksortmi(x,y/2);
    t=t*t%mo;
    if (y%2) t=t*(x%mo)%mo;
    return t;
}
void DFT(ll *a,ll sig){
    fo(i,0,len-1){
        ll p=0;
        for(ll j=0,tp=i;j<ce;j++,tp/=2) p=(p<<1)+(tp%2);
        tt[p]=a[i];
    }
    for(ll m=2;m<=len;m*=2){
        ll half=m/2,bei=len/m;
        fo(i,0,half-1){
            ll wi=(sig>0)?w[i*bei]:w[len-i*bei];
            for(ll j=i;j<len;j+=m){
                ll u=tt[j],v=tt[j+half]*wi%mo;
                tt[j]=(u+v)%mo;
                tt[j+half]=((u-v)%mo+mo)%mo;
            }
        }
    }
    if (sig==-1)
        fo(i,0,len-1) tt[i]=tt[i]*ni%mo;
    fo(i,0,len-1) a[i]=tt[i];
}
void FFT(ll *a,ll *b,ll *c,ll x){
    ll i;
    len=x*2;
    fo(i,0,len) e[i]=f[i]=0;
    fo(i,0,x) e[i]=a[i],f[i]=b[i];
    ni=quicksortmi(len,mo-2);
    w[1]=quicksortmi(GG,(mo-1)/len);
    w[0]=1;
    fo(i,2,len) w[i]=w[i-1]*w[1]%mo;
    ce=double(log(len)/log(2));
    DFT(e,1);DFT(f,1);
    fo(i,0,len-1) e[i]=e[i]*f[i]%mo;
    DFT(e,-1);
    fo(i,0,len-1) c[i]=e[i];
}
void getny(ll *a,ll x){
    if (x==1){
        ny[0]=1;
        return;
    }
    getny(a,x/2);
    FFT(ny,ny,b,x/2);
    FFT(b,a,b,x);
    fo(i,0,x-1) ny[i]=((2*ny[i]%mo-b[i])%mo+mo)%mo;
    //fo(i,0,x-1) printf("%lld ",ny[i]);
    //printf("\n");
}
int main(){
    //freopen("D:/out1.txt","w",stdout);
    scanf("%lld",&n);
    ff[0]=1;
    fo(i,1,n) ff[i]=ff[i-1]*i%mo;
    gg[n]=quicksortmi(ff[n],mo-2);
    fd(i,n-1,0) gg[i]=gg[i+1]*(i+1)%mo;
    fo(i,0,n) v[i]=quicksortmi(2,i*(i-1)/2);
    fo(i,0,n) u[i]=v[i]*gg[i]%mo;
    fo(i,1,n) h[i]=v[i]*gg[i-1]%mo;
    len=1;
    while (len<n*2) len*=2;
    wdc=len;
    getny(u,wdc/2);
    FFT(h,ny,g,wdc/2);
    printf("%lld\n",g[n]*ff[n-1]%mo);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值