bzoj 3456: 城市规划 （NTT+多项式求逆）

题目描述

传送门

题目大意：求n个点简单无向连通图数，其中任意点之间可以随意连边，不存在重边和自环。

题解

设f[n]表示n个点简单连通图个数（即1所属的连通块内有n个点）
f[n]=2(n−1)∗n/2−∑i=0n−1Ci−1n−1∗f[i]∗2(n−i)∗(n−i−1)/2
这个怎么理解呢？对于每一条边来说，有连与不连两种状态，但是这样直接计算会包含很多不连通的图。假设1所属的连通块内有i个点，那么满足该条件的不合法的图的数量就是Ci−1n−1∗f[i]∗2(n−i)∗(n−i−1)/2,首先从n-1个点中选出i-1个点，然后这i-1个点和1构成合法的连通图的方案数是f[i]，剩下的n-i个点之间随便连。

我们将上面的式子进行化简和变形：
f[n]=2(n−1)∗n/2−∑i=0n−1Ci−1n−1∗f[i]∗2(n−i)∗(n−i−1)/2
f[n]=2(n−1)∗n/2−∑i=0n−1(n−1)!∗f[i](i−1)!∗2(n−i)∗(n−i−1)/2(n−i)!
然后等式两边同时除去(n−1)!,得到
f[n](n−1)!=2(n−1)∗n/2(n−1)!−∑i=0n−1f[i](i−1)!∗2(n−i)∗(n−i−1)/2(n−i)!
将含有f 的项合并，得到
设a[i]=f[i](i−1)!,b[i]=2i∗(i−1)/2i! 其中b[0]=1
2(n−1)∗n/2(n−1)!=∑i=0na[i]∗b[n−i]
设c[i]=2(n−1)∗n/2(n−1)!,那么我们就得到了熟悉的卷积
c[j]=∑i=0na[i]∗b[n−i]
C(x)=A(x)B(x)
如果我们可以解出多项式A(x)的系数那么就可以求出f[n]，那么问题就迎刃而解了。
A(x)≡C(x)∗B(x)−1 mod xn+1
然后问题就变成了多项式求逆。

PS: 对于一个多项式A(x),称其最高项的次数为这个多项式的度，记作degA

多项式的逆元

对于一个多项式A(x)，如果存在B(x)满足degB≤degA并且

A(x)B(x)≡1 mod xn

那么称B(x)为A(x)在mod xn意义下的逆元

求解过程

现在考虑如何求A−1(x),当n=1时，A(x)≡c mod x ,c是常数，这样A−1(x) 就是c−1
对于n>1的情况，设B(x)=A−1(x) ,由定义可知

A(x)B(x)≡1 mod xn     (1)

假设在mod x⌈n2⌉意义下A(x)的逆元是B′(x),并且我们已经求出，那么

A(x)B′(x)≡1 mod x⌈n2⌉     (2)

再将(1)放到mod x⌈n2⌉意义下

A(x)B(x)≡1 mod x⌈n2⌉     (3)

然后(2)−(3),得

B(x)−B(x)′≡0 mod x⌈n2⌉

两边平方

B2(x)−2B′(x)B(x)+B′2(x)≡0 mod xn

为啥是对的？哈
因为多项式在mod xn意义下为0 ，说明其0到n−1次项的系数都是0，那么平方后对于第0≤i≤2n−1项，其系数ai为∑i=0iaiaj−i,那么i,j−i中一定有一项的是属于0≤i<n范围的，所以平方后在mod x2n的意义下仍然是0
然后同时两边同乘A(x),得到

B2(x)A(x)≡2B′(x)B(x)A(x)−A(x)B′2(x)  mod xn

将B(x)A(x)≡1 modxn带入消去，得

B(x)≡B′(x)∗(2−B′(x)∗A(x)) mod xn

这一步运算可以用NTT进行加速，所以最后的时间复杂度为O(nlogn)

回到这道题，我们只需要用多项式求逆求出B−1(x),然后用NTT加入B−1(x)C(x)
得到最终的A(x),答案就是a[n]∗(n−1)!

代码

#include<iostream>#include<cstdio>#include<algorithm>#include<cstring>#include<cmath>#define N 5000003#define LL long long#define p 1004535809using namespace std;int n,m;int a[N],b[N],c[N],jc[N],inv_j[N],wn[N];LL quickpow(LL num,LL x){    LL base=num%p; LL ans=1;    while (x) {        if (x&1) ans=ans*base%p;        x>>=1;        base=base*base%p;    }    return ans;}void init(){    jc[0]=1; inv_j[0]=quickpow(jc[0],p-2);    for (int i=1;i<=n;i++)      jc[i]=(LL)jc[i-1]*i%p,inv_j[i]=quickpow(jc[i],p-2);    for (int i=1;i<=n*8;i<<=1)     wn[i]=quickpow(3,(p-1)/(i<<1));}void NTT(int n,int *a,int opt){    for (int i=0,j=0;i<n;i++) {        if (i>j) swap(a[i],a[j]);        for (int l=n>>1;(j^=l)<l;l>>=1);    }    for (int i=1;i<n;i<<=1) {        LL wn1=wn[i];         for (int p1=i<<1,j=0;j<n;j+=p1) {            LL w=1;            for (int k=0;k<i;k++,w=(LL)w*wn1%p) {                int x=a[j+k]; int y=(LL)a[j+k+i]*w%p;                a[j+k]=(x+y)%p; a[j+k+i]=(x-y+p)%p;            }        }    }    if (opt==-1) reverse(a+1,a+n);}void inverse(int n,int *a,int *b,int *c){    if (n==1) b[0]=quickpow(a[0],p-2);    else {        inverse((n+1)>>1,a,b,c);        int k=0;        for (k=1;k<=(n<<1);k<<=1);        for (int i=0;i<n;i++) c[i]=a[i];        for (int i=n;i<k;i++) c[i]=0;        NTT(k,c,1);        NTT(k,b,1);        for (int i=0;i<k;i++) {            b[i]=(LL)(2-(LL)c[i]*b[i]%p)*b[i]%p;            if (b[i]<0) b[i]+=p;        }        NTT(k,b,-1);        int inv=quickpow(k,p-2);        for (int i=0;i<k;i++) b[i]=(LL)b[i]*inv%p;        for (int i=n;i<k;i++) b[i]=0;    }}int main(){    freopen("a.in","r",stdin);    freopen("my.out","w",stdout);    scanf("%d",&n); init();    int n1=0;     for (n1=1;n1<=n*2;n1<<=1);    a[0]=1;    for (int i=1;i<=n;i++) a[i]=(LL)quickpow(2,(LL)i*(i-1)/2)*inv_j[i]%p;    inverse(n1,a,b,c);    memset(c,0,sizeof(c));     for (int i=1;i<=n;i++) c[i]=(LL)quickpow(2,(LL)i*(i-1)/2)*inv_j[i-1]%p;    NTT(n1,b,1); NTT(n1,c,1);    for (int i=0;i<=n1;i++) b[i]=(LL)b[i]*c[i]%p;    NTT(n1,b,-1);    LL inv=quickpow(n1,p-2);    for (int i=0;i<=n1;i++) b[i]=(LL)b[i]*inv%p;    printf("%d\n",(LL)b[n]*jc[n-1]%p);}