Tibbar的后花园(生成函数，多项式exp)
Tibbar的后花园(多项式exp)
这篇文章仅限于没有入门生成函数的蒟蒻读，dalao勿喷
题目的数据范围是\(n< 2^{20}\)
对于联通块求
题目给出的限制，其实就是对于每一个联通块
1.不存在一个点度数\(\ge 3\)
2.不存在一个长度为\(3\)的倍数的环
可以看到，一个大小为\(n\)的联通块，合法方案只有为一条链，或者一个长度不为3的倍数的环
设\(f_n\)为大小为\(n\)的联通块方案数，则
\(f_n=\left\{\begin{aligned}1 && n\leq 2 \\ \frac{n!}{2}+[n \mod 3\ne 0](\frac{(n-1)!}{2}) && n\ge 3\end{aligned}\right.\)
即链和环的个数
对于图求
如果采用普通的分治+NTT枚举每个联通块大小，复杂度为\(O(n\log ^2n)\)，显然无法通过这个题
考虑构造生成函数快速求解
可以看到，如果排列\(n\)个点，每个点生成\(m\)个大小为\(a_i(\sum a_i=n)\)的联通块，那么答案是
\(\sum \frac{n!}{m!} \Pi \frac{f_{a_i}}{a_i!}\)
即除去联通块内的排列，联通块之间的排列
令\(g_n=\sum \frac{f_i}{i!}x^i\)
那么答案可以转化为
\[n
\]
就是\(i\)次累乘(叠加了\(i\)个联通块)之后，第\(n\)项的系数再乘上少掉的阶乘
然后由于
\(\sum \frac{x^i}{i!}=e^x\)
所以就是求出\(n![x^n]e^{g(x)}\)
多项式exp我跑得血慢啊
#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cctype>
#include<vector>
#include<cassert>
using namespace std;
//#pragma GCC optimize(2)
#pragma GCC optimize(3)
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)
#define pb push_back
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
char IO;
int rd(){
int s=0;
int f=0;
while(!isdigit(IO=getchar())) if(IO=='-') f=1;
do s=(s<<1)+(s<<3)+(IO^'0');
while(isdigit(IO=getchar()));
return f?-s:s;
}
const int N=(1<<21)+10,P=1004535809;
int n;
ll qpow(ll x,ll k){
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
int Mod_Inv[N];
int rev[N],w[N],iw[N],tmp[N];
int PreMake(int n){
reg int R=1,cc=-1;
while(R<n) R<<=1,cc++;
rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
return R;
}
void NTT(int n,int *a,int f){
rep(i,1,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
if(f==1) {
for(reg int i=1;i<n;i<<=1) {
reg int len=(1<<20)/i;
for(int j=0,t=0;j<i;++j,t+=len) tmp[j]=w[t];
for(reg int l=0;l<n;l+=i*2) {
for(reg int j=l;j<l+i;j++) {
reg int t=1ll*a[j+i]*tmp[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
} else {
for(reg int i=1;i<n;i<<=1) {
reg int len=(1<<20)/i;
for(int j=0,t=0;j<i;++j,t+=len) tmp[j]=iw[t];
for(reg int l=0;l<n;l+=i*2) {
for(reg int j=l;j<l+i;j++) {
reg int t=1ll*a[j+i]*tmp[j-l]%P;
a[j+i]=a[j]-t,Mod2(a[j+i]);
a[j]+=t,Mod1(a[j]);
}
}
}
int base=Mod_Inv[n];
rep(i,0,n-1) a[i]=1ll*a[i]*base%P;
}
}
namespace Inv{
//F(x)G(x)=1
// G(x)^2+H(x)^2-2G(x)H(x)=0
// G(x)+H(x)^2*F(x)-2H(x)=0
// G(x)=2H(x)-H(x)^2*F(x)
int A[N];
void Inv(int *a,int *b,int n) {
int len=1;
b[0]=qpow(a[0],P-2);
for(len=2;;len<<=1) {
int R=PreMake(len<<1);
rep(i,0,len-1) A[i]=a[i];
rep(i,len,R-1) A[i]=b[i]=0;
NTT(R,A,1),NTT(R,b,1);
rep(i,0,R-1) b[i]=(2-1ll*b[i]*A[i]%P+P)*b[i]%P;
NTT(R,b,-1);
rep(i,len,R) b[i]=0;
if(len>=n) break;
}
rep(i,n,len) b[i]=0;
}
}
namespace Ln{
// G(x)=ln F(x);
// G'(x) = F'(x)/F(x)
int c[N],d[N];
void Ln(int *a,int *b,int n) {
rep(i,0,n) b[i]=c[i]=0;
rep(i,1,n-1) c[i-1]=1ll*i*a[i]%P;
Inv::Inv(a,b,n-1);
int R=PreMake(n<<1);
rep(i,n,R) c[i]=b[i]=0;
NTT(R,c,1),NTT(R,b,1);
rep(i,0,R-1) b[i]=1ll*b[i]*c[i]%P;
NTT(R,b,-1);
rep(i,n,R) b[i]=0;
drep(i,n-1,1) b[i]=1ll*b[i-1]*Mod_Inv[i]%P;
b[0]=0;
}
}
namespace Exp{
int c[N];
void Exp(int *a,int *b,int n){
b[0]=1;
int len;
for(len=2;;len<<=1) {
rep(i,0,len) c[i]=0;
Ln::Ln(b,c,len);
int R=PreMake(len<<1);
rep(i,0,len-1) {
c[i]=(!i)-c[i]+a[i];
Mod2(c[i]);
}
rep(i,len,R) c[i]=0;
NTT(R,b,1),NTT(R,c,1);
rep(i,0,R-1) b[i]=1ll*b[i]*c[i]%P;
NTT(R,b,-1);
rep(i,len,R) b[i]=0;
if(len>=n) break;
}
rep(i,n,len) b[i]=0;
}
}
int Fac[N],g[N],f[N];
int main(){
n=1<<21;
w[0]=1,w[1]=qpow(3,(P-1)/n);
rep(i,2,n-1) w[i]=1ll*w[i-1]*w[1]%P;
iw[0]=1,iw[1]=qpow((P+1)/3,(P-1)/n);
rep(i,2,n-1) iw[i]=1ll*iw[i-1]*iw[1]%P;
Fac[0]=1;
rep(i,1,N-1) Fac[i]=1ll*Fac[i-1]*i%P;
Mod_Inv[1]=1;
rep(i,2,N-1) Mod_Inv[i]=1ll*(P-P/i)*Mod_Inv[P%i]%P;
n=rd();
f[1]=f[2]=1;
rep(i,3,n) {
f[i]=1ll*Fac[i]*(P+1)/2%P;
if(i%3) f[i]=(f[i]+1ll*Fac[i-1]*(P+1)/2)%P;
} // 预处理转移系数，注意还要额外除去阶乘
for(reg int i=1,t=1;i<=n;t=1ll*t*Mod_Inv[++i]%P) f[i]=1ll*f[i]*t%P;
Exp::Exp(f,g,n+1);
printf("%lld\n",1ll*g[n]*Fac[n]%P);
}
