[BZOJ3456] 城市规划 - 快速傅里叶变换 - 快速数论变换 - 卷积 - 多项式求逆

-太穷买不起权限号只能嘴巴+其他网站测评AC

----------以下部分为纸张蒟蒻zzt想到的部分，神犇可忽略不计----------

    啥子？要求代标号无向连通图？这不是思博题？不是一眼就转化为图总数-不连通的数目？蛤？不连通数目是啥？QAQ窝觉得是∑大小为i的连通图数目*大小为n-i的非连通图数目…（旁边：这不是一脸错误？泥算算小数据对吗昂）太弱被狂D。

----------以上部分就是纸张蒟蒻zzt想到的显而易见的性质，神犇可忽略不计----------

题解：

设g(i)=大小为i的图总数=2^(i*(i-1)/2)，f(i)=大小为i的连通图数目

枚举第一个点的连通块大小，则g(n)=∑C(n-1,j-1)f(j)g(n-j),j=1...n

显然我们要求的是f，那么我们怎么求f(n)呢

将g(i)的表达式带入原式，得2^(n*(n-1)/2)=∑f(j)*(n-1)!/(j-1)!/(n-j)!*2^((n-j)*(n-j-1)/2)

那么将两边同处以(n-1)!，发现式子对应生成函数的卷积，令A=∑2^(i*(i-1)/2)/(i-1)!*x^i，B=∑2^(i*(i-1)/2)/i!*x^i，C=∑f(i)*x^i

得A=BC 即C=A*(B^-1)

那么我们只要求出B在mod x^(n+1)意义下的逆元再用NTT与A相乘即可。

求逆元可以分块FFT，先求Bmod x^floor(n/2)意义下的逆元为S'，则Bmod x^(n+1)意义下的逆元S满足：

BS≡1 BS'≡1(mod x^floor(n/2)) 即B(S-S')≡0，显然B!=0(i=0...∞) 所以S-S'≡0(mod x^(floor(n/2))

即S^2-2SS'+S'^2≡0(mod x^(n+1))

同乘B降次，得S=2S'-BS'^2，后面的部分NTT即可

#include"bits/stdc++.h"using namespace std;typedef long long ll;const ll P=(479<<21)+1,g=3;const int N=262150,LG=20;ll w[LG];ll ksm(ll a,ll t,ll mod){ll re=1;while(t){if(t%2)re=re*a%mod;a=a*a%mod;t>>=1;}return re;}void build(int len){for(int i=1;(1<<i)<=len;i++)w[i]=ksm(g,(P-1)>>i,P);}template <class __t>void Rader(__t a[],int len){int i,j=len>>1,k;for(i=1;i<len-1;i++){if(i<j)swap(a[i],a[j]);for(k=len>>1;j>=k;j-=k,k>>=1);if(j<k)j+=k;}}void NTT(ll a[],int len,int DFT){int i,j,k,l,t;Rader(a,len);for(t=1;(1<<t)<=len;t++){l=1<<t;i=l>>1;for(j=0;j<len;j+=l){ll wn=1;for(k=j;k<j+i;k++){ll u=a[k],v=a[k+i]*wn%P;a[k]=(u+v)%P;a[k+i]=(u-v+P)%P;wn=wn*w[t]%P;}}}if(DFT==-1){for(i=1;i<len>>1;i++)swap(a[i],a[len-i]);ll mul=ksm(len,P-2,P);for(i=0;i<len;i++)a[i]=a[i]*mul%P;}}void Convol(ll a[],ll b[],ll in[],int len){NTT(a,len,1);NTT(b,len,1);for(int i=0;i<len;i++)in[i]=a[i]*b[i]%P;NTT(in,len,-1);}void SQR(ll a[],int len){NTT(a,len,1);for(int i=0;i<len;i++)a[i]=a[i]*a[i]%P;NTT(a,len,-1);}ll s[N],x[N],y[N];void inverse(ll a[],int deg){if(deg==1){s[0]=ksm(a[0],P-2,P);return;}int mid=(deg+1)>>1,i,L;inverse(a,mid);memcpy(x,s,sizeof(s));memset(y,0,sizeof(y));for(i=0;i<deg;i++)y[i]=a[i];for(L=1;L<=mid*2;L<<=1);SQR(x,L);for(;L<=deg*2;L<<=1);Convol(x,y,x,L);for(i=0;i<deg;i++)s[i]=(2*s[i]-x[i]+P)%P;}ll F[N],G[N],ans[N];ll fac[N>>1],invfac[N>>1],inv[N>>1];int pri[N>>1],flag[N>>1],cnt,n,l;void prepare(){scanf("%d",&n);int i,j;for(l=1;l<=(n+1)<<1;l<<=1);build(l);for(fac[1]=invfac[1]=inv[1]=1,i=2;i<=n;i++){fac[i]=fac[i-1]*i%P;if(flag[i]==0){pri[++cnt]=i;inv[i]=ksm(i,P-2,P);}for(j=1;pri[j]*i<=n;j++){inv[pri[j]*i]=inv[pri[j]]*inv[i]%P;if(i%pri[j]==0) break;}invfac[i]=invfac[i-1]*inv[i]%P;}F[1]=G[1]=G[0]=1;for(i=2;i<=n;i++){ll kuai=ksm(2,i-1,P);G[i]=G[i-1]*kuai%P*inv[i]%P;F[i]=F[i-1]*kuai%P*inv[i-1]%P;}}void work(){inverse(G,n+1);Convol(F,s,ans,l);printf("%lld\n",ans[n]*fac[n-1]%P);}int main(){prepare();work();return 0;}