bzoj 3456 城市规划 —— 分治FFT / 多项式求逆 / 指数型生成函数(多项式求ln)
题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3456
首先考虑DP做法,正难则反,考虑所有情况减去不连通的情况;
而不连通的情况就是那个经典做法:选定一个划分点,枚举包含它的连通块,连通块以外的部分随便连(但不和连通块连通),合起来就是不连通的方案数;
设 \( f[i] \) 表示一共 \( i \) 个点时的连通方案数,\( g[i] \) 表示 \( i \) 个点随便连的方案数,即 \( g[i] = 2^{C_{i}^{2}} \),则:
\( f[i] = 2^{C_{i}^{2}} - \sum\limits_{j=1}^{i-1} C_{i-1}^{j-1} * f_{j} * g_{i-j} \)
只要令 \( g[0] = 0 \),\( j \) 就可以枚举到 \( i \);
把 \( C_{i-1}^{j-1} \) 拆开,对应分配到各处,得到:
\( f[i] = 2^{C_{i}^{2}} - (i-1)!*\sum\limits_{j=1}^{i}\frac{f_{j}}{(j-1)!} * \frac{g_{i-j}}{(i-j)!} \)
所以把 \( g[i] \) 的定义改成 \( g[i] = \frac{2^{C_{i}^{2}}}{i!} \)
于是 \( f[i] = 2^{C_{i}^{2}} - (i-1)!*\sum\limits_{j=1}^{i}\frac{f_{j}}{(j-1)!} * g_{i-j} \)
然后就可以分治FFT了;
很容易写错的地方是那个 \( 2^{C_{i}^{2}} \),因为是指数,所以应该对 \( mod-1 \) 取模,而不是 \( mod \) !!
而且除以 \( 2 \) 也不能预处理 \( inv2 \) 了,因为 \( mod-1 \) 不是质数!!
代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=(1<<18),mod=1004535809; int n,f[xn],g[xn],jc[xn],jcn[xn],a[xn],b[xn],rev[xn],inv2,in[xn]; int rd() { int ret=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return f?ret:-ret; } ll pw(ll a,int b,int p=mod) { ll ret=1; for(;b;b>>=1,a=(a*a)%p)if(b&1)ret=(ret*a)%p; return ret; } int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} int C(int x){return ((ll)x*(x-1)/2)%(mod-1);}//mod-1!!! //not inv2 -- !pri void init() { jc[0]=1; for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod; jcn[n]=pw(jc[n],mod-2); for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod; g[0]=0; for(int i=1;i<=n;i++)g[i]=(ll)pw(2,C(i))*jcn[i]%mod; } void ntt(int *a,int tp,int lim) { for(int i=0;i<lim;i++) if(i<rev[i])swap(a[i],a[rev[i]]); for(int mid=1;mid<lim;mid<<=1) { int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len); for(int j=0;j<lim;j+=len) for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod) { int x=a[j+k],y=(ll)w*a[j+mid+k]%mod; a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y); } } if(tp==1)return; int inv=pw(lim,mod-2); for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod; } void work(int l,int r) { if(l==r){f[l]=upt((ll)pw(2,C(l))-(ll)jc[l-1]*f[l]%mod); return;} int len=r-l+1,mid=((l+r)>>1); work(l,mid); int lim=1,L=0; while(lim<=len)lim<<=1,L++;// for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(L-1))); for(int i=0;i<lim;i++)a[i]=b[i]=0; for(int i=l;i<=mid;i++)a[i-l]=(ll)f[i]*jcn[i-1]%mod;//jcn!! for(int i=0;i<=r-l;i++)b[i]=g[i];// ntt(a,1,lim); ntt(b,1,lim); for(int i=0;i<lim;i++)a[i]=(ll)a[i]*b[i]%mod; ntt(a,-1,lim); for(int i=mid+1;i<=r;i++)f[i]=upt(f[i]+a[i-l]); work(mid+1,r); } int main() { n=rd(); init(); work(1,n); printf("%d\n",f[n]); return 0; }
也可以用多项式求逆做:
设 \( g[i] = 2^{C_{i}^{2}} \),这里 \( g[0] = 1 \),\( f[i] \) 定义同上;
得到 \( g[n] = \sum\limits_{i=1}^{n} C_{n-1}^{i-1} * f[i] * g[n-i] \)
拆 \( C_{n-1}^{i-1} \),得到
\( \frac{g[n]}{(n-1)!} = \sum\limits_{i=1}^{n} \frac{f[i]}{(i-1)!} * \frac{g[n-i]}{(n-i)!} \)
设 \( A[i] = \frac{g[i]}{(i-1)!} \) , \( B[i] = \frac{f[i]}{(i-1)!} \) , \( C[i] = \frac{g[i]}{i!} \)
注意这里 \( A[0] = 0 \) , \( B[0] = 0 \) , \( C[0] = 1 \)
所以 \( A(x) = B(x) * C(x) \)
即 \( B(x) = A(x) * C^{-1}(x) ( mod x^{n+1} ) \)
多项式求逆即可;
别忘了外层的乘法还要处理一下 \( rev \) 数组。
代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=(1<<18),mod=1004535809; int n,a[xn],b[xn],c[xn],t[xn],rev[xn],jc[xn],jcn[xn]; int rd() { int ret=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return f?ret:-ret; } ll pw(ll a,int b) { ll ret=1; for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod; return ret; } int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} int C(int x){return ((ll)x*(x-1)/2)%(mod-1);} void init() { jc[0]=1; for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod; jcn[n]=pw(jc[n],mod-2); for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod; } void ntt(int *a,int tp,int lim) { for(int i=0;i<lim;i++) if(i<rev[i])swap(a[i],a[rev[i]]); for(int mid=1;mid<lim;mid<<=1) { int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len); for(int j=0;j<lim;j+=len) for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod) { int x=a[j+k],y=(ll)w*a[j+mid+k]%mod; a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y); } } if(tp==1)return; int inv=pw(lim,mod-2); for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod; } void inv(int *a,int *b,int n) { if(n==1){b[0]=pw(a[0],mod-2); return;} inv(a,b,(n+1)>>1); int lim=1,l=0; while(lim<n+n)lim<<=1,l++; for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1))); for(int i=0;i<n;i++)t[i]=a[i]; for(int i=n;i<lim;i++)t[i]=0; ntt(t,1,lim); ntt(b,1,lim); for(int i=0;i<lim;i++)b[i]=upt((2-(ll)t[i]*b[i])%mod*b[i]%mod); ntt(b,-1,lim); for(int i=n;i<lim;i++)b[i]=0; } int main() { n=rd(); init(); int lim=1,l=0; while(lim<=n)lim<<=1,l++; for(int i=1;i<lim;i++)a[i]=(ll)pw(2,C(i))*jcn[i-1]%mod; for(int i=1;i<lim;i++)c[i]=(ll)pw(2,C(i))*jcn[i]%mod; c[0]=1; inv(c,b,n+1); for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1)));//!!! ntt(a,1,lim); ntt(b,1,lim); for(int i=0;i<lim;i++)b[i]=(ll)a[i]*b[i]%mod; ntt(b,-1,lim); printf("%lld\n",(ll)b[n]*jc[n-1]%mod); return 0; }
也可以用指数型生成函数,具体做法可以看这篇博客:https://blog.csdn.net/wzq_qwq/article/details/48435621
关于为什么 \( G(x) = e^{F(x)} \)
因为从组合意义上来看,任意无向图实际上可以分为:由0个连通图组成(没有点),由1个连通图组成,由2个连通图组成(这2个之间不连通),由3个连通图组成...
所以 \( G(x) = 1 + \frac{F(x)}{1!} + \frac{F(x)^{2}}{2!} + \frac{F(x)^{3}}{3!} + ... \)
即 \( G(x) = e^{F(x)} \)
而 \( G(x) \) 很好构造,求 \( F(x) = lnG(x) \)
求 \( ln \) 的方法可以看这两篇博客:
https://blog.csdn.net/litble/article/details/81749788
https://blog.csdn.net/ezoixx118/article/details/81235586
直接用上公式...这题原来是多项式求 \( ln \) 的裸题;
别忘了最后乘上 \( n! \)
代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; typedef long long ll; int const xn=(1<<18),mod=1004535809; int n,g[xn],dg[xn],ig[xn],t[xn],rev[xn],jc[xn],jcn[xn]; int rd() { int ret=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();} while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar(); return f?ret:-ret; } ll pw(ll a,ll b) { ll ret=1; b=b%(mod-1); for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod; return ret; } int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} void init() { jc[0]=1; for(int i=1;i<=n;i++)jc[i]=(ll)jc[i-1]*i%mod; jcn[n]=pw(jc[n],mod-2); for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod; } void ntt(int *a,int tp,int lim) { for(int i=0;i<lim;i++) if(i<rev[i])swap(a[i],a[rev[i]]); for(int mid=1;mid<lim;mid<<=1) { int len=(mid<<1),wn=pw(3,tp==1?(mod-1)/len:(mod-1)-(mod-1)/len); for(int j=0;j<lim;j+=len) for(int k=0,w=1;k<mid;k++,w=(ll)w*wn%mod) { int x=a[j+k],y=(ll)w*a[j+mid+k]%mod; a[j+k]=upt(x+y); a[j+mid+k]=upt(x-y); } } if(tp==1)return; int inv=pw(lim,mod-2); for(int i=0;i<lim;i++)a[i]=(ll)a[i]*inv%mod; } void inv(int *a,int *b,int n) { if(n==1){b[0]=pw(a[0],mod-2); return;} inv(a,b,(n+1)>>1); int lim=1,l=0; while(lim<n+n)lim<<=1,l++; for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1))); for(int i=0;i<n;i++)t[i]=a[i]; for(int i=n;i<lim;i++)t[i]=0; ntt(t,1,lim); ntt(b,1,lim); for(int i=0;i<lim;i++)b[i]=upt((2-(ll)t[i]*b[i])%mod*b[i]%mod); ntt(b,-1,lim); for(int i=n;i<lim;i++)b[i]=0; } int main() { n=rd(); init(); for(int i=0;i<=n;i++) { if(i<2)g[i]=1; g[i]=(ll)pw(2,(ll)i*(i-1)/2)*jcn[i]%mod; } for(int i=1;i<=n;i++)dg[i-1]=(ll)i*g[i]%mod; dg[n]=0;// inv(g,ig,n+1); int lim=1,l=0; while(lim<n+n)lim<<=1,l++; for(int i=0;i<lim;i++)rev[i]=((rev[i>>1]>>1)|((i&1)<<(l-1))); ntt(ig,1,lim); ntt(dg,1,lim); for(int i=0;i<lim;i++)dg[i]=(ll)ig[i]*dg[i]%mod; ntt(dg,-1,lim); printf("%lld\n",(ll)dg[n-1]*pw(n,mod-2)%mod*jc[n]%mod); return 0; }