bzoj3456-城市规划
题意
求 \(n\) 个点的简单无向连通图个数。\(n\le 130000\) 。
分析
设 \(f(n)\) 为 \(n\) 个点的 带标号 简单无向连通图的个数,那么总的简单无向图个数 \(h(n)\) 就是
\[\begin{aligned}
h(n)&=\sum _{k=1}^n\sum _{x_1+x_2+\cdots +x_k=n}\binom n {x_1}f(x_1)\binom {n-x_1}{x_2}f(x_2)\cdots \binom {n-x_1-x_2\cdots x_{k-1}}{x_k}f(x_k) \\
&=n!\sum _{k=1}^n\prod _{\sum x=n}\frac{f(x_i)}{x_i!}
\end{aligned}
\]
即 \(n\) 个点的简单无向图个数就是所有连通块情况的和,其中 \(f\) 内部是带标号的,而连通块之间是无序的。显然,\(h(n)\) 的另一种计算方式是任选一些边连起来,即 \(h(n)=2^\binom n 2\) 。
显然这是一个生成函数的形式,有
\[\frac{h(n)}{n!}=\sum _{k=1}^n\prod _{\sum x=n}\frac{f(x_i)}{x_i!}
\]
令
\[H=\sum _{k=0}^\infty \frac{h(k)}{k!}x^k \\
F=\sum _{k=0}^\infty f(k)x^k
\]
那么就有
\[H=\sum _{k=1}^\infty \frac{F^k}{k!}
\]
即
\[H=\exp F-1
\]
我们不需要计算无限次数的生成函数,只需要计算 \(n\) 位就行了。因为 \(H\) 非常好求,所以我们用 \(H\) 来得到 \(F\) ,即
\[F=\ln (H+1)
\]
用计算 \(\ln\) 的方法 即可。复杂度为 \(O(n\log n)\) 。
代码
今天尝试找出一些好的写法,变成一个自己的多项式风格,然而几乎是失败了——封装一些东西导致它跑得很慢,差不多其他人的两倍时间,这是因为之前用了 vector
。后来发现复杂度没什么差别,而且 vector
在这种情况下可能还会慢一点。于是就改回了数组实现。
核心思想是用一个 P
指针实现递归。
#include<bits/stdc++.h>
using namespace std;
typedef long long giant;
const int maxn=1<<18;
const int q=479<<21|1;
const int toro=3;
typedef vector<int> vec;
inline int Plus(int x,int y) {return ((giant)x+(giant)y)%q;}
inline void Pe(int &x,int y) {x=Plus(x,y);}
inline int Multi(int x,int y) {return (giant)x*y%q;}
inline void Me(int &x,int y) {x=Multi(x,y);}
inline int Sub(int x,int y) {return Plus(x,q-y);}
inline void Se(int &x,int y) {x=Sub(x,y);}
template<typename T> inline int mi(int x,T y) {
int ret=1;
for (;y;y>>=1,Me(x,x)) if (y&1) Me(ret,x);
return ret;
}
inline int Inv(int x) {return x==1?1:Multi(q-q/x,Inv(q%x));}
//inline int Inv(int x) {return mi(x,q-2);}
inline int Div(int x,int y) {return Multi(x,Inv(y));}
int Wn[22][2],inv[maxn],fac[maxn],ifac[maxn],rev[maxn],n;
inline void getM(int n,int &M) {for (M=1;M<=(n<<1);M<<=1);}
struct P {
int a[maxn],M;
inline P (int m=0) {
memset(a,0,sizeof a);
give(m);
}
inline int& operator [] (int x) {return a[x];}
inline void set(int n) {getM(n,M);}
inline void give(int m) {M=m;}
inline void deri() {
for (int i=0;i<M-1;++i) a[i]=Multi(a[i+1],i+1);
a[M-1]=0;
}
inline void inte() {
for (int i=M-1;i>0;--i) a[i]=Multi(a[i-1],::inv[i]);
a[0]=0;
}
inline void ntt(int len,int op=0) {
for (int i=0,xj=__builtin_ctz(len);i<len;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(xj-1));
for (int i=0;i<len;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
for (int i=2,the=1;i<=len;i<<=1,++the)
for (int j=0,&wn=Wn[the][op];j<len;j+=i)
for (int k=j,w=1;k<j+i/2;++k,Me(w,wn)) {
const int u=a[k],v=Multi(a[k+i/2],w);
a[k]=Plus(u,v),a[k+i/2]=Sub(u,v);
}
if (op) for (int i=0,iv=::Inv(len);i<len;++i) Me(a[i],iv);
}
inline void operator *= (P &b) {
assert(M==b.M);
ntt(M),b.ntt(b.M);
for (int i=0;i<M;++i) Me(a[i],b[i]);
ntt(M,1),b.ntt(b.M,1);
}
P* inv(int m) {
P *ret;
if (m==2) {
ret=new P(m);
(*ret)[0]=::Inv(a[0]);
return ret;
}
const int hf=m>>1;
ret=inv(hf);
ret->give(m);
static P aux;
aux.give(m);
for (int i=0;i<m;++i) aux[i]=i<hf?a[i]:0;
aux.ntt(m),ret->ntt(m);
for (int i=0;i<m;++i) {
int &x=(*ret)[i];
Me(x,Sub(2,Multi(x,aux[i])));
}
ret->ntt(m,1);
for (int i=hf;i<m;++i) (*ret)[i]=0;
return ret;
}
P ln() {
static P pr,*iv;
(pr=*this).deri();
pr*=*(iv=inv(M));
pr.inte();
delete iv;
return pr;
}
} F,H;
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
#endif
scanf("%d",&n);
for (int i=1;i<22;++i) Wn[i][1]=Inv(Wn[i][0]=mi(toro,(q-1)>>i));
if (n<2) puts("1"),exit(0);
inv[1]=fac[0]=fac[1]=ifac[0]=ifac[1]=1;
for (int i=2;i<maxn;++i) {
inv[i]=Multi(q-q/i,inv[q%i]);
fac[i]=Multi(fac[i-1],i);
ifac[i]=Multi(ifac[i-1],inv[i]);
}
H.set(n);
H[0]=1;
for (int i=1;i<=n;++i) H[i]=Multi(mi(2,(giant)i*(i-1)/2),ifac[i]);
F=H.ln();
int ans=Multi(F[n],fac[n]);
printf("%d\n",ans);
return 0;
}