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 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=;
    while (x) {
        if (x&) ans=ans*base%p;
        x>>=;
        base=base*base%p;
    }
    return ans;
}
void init()
{
    jc[]=; inv_j[]=quickpow(jc[],p-);
    for (int i=;i<=n;i++) 
     jc[i]=(LL)jc[i-]*i%p,inv_j[i]=quickpow(jc[i],p-);
    for (int i=;i<=n*8;i<<=)
     wn[i]=quickpow(,(p-)/(i<<));
}
void NTT(int n,int *a,int opt)
{
    for (int i=,j=;i<n;i++) {
        if (i>j) swap(a[i],a[j]);
        for (int l=n>>;(j^=l)<l;l>>=);
    }
    for (int i=;i<n;i<<=) {
        LL wn1=wn[i]; 
        for (int p1=i<<,j=;j<n;j+=p1) {
            LL w=;
            for (int k=;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==-) reverse(a+,a+n);
}
void inverse(int n,int *a,int *b,int *c)
{
    if (n==) b[]=quickpow(a[],p-);
    else {
        inverse((n+)>>,a,b,c);
        int k=;
        for (k=;k<=(n<<);k<<=);
        for (int i=;i<n;i++) c[i]=a[i];
        for (int i=n;i<k;i++) c[i]=;
        NTT(k,c,);
        NTT(k,b,);
        for (int i=;i<k;i++) {
            b[i]=(LL)(-(LL)c[i]*b[i]%p)*b[i]%p;
            if (b[i]<) b[i]+=p;
        }
        NTT(k,b,-);
        int inv=quickpow(k,p-);
        for (int i=;i<k;i++) b[i]=(LL)b[i]*inv%p;
        for (int i=n;i<k;i++) b[i]=;
    }
}
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    scanf("%d",&n); init();
    int n1=; 
    for (n1=;n1<=n*2;n1<<=);
    a[]=;
    for (int i=;i<=n;i++) a[i]=(LL)quickpow(,(LL)i*(i-)/)*inv_j[i]%p;
    inverse(n1,a,b,c);
    memset(c,,sizeof(c)); 
    for (int i=;i<=n;i++) c[i]=(LL)quickpow(,(LL)i*(i-)/)*inv_j[i-]%p;
    NTT(n1,b,); NTT(n1,c,);
    for (int i=;i<=n1;i++) b[i]=(LL)b[i]*c[i]%p;
    NTT(n1,b,-);
    LL inv=quickpow(n1,p-);
    for (int i=;i<=n1;i++) b[i]=(LL)b[i]*inv%p;
    printf("%d\n",(LL)b[n]*jc[n-]%p);
}