【BZOJ】3456: 城市规划 动态规划+多项式求逆

【题意】求n个点的带标号无向连通图个数 mod 1004535809。n<=130000。

【算法】动态规划+多项式求逆

【题解】设$g_n$表示n个点的无向图个数，那么显然

$$g_n=2^{frac{n(n-1)}{2}}$$

设$f_n$表示n个点的无向连通图个数，通过枚举1号点所属连通块大小很容易得到$g_n$的等式：

$$g_n=sum_{i=1}^{n}inom{n-1}{i-1}*f_i*g_{n-i}$$

特别的，$g_0=1$。

将组合数拆分一下，即可得到：

$$frac{g_n}{(n-1)!}=sum_{i=1}^{n}frac{f_i}{(i-1)!}*frac{g_{n-i}}{(n-i)!}$$

多项式求逆即可，注意下标从1开始需要强制F(0)=0，以及由于$g_0=1$所以右边额外+1。（这个式子恰好是多项式乘法的形式）

复杂度O(n log n)。

#include<cstdio>
#include<cstring>
#include<cctype>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<vector>
#include<algorithm>
#define ll long long
#define lowbit(x) x&-x
using namespace std;
int read(){
    char c;int s=0,t=1;
    while(!isdigit(c=getchar()))if(c=='-')t=-1;
    do{s=s*10+c-'0';}while(isdigit(c=getchar()));
    return s*t;
}
int min(int a,int b){return a<b?a:b;}
int max(int a,int b){return a<b?b:a;}
int ab(int x){return x>0?x:-x;}
//int MO(int x){return x>=MOD?x-MOD:x;}
//void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;}
/*------------------------------------------------------------*/
const int inf=0x3f3f3f3f,maxn=270010,MOD=1004535809;

int n,F[maxn],G[maxn],H[maxn];

int power(int x,int k){int ans=1;while(k){if(k&1)ans=1ll*ans*x%MOD;x=1ll*x*x%MOD;k>>=1;}return ans;}
int inv(int x){return power(x,MOD-2);}
int M(int x){return x>=MOD?x-MOD:x;}

void ntt(int *a,int n,int f){
    int k=0;
    for(int i=0;i<n;i++){
        if(i<k)swap(a[i],a[k]);
        for(int j=n>>1;(k^=j)<j;j>>=1);
    }
    for(int l=2;l<=n;l<<=1){
        int m=l>>1,wn=(~f?power(3,(MOD-1)/l):power(3,MOD-1-(MOD-1)/l));
        for(int *p=a;p!=a+n;p+=l){
            int w=1;
            for(int i=0;i<m;i++){
                int t=1ll*w*p[i+m]%MOD;
                p[i+m]=M(p[i]-t+MOD);
                p[i]=M(p[i]+t);
                w=1ll*w*wn%MOD;
            }
        }
    }
    if(f==-1){
        int o=inv(n);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*o%MOD;
    }
}
int h[maxn];
void pinv(int *f,int *g,int n){
    if(n==1){g[0]=inv(f[0]);return;}
    pinv(f,g,n>>1);n<<=1;
    for(int i=0;i<n/2;i++)h[i]=f[i];
    //for(int i=n/2;i<n;i++)h[i]=g[i]=0;
    ntt(h,n,1);ntt(g,n,1);
    for(int i=0;i<n;i++)g[i]=1ll*g[i]*(2-1ll*h[i]*g[i]%MOD+MOD)%MOD;
    ntt(g,n,-1);
    for(int i=n/2;i<n;i++)g[i]=0;//!
}
int fac[maxn],fav[maxn];
int main(){
    scanf("%d",&n);n++;
    fac[0]=1;for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%MOD;
    fav[n]=inv(fac[n]);for(int i=n;i>=1;i--)fav[i-1]=1ll*fav[i]*i%MOD;
    for(int i=1;i<n;i++)H[i]=1ll*power(2,1ll*i*(i-1)/2%(MOD-1))*fav[i-1]%MOD;
    for(int i=0;i<n;i++)G[i]=1ll*power(2,1ll*i*(i-1)/2%(MOD-1))*fav[i]%MOD;//n change
    int N=1;while(N<n+n)N<<=1;
    pinv(G,F,N>>1);//n>>1
    ntt(H,N,1);ntt(F,N,1);
    for(int i=0;i<N;i++)F[i]=1ll*H[i]*F[i]%MOD;
    ntt(F,N,-1);
    printf("%lld",1ll*F[n-1]*fac[n-2]%MOD);
    return 0;
}

注意：求逆元的时候传进去一倍即可，里面会再翻倍。逆元每次最后g数组后半部分要清零。n++之后对应的n要改变。

NTT最好预处理omega[]，不然NTT太多次会变得很慢，即：

#include<cstdio>
#include<cstring>
#include<cctype>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<vector>
#include<algorithm>
#define ll long long
#define lowbit(x) x&-x
using namespace std;
int read(){
    char c;int s=0,t=1;
    while(!isdigit(c=getchar()))if(c=='-')t=-1;
    do{s=s*10+c-'0';}while(isdigit(c=getchar()));
    return s*t;
}
int min(int a,int b){return a<b?a:b;}
int max(int a,int b){return a<b?b:a;}
int ab(int x){return x>0?x:-x;}
//int MO(int x){return x>=MOD?x-MOD:x;}
//void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;}
/*------------------------------------------------------------*/
const int inf=0x3f3f3f3f,maxn=270010,MOD=1004535809;

int n,F[maxn],G[maxn],H[maxn];

int power(int x,int k){int ans=1;while(k){if(k&1)ans=1ll*ans*x%MOD;x=1ll*x*x%MOD;k>>=1;}return ans;}
int inv(int x){return power(x,MOD-2);}
int M(int x){return x>=MOD?x-MOD:x;}
int o[maxn],oi[maxn];
void init(int n){
    int x=1,y=power(3,(MOD-1)/n);
    for(int i=0;i<=n;i++){
        o[i]=oi[n-i]=x;
        x=1ll*x*y%MOD;
    }
}
void ntt(int *a,int n,int *o,int f){
    int k=0;
    for(int i=0;i<n;i++){
        if(i<k)swap(a[i],a[k]);
        for(int j=n>>1;(k^=j)<j;j>>=1);
    }
    for(int l=2;l<=n;l<<=1){
        int m=l>>1;
        for(int *p=a;p!=a+n;p+=l){
            for(int i=0;i<m;i++){
                int t=1ll*o[n/l*i]*p[i+m]%MOD;
                p[i+m]=M(p[i]-t+MOD);
                p[i]=M(p[i]+t);
            }
        }
    }
    if(f==-1){
        int o=inv(n);
        for(int i=0;i<n;i++)a[i]=1ll*a[i]*o%MOD;
    }
}
int h[maxn];
void pinv(int *f,int *g,int n){
    if(n==1){g[0]=inv(f[0]);return;}
    pinv(f,g,n>>1);n<<=1;
    init(n);
    for(int i=0;i<n/2;i++)h[i]=f[i];
    //for(int i=n/2;i<n;i++)h[i]=g[i]=0;//?
    ntt(h,n,o,1);ntt(g,n,o,1);
    for(int i=0;i<n;i++)g[i]=1ll*g[i]*(2-1ll*h[i]*g[i]%MOD+MOD)%MOD;
    ntt(g,n,oi,-1);
    for(int i=n/2;i<n;i++)g[i]=0;//?
}
int fac[maxn],fav[maxn];
int main(){
    scanf("%d",&n);n++;
    fac[0]=1;for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%MOD;
    fav[n]=inv(fac[n]);for(int i=n;i>=1;i--)fav[i-1]=1ll*fav[i]*i%MOD;
    for(int i=1;i<n;i++)H[i]=1ll*power(2,1ll*i*(i-1)/2%(MOD-1))*fav[i-1]%MOD;
    for(int i=0;i<n;i++)G[i]=1ll*power(2,1ll*i*(i-1)/2%(MOD-1))*fav[i]%MOD;//!
    int N=1;while(N<n+n)N<<=1;
    pinv(G,F,N>>1);
    init(N);
    ntt(H,N,o,1);ntt(F,N,o,1);
    for(int i=0;i<N;i++)F[i]=1ll*H[i]*F[i]%MOD;
    ntt(F,N,oi,-1);
    printf("%lld",1ll*F[n-1]*fac[n-2]%MOD);
    return 0;
}

这道题有个O(n^2)的递推做法，设$h_n$表示n个点无向不连通图个数，那么$h_n=g_n-f_n$，枚举1所在连通块大小，有：

$$h_n=sum_{i=1}^{n-1}inom{n-1}{i-1}*f_i*h_{n-i}$$

这个不包含n自己，就可以递推了。

在这里记个无关紧要的笔记……要拆分$2^{k(n-k)}$，乘法不好拆分，必须转化为加法。

其实就是要把nk转化为有关n或k或n-k的加减法的形式，由初中常用套路$nk=frac{(n-k)^2-n^2-k^2}{2}$就可以了。