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

题目描述

传送门

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

题解

设$f[n]$表示n个点简单连通图个数（即1所属的连通块内有n个点）
$f[n]=2^{(n-1)*n/2}-\sum\limits_{i=0}^{n-1}{C_{n-1}^{i-1}*f[i]*2^{(n-i)*(n-i-1)/2}}$
这个怎么理解呢？对于每一条边来说，有连与不连两种状态，但是这样直接计算会包含很多不连通的图。假设1所属的连通块内有i个点，那么满足该条件的不合法的图的数量就是$C_{n-1}^{i-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}-\sum\limits_{i=0}^{n-1}{C_{n-1}^{i-1}*f[i]*2^{(n-i)*(n-i-1)/2}}$
$f[n]=2^{(n-1)*n/2}-\sum\limits_{i=0}^{n-1}{(n-1)!*{f[i]\over{(i-1)!}}*{2^{(n-i)*(n-i-1)/2}\over{(n-i)!}}}$
然后等式两边同时除去$(n-1)!$,得到
${f[n]\over{(n-1)!}}={2^{(n-1)*n/2}\over (n-1)!}-\sum\limits_{i=0}^{n-1}{{f[i]\over{(i-1)!}}*{2^{(n-i)*(n-i-1)/2}\over{(n-i)!}}}$
将含有$f$ 的项合并，得到
设$a[i]={f[i]\over{(i-1)!}}$,$b[i]={{2^{i*(i-1)/2}}\over{i!}}$ 其中$b[0]=1$
${2^{(n-1)*n/2}\over (n-1)!}=\sum\limits_{i=0}^n{a[i]*b[n-i]}$
设$c[i]={2^{(n-1)*n/2}\over (n-1)!}$,那么我们就得到了熟悉的卷积
$c[j]=\sum\limits_{i=0}^n{a[i]*b[n-i]}$
$C(x)=A(x)B(x)$
如果我们可以解出多项式$A(x)$的系数那么就可以求出$f[n]$，那么问题就迎刃而解了。
$A(x)\equiv{C(x)*B(x)^{-1}} \space mod \space x^{n+1}$
然后问题就变成了多项式求逆。

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

多项式的逆元

对于一个多项式$A(x)$，如果存在$B(x)$满足$deg_B\leq{deg_A}$并且

$$A(x)B(x)\equiv1 \space mod \space x^n$$
那么称 $B(x)$为 $A(x)$在 $mod \space x^n$意义下的逆元

求解过程

现在考虑如何求$A^{-1}(x)$,当$n=1$时，$A(x)\equiv c \space mod \space x$ ,$c$是常数，这样$A^{-1}(x)$ 就是$c^{-1}$
对于$n>1$的情况，设$B(x)=A^{-1}(x)$ ,由定义可知

$$A(x)B(x)\equiv1 \space mod \space x^n \space \space \space \space \space (1)$$
假设在 $mod \space x^{\lceil{n\over 2}\rceil}$意义下 $A(x)$的逆元是 $B'(x)$,并且我们已经求出，那么 $$A(x)B'(x)\equiv1 \space mod \space x^{\lceil{n\over 2}\rceil} \space \space \space \space \space (2)$$
再将 $(1)$放到 $mod \space x^{\lceil{n\over 2}\rceil}$意义下 $$A(x)B(x)\equiv1 \space mod \space x^{\lceil{n\over 2}\rceil} \space \space \space \space \space (3)$$
然后 $(2)-(3)$,得
$$B(x)-B(x)'\equiv 0 \space mod \space x^{\lceil{n\over 2}\rceil}$$
两边平方
$$B^2(x)-2B'(x)B(x)+B'^2(x)\equiv 0 \space mod \space x^n$$

为啥是对的？哈
因为多项式在$mod \space x^n$意义下为0 ，说明其0到$n-1$次项的系数都是0，那么平方后对于第$0\leq i \leq 2n-1$项，其系数$a_i$为$\sum\limits_{i=0}^i {a_i a_{j-i}}$,那么$i$,$j-i$中一定有一项的是属于$0\leq i < n$范围的，所以平方后在$mod \space x^{2n}$的意义下仍然是0
然后同时两边同乘$A(x)$,得到

$${B^2(x)A(x)}\equiv{2B'(x)B(x)A(x)-A(x)B'^2(x)} \space\space mod \space x^n$$
将 ${B(x)A(x)}\equiv 1 \space mod x^n$带入消去，得
$$B(x)\equiv {B'(x)*(2-B'(x)*A(x))} \space mod \space x^n$$
这一步运算可以用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 1004535809
using 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);
}