[bzoj3456]城市规划

题目大意

给你n个点（存在顺序性），初始无边，你可以任意加边。求满足以下条件的连通图数量：无重边无自环。答案模479∗221+1

DP

我们设f[i]表示i个点对应的答案。

正难则反，我们可以用总数量减去不合法数量。

总数量为2C2i

如何不重复不遗漏统计不合法数量？

我们可以枚举i点所在连通块的数量，然后其他乱连，

也就是f[i]=2C2i−∑i−1j=12C2i−j∗Cj−1i−1∗f[j]

注意到C2i=i∗(i−1)/2

我们预处理v[i]=2i∗(i−1)/2

那么f[i]=v[i]−∑i−1j=1v[i−j]∗Cj−1i−1∗f[j]

可以降一个log n的复杂度。

总复杂度n^2。

往FFT方面进军

我们设ff[i]=i!,gg[i]=(i!)′即ff[i]的逆元

则Cxy=ff[y]∗gg[x]∗gg[y−x]

那么f[i]=v[i]−∑i−1j=1v[i−j]∗ff[i−1]∗gg[j−1]∗gg[i−j]∗f[j]

我们设u[i]=v[i]∗gg[i],g[i]=f[i]∗gg[i−1],h[i]=v[i]∗gg[i−1]

则f[i]=v[i]−ff[i−1]∗∑i−1j=1u[i−j]∗g[j]

注意到f[i]=f[i]∗u[0]=f[i]∗gg[i−1]∗ff[i−1]∗u[0]=u[0]∗g[i]

所以两边同时减去f[i]

v[i]−ff[i−1]∗∑ij=1u[i−j]∗g[j]=0

两边同乘gg[i-1]并移项，得到标准卷积形式。

h=u∗g

h与u我们是可以处理出来的，要求的多项式为g

g=h∗u−1

我们要做的是求多项式u的逆元，然后进行有取模运算的FFT。

带取模运算的FFT

这题的FFT需要进行取模，于是我们没必要再用复数运算，那该怎么办呢？

原来的w数组满足以下几个性质，我们只要满足这个性质，FFT的分治策略便依然正确：

1、w[0]=w
=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(mo−1)∗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，则

A∗B≡1(mod xn)

如果我们已经求出在A在模xn2的逆元B0

即A∗B0≡1(mod xn2)

因为A∗B≡1(mod xn)

所以A∗B≡1(mod xn2)

那么就有B−B0≡0(mod xn2)

两边平方并移项

B2≡B0(2B−B0)(mod xn)

两边同乘A

B≡B0(2−B0A)(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
=quicksortmi(ff
,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
*ff[n-1]%mo);
}