[NTT][JZOJ3303]【集训队互测2013】城市规划

前言

这道题既是我NTT的第一份模板，也是第一道FFT相关的应用题。

题目描述

刚才说过, 阿狸的国家有n 个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通.
为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n 个点的简单(无重边无自环)无向连通图数目.(同构也算不同的)
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^21 + 1)即可.

分析

直接求好像很难求，怎么办呢？
先考虑n个点时无向图有多少种，设为G[N]，则
$G[n]=2^{C_n^2}$。
得到了这个，我们可以推出n个点时联通无向图有多少种了，设为F[N]，那么我们知道，F[n]可以由G[n]减去不合法的图。直观上的想法是一个合法的图和一些混乱的点拼在一起，即f[n-i]*g[i]*组合数，不过要考虑怎么排除重复。考虑把第一个点定住，即一定属于合法的那部分，这样就可以算出所有情况了。即
$F[n]=G[n]-\sum_{i=1}^{n-1}C_{n-1}^{i-1}G[i]*F[n-i]$.
这样得到了n^2的方法。哎我们发现后面那一块特别像卷积，试着化一化。
令$f[n]=\sum_{i=1}^{n-1}C_{n-1}^{i-1}G[i]*F[n-i]$
把组合数拆开
$f[n]=\sum_{i=1}^{n-1}\frac{(n-1)!}{(i-1)!(n-i)!}G[i]*F[n-i]$
$f[n]=(n-1)!\sum_{i=1}^{n-1}\frac{1}{(i-1)!(n-i)!}G[i]*F[n-i]$
$f[n]=(n-1)!\sum_{i=1}^{n-1}\frac{G[i]}{(i-1)!}·\frac{F[n-i]}{(n-i)!}$
拉开$n!$，设$h[n]=\sum_{i=1}^{n-1}\frac{G[i]}{(i-1)!}·\frac{F[n-i]}{(n-i)!}$
那么$h[n]$就是一个卷积了嘛。考虑上带模运算的FFT，但是对$h[n]$而言，前面的所有$F$都要先算出来，那么我们可以使用分治FFT。
在分治之前，我们先考虑一些东西
为了方便叙述，把G[n]/(i-1)!叫做G’[n]，F也一样，F’。
边界条件。由于FFT快速处理的是形如$C[n]=\sum_{i=0}^{n}a[i]*b[n-i]$的标准形式，那么我们得考虑上面式子的i=0和i=n时的情况。可以发现对于i=0，我们只要令G’[0]=0即可；对于i=n，我们只要不把$F'[n]$赋值给搞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(nlog^2n)$的。

NTT的实现

(中文不记得了···
其实跟FFT没什么区别，就是把单位复数根换成了模数的原根，精度是绝对保证的，根本没有实数运算。

指数与原根

给定一个模数p，一个与p互质的数a，（一般做题p是质数）。使得$a^d\%p=1$满足的最小的正整数d，称为a对模p的指数，若d等于phi(p)，那么把a称作p的一个原根。

常用性质

我们知道$a^{phi(p)}\%p=1$，那么指数$d|phi(p)$；
对于质数p而言，它的任意原根g，都可以通过指数运算表示p的整个剩余系，当然0除外。即$\{g^x\%p|x\in[1,p-1]\}=[1,p-1]$

寻找原根

一个数可能有多个原根，也可能没有。
原根存在性：一个数有原根的充要条件是，他形如：$1,2,4,p,2*p,p^n$，其中p是奇质数。
一般找原根都很快，反正一道题找到了就直接打个表就好了。
一般题目都是找质数的原根，那就假设找质数p的原根就好了。
找原根的方法：暴力枚举，合理判断。
1，首先暴力从2到p-1枚举；
2，验证是否原根。
而a是原根的条件是对于任意的$x\in[1,p-2],(a^x\%p) !=1$，而实际上有上面的性质，我们只要看看p-1的因数（除他自己）是不是全都满足上面的式子就好了，那么只用判断所有x=p-1/(p-1的质因数)即可，因为如果比这些数小的x0不满足，即$a^{x0}$同余0，那么我们枚举的x里面肯定至少有一个是x0的倍数，那$a^x$也同余0。

回到NTT

我们的原根g相当于什么呢？就是$w_{mo-1}^1$，群的性质，折半引理和求和引理都是符合的。那我们的n次主单位根就很容易弄了嘛，我们要$w_n^1$，就是$w_{mo-1}^{mo-1/n}$，这用消去引理可以推出来。当然前提是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]);
}