前言
这道题既是我NTT的第一份模板,也是第一道FFT相关的应用题。
题目描述
刚才说过, 阿狸的国家有n 个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通.
为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n 个点的简单(无重边无自环)无向连通图数目.(同构也算不同的)
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^21 + 1)即可.
分析
直接求好像很难求,怎么办呢?
先考虑n个点时无向图有多少种,设为G[N],则
G[n]=2C2n
。
得到了这个,我们可以推出n个点时联通无向图有多少种了,设为F[N],那么我们知道,F[n]可以由G[n]减去不合法的图。直观上的想法是一个合法的图和一些混乱的点拼在一起,即f[n-i]*g[i]*组合数,不过要考虑怎么排除重复。考虑把第一个点定住,即一定属于合法的那部分,这样就可以算出所有情况了。即
F[n]=G[n]−∑n−1i=1Ci−1n−1G[i]∗F[n−i]
.
这样得到了n^2的方法。哎我们发现后面那一块特别像卷积,试着化一化。
令
f[n]=∑n−1i=1Ci−1n−1G[i]∗F[n−i]
把组合数拆开
f[n]=∑n−1i=1(n−1)!(i−1)!(n−i)!G[i]∗F[n−i]
f[n]=(n−1)!∑n−1i=11(i−1)!(n−i)!G[i]∗F[n−i]
f[n]=(n−1)!∑n−1i=1G[i](i−1)!⋅F[n−i](n−i)!
拉开
n!
,设
h[n]=∑n−1i=1G[i](i−1)!⋅F[n−i](n−i)!
那么
h[n]
就是一个卷积了嘛。考虑上带模运算的FFT,但是对
h[n]
而言,前面的所有
F
都要先算出来,那么我们可以使用分治FFT。
在分治之前,我们先考虑一些东西
为了方便叙述,把G[n]/(i-1)!叫做G’[n],F也一样,F’。
边界条件。由于FFT快速处理的是形如
分治FFT
我们的目标是,在分治区间[l,r]中,先分治算出[l,mid]的所有F,再把这区间的F’的贡献加给[mid+1,r],再分治处理[mid+1,r]这个区间。
考虑给F’[l~r]重编号为a[1~r-l+1];而G’[]恰好用到的是1~r-l+1,直接把对应位赋值给b[1~r-l+1],那么我们只需要把a[1~mid-l+1]和b[1~r-l+1],FFT起来就好。再把位置对应回去,即F[l+i-1]+=a[i](l+i-1>mid),就好了。
注意一点,FFT中的数组长度要是a和b次数界较大值的两倍,不然逆DFT的时候会错。
这样时间复杂度是
O(nlog2n)
的。
NTT的实现
(中文不记得了···
其实跟FFT没什么区别,就是把单位复数根换成了模数的原根,精度是绝对保证的,根本没有实数运算。
指数与原根
给定一个模数p,一个与p互质的数a,(一般做题p是质数)。使得 ad%p=1 满足的最小的正整数d,称为a对模p的指数,若d等于phi(p),那么把a称作p的一个原根。
常用性质
我们知道
aphi(p)%p=1
,那么指数
d|phi(p)
;
对于质数p而言,它的任意原根g,都可以通过指数运算表示p的整个剩余系,当然0除外。即
{gx%p|x∈[1,p−1]}=[1,p−1]
寻找原根
一个数可能有多个原根,也可能没有。
原根存在性:一个数有原根的充要条件是,他形如:
1,2,4,p,2∗p,pn
,其中p是奇质数。
一般找原根都很快,反正一道题找到了就直接打个表就好了。
一般题目都是找质数的原根,那就假设找质数p的原根就好了。
找原根的方法:暴力枚举,合理判断。
1,首先暴力从2到p-1枚举;
2,验证是否原根。
而a是原根的条件是对于任意的
x∈[1,p−2],(ax%p)!=1
,而实际上有上面的性质,我们只要看看p-1的因数(除他自己)是不是全都满足上面的式子就好了,那么只用判断所有x=p-1/(p-1的质因数)即可,因为如果比这些数小的x0不满足,即
ax0
同余0,那么我们枚举的x里面肯定至少有一个是x0的倍数,那
ax
也同余0。
回到NTT
我们的原根g相当于什么呢?就是 w1mo−1 ,群的性质,折半引理和求和引理都是符合的。那我们的n次主单位根就很容易弄了嘛,我们要 w1n ,就是 wmo−1/nmo−1 ,这用消去引理可以推出来。当然前提是mo-1包含的2的次数要足够大,不然没法NTT。而其他部分是一样的。实现的时候有两种方法,一种是每次NTT前把会用到的w全部存起来;另外一种是每次做的时候算出主单位根,逆DFT的时候要弄成逆元。我习惯第二种。
代码
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
using namespace std;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
typedef long long ll;
typedef double db;
const int N=550005;
const int root=3;
const int mo=1004535809;
int g[N],f[N],a[N],b[N],t[N],fac[N],rev[N],w[20],halfw[20],rw[20],rhalfw[20];
int len,Log,n;
int ksm(int x,ll y)
{
int ret=1;
while (y)
{
if (y&1) ret=1ll*ret*x%mo;
y>>=1;
x=1ll*x*x%mo;
}
return ret;
}
void predo(int n)
{
//fac,rev,w,ksm,halfw
int i,siz;
fac[0]=1;
fo(i,1,n) fac[i]=1ll*fac[i-1]*i%mo;
rev[n]=ksm(fac[n],mo-2);
fd(i,n,1)
{
rev[i-1]=1ll*rev[i]*i%mo;
g[i]=ksm(2,(ll)i*(i-1)/2);
}
siz=1;
fo(i,0,n)
{
w[i]=ksm(root,(mo-1)/siz);
rw[i]=ksm(w[i],mo-2);
siz<<=1;
if (siz>n) break;
}
}
void DFT(int *a,int n,int sig)
{
int i,siz,half,tw,hw,W,u,v,j,k,tmp,ws;
fo(i,0,n-1)
{
int p=0;
for(tmp=i,ws=0;ws<Log;tmp/=2,ws++) p=(p<<1)+(tmp&1);
t[p]=a[i];
}
for(ws=1,siz=2;siz<=n;siz*=2,ws++)
{
half=siz/2;
if (sig==1) W=w[ws];else W=rw[ws];//其实可以把剩余系按原根弄出来的顺序排
tw=1;
fo(i,0,half-1)
{
for(j=i;j<=n;j+=siz)
{
k=j+half;
u=t[j];v=1ll*t[k]*tw%mo;
t[j]=(u+v)%mo;
t[k]=(u-v+mo)%mo;
}
tw=1ll*tw*W%mo;
}
}
fo(i,0,n-1) a[i]=t[i];
}
int max2(int n)
{
int ret=1;
for(;ret<n;ret<<=1);
return ret<<1;
}
int ci(int n){return (n==1)?0:ci(n/2)+1;}
void FFT(int l,int r,int rr)
{
int i,inv;
len=max2(rr-l+1);
// if (len==32) len=16;
Log=ci(len);
inv=ksm(len,mo-2);
// get the rem
fo(i,0,len-1) {a[i]=0;b[i]=0;}
fo(i,1,r-l+1)//必须保证位置正确,不然没法算
a[i]=1ll*f[i+l-1]*rev[i+l-2]%mo;
fo(i,1,rr-l)
b[i]=1ll*g[i]*rev[i]%mo;
DFT(a,len,1);
DFT(b,len,1);
fo(i,0,len-1) a[i]=1ll*a[i]*b[i]%mo;
DFT(a,len,-1);
fo(i,0,len-1) a[i]=1ll*a[i]*inv%mo;
}
void solve(int l,int r)
{
if (l==r)
{
f[l]=(g[l]-1ll*f[l]*fac[l-1]%mo+mo)%mo;
return;
}
int mid=(l+r)/2;
solve(l,mid);
FFT(l,mid,r);
int i;
fo(i,mid+1,r)
f[i]=(f[i]+a[i-l+1])%mo;
solve(mid+1,r);
}
int main()
{
freopen("3303.in","r",stdin);
freopen("3303.out","w",stdout);
scanf("%d",&n);
predo(130000*4);
solve(1,max2(n)/2);
printf("%d\n",f[n]);
}