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

时间:2022-06-23 18:42:54

-太穷买不起权限号只能嘴巴+其他网站测评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;
}