bzoj 3456 城市规划——分治FFT / 多项式求逆 / 多项式求ln
题目：https://www.lydsy.com/JudgeOnline/problem.php?id=3456
分治FFT：
设 dp[ i ] 表示 i 个点时连通的方案数。
考虑算补集：连通的方案数 == 随便连方案数 - 不连通方案数
不连通方案数就和很久之前做过的“地震后的幻想乡”一样，枚举一个连通的点集，其中需要一直包含一个“划分点”保证不重复；其余部分随便连。注意还有从 i 个点里选 j 个点作为连通点集的那个组合数。
\( dp[i]=2^{C^{2}_{i}} - \sum\limits^{i-1}_{j=1} dp[j]*C^{j-1}_{i-1}*2^{C^{2}_{i-j}} \)
\( dp[i]=2^{C^{2}_{i}} - (i-1)!\sum\limits^{i-1}_{j=1} ( dp[j]*\frac{1}{(j-1)!} )( 2^{C^{2}_{i-j}}*\frac{1}{(i-j)!} ) \)
就可以分治FFT啦！
一开始还记得指数是对 mod-1 取模，后来就忘了，调了好久……注意 mod-1 不是质数，且是偶数，所以2没有逆元，算 C 的时候直接 /2 就行了。
在 L==R 的地方把 f [ ] 弄好。虽然也可以给 f 赋上 \( C^{2}_{i} \) 的初值，然后卷积的时候每次 -=\( (i-1)!*a[j] \)，就不用在 L==R 的时候特殊判断了，但那样可能因为总要乘 \( (i-1)! \)，会变慢。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=130005,M=N<<1,mod=1004535809,m2=mod-1;//N<<1 int n,len,r[M],f[M],g[M],a[M],b[M],jc[N],jcn[N],inv2; int rdn() { int ret=0;bool fx=1;char ch=getchar(); while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();} while(ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar(); return fx?ret:-ret; } void upd(int &x){x>=mod?x-=mod:0;} int pw(int x,int k,int md=mod) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%md;x=(ll)x*x%md;k>>=1;}return ret;} int C(int n)///////%mod-1//// {return (ll)n*(n-1)/2%m2;}// /2 void init() { ////// inv2=pw(2,m2-2,m2);//////m2!!!! m2 is not prime!!!!! jc[1]=1;for(int i=2;i<n;i++)jc[i]=(ll)jc[i-1]*i%mod; jcn[n-1]=pw(jc[n-1],mod-2); for(int i=n-2;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod;//>=0 for(int i=1;i<n;i++)g[i]=(ll)pw(2,C(i))*jcn[i]%mod; } void ntt(int *a,bool fx) { for(int i=0;i<len;i++) if(i<r[i])swap(a[i],a[r[i]]); for(int R=2;R<=len;R<<=1) { int Wn=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R ); for(int i=0,m=R>>1;i<len;i+=R) for(int j=0,w=1;j<m;j++,w=(ll)w*Wn%mod) { int x=a[i+j], y=(ll)w*a[i+m+j]%mod; a[i+j]=x+y; upd(a[i+j]); a[i+m+j]=x+mod-y; upd(a[i+m+j]); } } if(!fx)return; int inv=pw(len,mod-2); for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%mod; } void solve(int L,int R) { if(L==R){f[L]=(pw(2,C(L))-(ll)jc[L-1]*f[L])%mod+mod;upd(f[L]);return;} int mid=L+R>>1; solve(L,mid); int d=R-L,i,j; for(len=1;len<d;len<<=1); for(i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0); for(i=0,j=L;j<=mid;i++,j++)a[i]=(ll)f[j]*jcn[j-1]%mod;//jcn for(;i<len;i++)a[i]=0; for(i=0,j=R-L;i<=j;i++)b[i]=g[i+1]; for(;i<len;i++)b[i]=0; ntt(a,0); ntt(b,0); for(i=0;i<len;i++)a[i]=(ll)a[i]*b[i]%mod; ntt(a,1); for(int i=mid+1,j=i-L-1;i<=R;i++,j++)f[i]+=a[j],upd(f[i]); solve(mid+1,R); } int main() { n=rdn(); init(); solve(1,n); printf("%d\n",f[n]); return 0; }
多项式求逆：
http://blog.miskcoo.com/2015/05/bzoj-3456
应该只有G(x)有0次项。在阶乘的地方看，0!应该等于1才对（预处理组合数时就是这样，想一想，\( C^{n}_{n} \)就是\( \frac{n!}{n!*0!} \)）。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=1e5+5,M=N<<2,mod=1004535809; int n,a[M],c[M],cn[M],A[M],jcn[M],len,r[M]; void upd(int &x){x>=mod?x-=mod:0;} int C(int n){return (ll)n*(n-1)/2%(mod-1);} int pw(int x,int k) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} void ntt(int *a,bool fx) { for(int i=0;i<len;i++) if(i<r[i])swap(a[i],a[r[i]]); for(int R=2;R<=len;R<<=1) { int Wn=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R ); for(int i=0,m=R>>1;i<len;i+=R) for(int j=0,w=1;j<m;j++,w=(ll)w*Wn%mod) { int x=a[i+j], y=(ll)w*a[i+m+j]%mod; a[i+j]=x+y; upd(a[i+j]); a[i+m+j]=x+mod-y; upd(a[i+m+j]); } } if(!fx)return ; int inv=pw( len,mod-2 ); for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%mod; } void getinv(int n,int *a,int *b) { if(n==1){b[0]=1;return;} getinv(n+1>>1,a,b); int i; for(len=1;len<n<<1;len<<=1); for(i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0); for(i=0;i<n;i++)A[i]=a[i]; for(;i<len;i++)A[i]=0; ntt(A,0); ntt(b,0); for(i=0;i<len;i++)b[i]=(2-(ll)A[i]*b[i]%mod)*b[i]%mod+mod,upd(b[i]); ntt(b,1); for(int i=n;i<len;i++)b[i]=0; } int main() { scanf("%d",&n); for(int i=1;i<=n;i++)a[i]=pw(2,C(i)); jcn[0]=1;for(int i=1;i<=n;i++)jcn[i]=(ll)jcn[i-1]*i%mod; jcn[n]=pw(jcn[n],mod-2);for(int i=n-1;i>=0;i--)jcn[i]=(ll)jcn[i+1]*(i+1)%mod; c[0]=1; for(int i=1;i<=n;i++)c[i]=(ll)a[i]*jcn[i]%mod; for(int i=1;i<=n;i++)a[i]=(ll)a[i]*jcn[i-1]%mod; getinv(n+1,c,cn); for(len=1;len<=n<<1;len<<=1); for(int i=0;i<len;i++)r[i]=(r[i>>1]>>1)+((i&1)?len>>1:0); ntt(a,0); ntt(cn,0); for(int i=0;i<len;i++)a[i]=(ll)a[i]*cn[i]%mod; ntt(a,1); printf("%lld\n",(ll)a[n]*pw(jcn[n-1],mod-2)%mod); return 0; }
多项式求 ln ：
关于指数型生成函数：https://max.book118.com/html/2015/1111/29175835.shtm
想求两种物品一共取 i 个的排列方案，设取第一种物品的方案是 f[ ] ，第二种物品的方案是 g[ ] ，那么就是想求
\( \sum\limits_{j=0}^{i}C_{i}^{j}f[i]*g[i-j] \)
令 F(x) 是 f[ ] 的指数型生成函数、G(x) 是 g[ ] 的生成函数的话，要求的就恰好是 F(x) * G(x) 的第 i 次项系数乘上 i! 。
关于本题：https://www.cnblogs.com/ivorysi/p/9756961.html
至于为什么式子 \( G(x)=\sum\limits_{i=1}^{\infty}\frac{F(i)}{i!} \) 里要除以 i 的阶乘，可以从 n=2 并希望有两个连通块的情况考虑：
令 \( g[ i ] = 2^{C_{n}^{2}} \) 表示 i 个点的无向图个数，这个算出来是有标号的；令 f[ i ] 表示 i 个点有标号无向连通图个数。
\( F^{2}(x) = \sum\limits_{i=0}^{n}\frac{ \sum\limits_{j=0}^{i}C_{i}^{j}f[j]*f[i-j] }{i!}x^i \)
直接这样的话，就有 \( g[2]=C_{2}^{0}*f[0]*f[2]+C_{2}^{1}*f[1]*f[1]+C_{2}^{2}f[2]*f[0] \)
可以发现，因为有了那个 C 的系数，所以即使是卷积的时候只会算一遍的 j == i-j 的情况，也会有重复；即那个博客里说的本质一样的情况。
关于求ln：http://www.cnblogs.com/zhoushuyu/p/8763215.html
#include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const int N=(1<<19)+5,mod=1004535809; int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;} int pw(int x,int k) {int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;} int n,len,r[N],f[N],g[N],A[N],B[N]; int jc[N],jcn[N],inv[N]; void init(int n) { jc[0]=jcn[0]=inv[0]=1; jc[1]=jcn[1]=inv[1]=1; for(int i=2;i<=n;i++) { jc[i]=(ll)jc[i-1]*i%mod; inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod; jcn[i]=(ll)jcn[i-1]*inv[i]%mod; } } void ntt(int *a,bool fx,bool deb=0) { for(int i=0;i<len;i++) if(i<r[i])swap(a[i],a[r[i]]); for(int R=2;R<=len;R<<=1) { int wn=pw( 3,fx?(mod-1)-(mod-1)/R:(mod-1)/R ); for(int i=0,m=R>>1;i<len;i+=R) for(int j=0,w=1;j<m;j++,w=(ll)w*wn%mod) { int x=a[i+j],y=(ll)w*a[i+m+j]%mod; a[i+j]=upt(x+y); a[i+m+j]=upt(x-y); } } if(!fx)return; int inv=pw(len,mod-2); for(int i=0;i<len;i++)a[i]=(ll)a[i]*inv%mod; } void get_dao(int n,int *a,int *b) { for(int i=1;i<n;i++)b[i-1]=(ll)a[i]*i%mod; b[n]=b[n-1]=0; } void get_jf(int n,int *a) { for(int i=n-1;i;i--)a[i]=(ll)a[i-1]*inv[i]%mod;//i-- a[0]=0; } void get_inv(int n,int *a,int *b) { if(n==1){b[0]=pw(a[0],mod-2);return;} get_inv(n>>1,a,b); len=n<<1; for(int i=0,j=len>>1;i<len;i++) r[i]=(r[i>>1]>>1)+((i&1)?j:0); for(int i=0;i<n;i++)A[i]=a[i],A[i+n]=0; ntt(A,0); ntt(b,0); for(int i=0;i<len;i++) b[i]=(ll)b[i]*(2-(ll)A[i]*b[i]%mod+mod)%mod; ntt(b,1); for(int i=n;i<len;i++)b[i]=0; } void get_ln(int n,int *a,int *b) { get_inv(n,a,B); get_dao(n,a,A);//get_inv will use A len=n<<1; for(int i=n;i<len;i++)A[i]=0;// for(int i=0,j=len>>1;i<len;i++) r[i]=(r[i>>1]>>1)+((i&1)?j:0); ntt(A,0,1); ntt(B,0,1); for(int i=0;i<len;i++)b[i]=(ll)A[i]*B[i]%mod; ntt(b,1); get_jf(n,b); for(int i=n;i<len;i++)b[i]=0;// } int main() { scanf("%d",&n);int l;for(l=1;l<=n;l<<=1);//l>n init(l); g[0]=1;/// for(int i=1;i<=n;i++) g[i]=(ll)pw(2,(ll)i*(i-1)/2%(mod-1))*jcn[i]%mod; get_ln(l,g,f); printf("%lld\n",(ll)f[n]*jc[n]%mod); return 0; }
