BZOJ3456: 城市规划
Description
刚刚解决完电力网络的问题, 阿狸又被领导的任务给难住了.
刚才说过, 阿狸的国家有n个城市, 现在国家需要在某些城市对之间建立一些贸易路线, 使得整个国家的任意两个城市都直接或间接的连通. 为了省钱, 每两个城市之间最多只能有一条直接的贸易路径. 对于两个建立路线的方案, 如果存在一个城市对, 在两个方案中是否建立路线不一样, 那么这两个方案就是不同的, 否则就是相同的. 现在你需要求出一共有多少不同的方案.
好了, 这就是困扰阿狸的问题. 换句话说, 你需要求出n个点的简单(无重边无自环)无向连通图数目.
由于这个数字可能非常大, 你只需要输出方案数mod 1004535809(479 * 2 ^ 21 + 1)即可.
Input
仅一行一个整数n(<=130000)
Output
仅一行一个整数, 为方案数 mod 1004535809.
Sample Input
3
Sample Output
4
HINT
对于 100%的数据, n <= 130000
O(N^2)的做法：
设f[i]表示n个点的无向连通图数目，考虑容斥原理计算n个点的无向非连通图数目。
枚举包含第1个点的连通分量大小j，不难得出f[i]=2C(i,2)-∑f[j]*C(i-1,j-1)*2C(i-j,2)。
然后把转移方程拆开：f[i]=2C(i,2)-(i-1)!*∑(f[j]*(j-1)!)*((i-j)!*2C(i-j,2))
设A[i]=f[j]*(j-1)!，B[i]=(i-j)!*2C(i-j,2)，那么式子就是A与B的卷积了。
然后我们就可以分治+NTT辣！
具体细节可见code：
#include<cstdio>
#include<cctype>
#include<queue>
#include<cstring>
#include<algorithm>
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define dwn(i,s,t) for(int i=s;i>=t;i--)
#define ren for(int i=first[x];i;i=next[i])
using namespace std;
inline int read() {
int x=0,f=1;char c=getchar();
for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
for(;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const int maxn=132000;
const int G=3;
const int p=1004535809;
typedef long long ll;
ll pow(ll n,ll m,ll mod=p) {
ll ans=1;
for(;m;m>>=1,(n*=n)%=mod) if(m&1) (ans*=n)%=p;
return ans;
}
ll wn[20];
void NTT(ll* A,int len,int tp) {
int j=len>>1,c=0;
rep(i,1,len-2) {
if(i<j) swap(A[i],A[j]);int k=len>>1;
while(j>=k) j-=k,k>>=1;j+=k;
}
for(int i=2;i<=len;i<<=1) {
c++;
for(int j=0;j<len;j+=i) {
ll w=1;
for(int k=j;k<j+(i>>1);k++) {
ll u=A[k],t=w*A[k+(i>>1)]%p;
A[k]=(u+t)%p;A[k+(i>>1)]=(u-t+p)%p;
w=(w*wn[c])%p;
}
}
}
if(tp<0) {
ll inv=pow(len,p-2);
rep(i,1,len/2-1) swap(A[i],A[len-i]);
rep(i,0,len-1) A[i]=(A[i]*inv)%p;
}
}
ll xp[maxn],inv[maxn];
ll f[maxn],T[maxn],A[maxn],B[maxn];
void solve(int l,int r) {
if(l==r) return;
int mid=l+r>>1,len=1;solve(l,mid);
while(len<=(max(mid-l+1,r-mid)<<1)) len<<=1;
rep(i,0,len-1) A[i]=B[i]=0;
rep(i,l,mid) A[i-l]=f[i]*inv[i-1]%p;
rep(i,1,r-l) B[i]=inv[i]*T[i]%p;
NTT(A,len,1);NTT(B,len,1);
rep(i,0,len-1) A[i]=(A[i]*B[i])%p;
NTT(A,len,-1);
rep(i,mid+1,r) f[i]=((f[i]-xp[i-1]*A[i-l])%p+p)%p;
solve(mid+1,r);
}
int main() {
xp[0]=inv[0]=1;int n=read();
rep(i,1,19) wn[i]=pow(G,(p-1)/(1<<i));
rep(i,1,n) xp[i]=(xp[i-1]*i)%p,inv[i]=pow(xp[i],p-2),f[i]=T[i]=pow(2,(ll)i*(i-1)/2);
solve(1,n);printf("%lld\n",f[n]);
return 0;
}
(UPD)从Po姐那里学来了O(NlogN)的多项式逆元做法。
本题题解：http://blog.csdn.net/popoqqq/article/details/46049331
多项式逆元：http://picks.logdown.com/posts/189620-the-inverse-element-of-polynomial%20%E8%B7%AApicks%E5%A4%A7%E6%AF%92%E7%98%A4
#include<cstdio>
#include<cctype>
#include<queue>
#include<cstring>
#include<algorithm>
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define dwn(i,s,t) for(int i=s;i>=t;i--)
#define ren for(int i=first[x];i;i=next[i])
using namespace std;
inline int read() {
int x=0,f=1;char c=getchar();
for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
for(;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
typedef long long ll;
const int p=1004535809;
const int G=3;
const int maxn=270000;
ll wn[20];
ll pow(ll n,ll m,ll mod=p) {
ll ans=1;
for(;m;m>>=1,(n*=n)%=mod) if(m&1) (ans*=n)%=mod;
return ans;
}
void NTT(ll* A,int len,int tp) {
int j=len>>1,c=0;
rep(i,1,len-2) {
if(i<j) swap(A[i],A[j]);int k=len>>1;
while(j>=k) j-=k,k>>=1;j+=k;
}
for(int i=2;i<=len;i<<=1) {
c++;
for(int j=0;j<len;j+=i) {
ll w=1;
for(int k=j;k<j+(i>>1);k++) {
ll u=A[k],v=w*A[k+(i>>1)]%p;
A[k]=(u+v)%p;A[k+(i>>1)]=(u-v+p)%p;
w=(w*wn[c])%p;
}
}
}
if(tp<0) {
ll inv=pow(len,p-2);
rep(i,1,len/2-1) swap(A[i],A[len-i]);
rep(i,0,len-1) (A[i]*=inv)%=p;
}
}
ll tmp[maxn];
void getinv(ll* A,ll* B,int n) {
if(n==1) {B[0]=pow(A[0],p-2);return;}
getinv(A,B,n>>1);
rep(i,0,n-1) tmp[i]=A[i],tmp[i+n]=0;
NTT(B,n<<1,1);NTT(tmp,n<<1,1);
rep(i,0,(n<<1)-1) tmp[i]=(2-tmp[i]*B[i]%p+p)%p;
rep(i,0,(n<<1)-1) (B[i]*=tmp[i])%=p;
NTT(B,n<<1,-1);
rep(i,n,(n<<1)-1) B[i]=0;
}
ll A[maxn],B[maxn],C[maxn],B2[maxn];
ll xp[maxn],invxp[maxn],xp2[maxn];
int main() {
rep(i,0,19) wn[i]=pow(G,(p-1)/(1<<i));
int n=read(),len=1;while(len<=(n<<1))len<<=1;
xp[0]=invxp[0]=xp2[0]=1;
rep(i,1,n) xp2[i]=pow(2,(ll)i*(i-1)/2),xp[i]=(xp[i-1]*i)%p,invxp[i]=pow(xp[i],p-2);
rep(i,0,n) B[i]=xp2[i]*invxp[i]%p;
rep(i,1,n) C[i]=xp2[i]*invxp[i-1]%p;
getinv(B,B2,len>>1);
NTT(B2,len,1);NTT(C,len,1);
rep(i,0,len-1) A[i]=(B2[i]*C[i])%p;
NTT(A,len,-1);
printf("%lld\n",A[n]*xp[n-1]%p);
return 0;
}
