[HEOI2016]求和 sum

[HEOI2016]求和 sum

标签： NTT cdq分治 多项式求逆 第二类斯特林数

---

Description

求$$sum_{i=0}^nsum_{j=0}^i S(i,j)×2^j×(j!)$$
其中S(i,j)代表第二类斯特林数。

Solution

解法一

记Bell数(B(n)=sum_{i=0}^nS(n,i))
根据第二类斯特林数的组合意义,(B(i))代表把n个球放进任意个相同的盒子的方案数。
那么有$$B(n)=sum_{i=0}^{n-1} C(n-1,i) B(i)$$
讨论一下有多少个球与第一个球同盒即可。

记(B'(n)=sum_{i=0}^nS(n,i)×i!)
这个的组合意义是把n个求放进任意个不同盒子的方案数。
有$$B'(n)=sum_{i=0}^{n-1} C(n,i)B'(i)$$
讨论一下第一个盒子有多少个球

现在记(f(n)=sum_{i=0}^nS(n,i)×2^i×i!)
这个的组合意义是把n个球放进任意个不同且有两种颜色盒子的方案数。
显然$$f(n)=2sum_{i=0}^{n-1} C(n,i)f(i)$$
相当于每个盒子有两种不同的颜色选。

化简式子

[{f(n)over n!}=sum_{i=0}^{n-1}{f(i)over i!} {2 over (n-i)!} ]

设(g(i)={2 over i!}x^i),特别地,(g(0)=0).
这就是一个卷积的形式了。
然后cdq+NTT即可。

解法2

以上面的式子为基础做多项式求逆。
构造一个生成函数(F(x)=sum_{i=0}^{infty}{f(i) over i!}x^i)与(G(x)=sum _{i=1}^{infty}g(i))

那么根据上面的递推式子,有(F(x)=F(x)×G(x)+1)
解得(F(x)={1 over {1-G(x)}})

解法3

考虑斯特林数的容斥求法。
我们令(f(i)=sum_{j=i}^nS(j,i))
那么答案显然是(sum_{i=0}^nf(i)2^ii!)
接下来考虑怎么求(f(i))

[f(i)=sum_ {j=i}^nS(j,i)=sum_ {j=0}^nS(j,i) \ =sum_ {j=0}^n {frac 1 {i!}}sum_{k=0}^i (-1)^kC(i,k)(i-k)^j \ =sum_ {k=0} {(-1)^k over k!} {sum_{j=0}^n (i-k)^j over (i-k)! } ]

一遍NTT即可。

Code

解法一

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define REP(i,a,b) for(int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
	int sum=0,p=1;char ch=getchar();
	while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
	if(ch=='-')p=-1,ch=getchar();
	while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
	return sum*p;
}

const int maxn=4e5+20;
const int mod=998244353;

inline int power(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b & 1)ans=(ll)ans*a%mod;
		b>>=1;
		a=(ll)a*a%mod;
	}
	return ans;
}

int n,f[maxn],g[maxn],jc[maxn],jcn[maxn],inv[maxn],rev[maxn];

void init()
{
	n=read();
	jc[0]=jcn[0]=jc[1]=jcn[1]=inv[1]=1;
	REP(i,2,n)inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod,jc[i]=(ll)i*jc[i-1]%mod,jcn[i]=(ll)inv[i]*jcn[i-1]%mod;
	REP(i,1,n)g[i]=jcn[i]*2%mod;
}

int A[maxn],B[maxn];

inline void NTT(int *p,int N,int op)
{
	int n=1,l=0;while(n<N)n<<=1,l++;
	REP(i,1,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	REP(i,0,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
	for(int i=1;i<n;i<<=1)
	{
		int W=power(3,(mod-1)/(i<<1));
		for(int j=0;j<n;j+=i<<1)
		{
			int w=1;
			for(int k=j;k<i+j;k++,w=(ll)w*W%mod)
			{
				int x=p[k],y=(ll)p[k+i]*w%mod;
				p[k]=x+y;p[k+i]=x-y;
				if(p[k]>=mod)p[k]-=mod;
				if(p[k+i]<0)p[k+i]+=mod;
			}
		}
	}
	if(op==-1)
	{
		int inv=power(n,mod-2);
		REP(i,0,n-1)p[i]=(ll)p[i]*inv%mod;
		reverse(p+1,p+n);
	}
}

void solve(int l,int r)
{
	if(l==r)return;
	int mid=(l+r)>>1;
	solve(l,mid);
	REP(i,0,mid-l)A[i]=f[i+l];
	REP(i,0,r-l)B[i]=g[i];
	int N=1;
	while(N<=(mid-l)+(r-l)+1)N<<=1;
	NTT(A,N,1);NTT(B,N,1);
	REP(i,0,N-1)A[i]=(ll)A[i]*B[i]%mod;
	NTT(A,N,-1);
	REP(i,mid-l+1,r-l)f[i+l]=(A[i]+f[i+l])%mod;
	REP(i,0,N-1)A[i]=B[i]=0;
	solve(mid+1,r);
}

void doing()
{
	f[0]=1;
	solve(0,n);
	int ans=0;
	REP(i,0,n)ans=(ans+(ll)f[i]*jc[i])%mod;
	printf("%d
",ans);
}

int main()
{
	init();
	doing();
	return 0;
}

解法二

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define REP(i,a,b) for(int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
	int sum=0,p=1;char ch=getchar();
	while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
	if(ch=='-')p=-1,ch=getchar();
	while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
	return sum*p;
}

const int maxn=3e5+20;
const int mod=998244353;

int n,jc[maxn],jcn[maxn],inv[maxn],rev[maxn];

inline int power(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b & 1)ans=(ll)ans*a%mod;
		b>>=1;
		a=(ll)a*a%mod;
	}
	return ans;
}

inline void init()
{
	n=read();jc[0]=jcn[0]=jc[1]=jcn[1]=inv[1]=1;
	REP(i,2,n)inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod,jc[i]=(ll)i*jc[i-1]%mod,jcn[i]=(ll)inv[i]*jcn[i-1]%mod;
}

int g[maxn],ans[maxn],A[maxn],B[maxn],C[maxn];

inline void NTT(int *p,int N,int op)
{
	int n=1,l=0;while(n<N)n<<=1,l++;
	REP(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	REP(i,0,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
	for(int i=1;i<n;i<<=1)
	{
		int W=power(3,(mod-1)/(i<<1));
		for(int j=0;j<n;j+=i<<1)
		{
			int w=1;
			for(int k=j;k<i+j;k++,w=(ll)W*w%mod)
			{
				int x=p[k],y=(ll)w*p[k+i]%mod;
				p[k]=x+y;p[k+i]=x-y;
				if(p[k]>=mod)p[k]-=mod;
				if(p[k+i]<0)p[k+i]+=mod;
			}
		}
	}
	if(op==-1)
	{
		int inv=power(n,mod-2);
		REP(i,0,n-1)p[i]=(ll)p[i]*inv%mod;
		reverse(p+1,p+n);
	}
}
	
void Inv(int *p,int *q,int len)
{
	if(len==1){ q[0]=power(p[0],mod-2);return;}
	Inv(p,q,len>>1);
	REP(i,0,len-1)A[i]=p[i],B[i]=q[i];
	NTT(A,len<<1,1);NTT(B,len<<1,1);
	REP(i,0,(len<<1)-1)A[i]=(ll)A[i]*B[i]%mod*B[i]%mod;
	NTT(A,len<<1,-1);
	REP(i,0,len-1)q[i]=((2*q[i]-A[i])%mod+mod)%mod;//len不需要左移
	REP(i,0,(len<<1)-1)A[i]=B[i]=0;
}

inline void doing()
{
	int N=1,l=0;
	REP(i,1,n)g[i]=(-jcn[i]*2+2*mod)%mod;
	g[0]=1;
	while(N<=n)N<<=1,l++;
	Inv(g,ans,N);
	int Ans=0;
	REP(i,0,n)Ans=(Ans+(ll)ans[i]*jc[i])%mod;
	printf("%d
",Ans);
}

int main()
{
	init();
	doing();
	return 0;
}

解法三

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<map>
using namespace std;
#define ll long long
#define REP(i,a,b) for(int i=(a),_end_=(b);i<=_end_;i++)
#define DREP(i,a,b) for(int i=(a),_end_=(b);i>=_end_;i--)
#define EREP(i,a) for(int i=start[(a)];i;i=e[i].next)
inline int read()
{
	int sum=0,p=1;char ch=getchar();
	while(!(('0'<=ch && ch<='9') || ch=='-'))ch=getchar();
	if(ch=='-')p=-1,ch=getchar();
	while('0'<=ch && ch<='9')sum=sum*10+ch-48,ch=getchar();
	return sum*p;
}

const int maxn=4e5+20;
const int mod=998244353;

int n,inv[maxn],jc[maxn],jcn[maxn];

inline int power(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b & 1)ans=(ll)ans*a%mod;
		b>>=1;
		a=(ll)a*a%mod;
	}
	return ans;
}

int f[maxn],g[maxn],rev[maxn];

inline void NTT(int *p,int n,int opt)
{
	REP(i,0,n-1)if(i<rev[i])swap(p[i],p[rev[i]]);
	for(int i=1;i<n;i<<=1)
	{
		int W=power(3,(mod-1)/(i<<1));
		for(int j=0;j<n;j+=i<<1)
		{
			int w=1;
			for(int k=j;k<i+j;k++,w=(ll)w*W%mod)
			{
				int x=p[k],y=(ll)p[k+i]*w%mod;
				p[k]=x+y;
				p[k+i]=x-y;
				if(p[k]>=mod)p[k]-=mod;
				if(p[k+i]<0)p[k+i]+=mod;
			}
		}
	}
	if(opt==-1)
	{
		int inv=power(n,mod-2);
		REP(i,0,n-1)p[i]=(ll)p[i]*inv%mod;
		reverse(p+1,p+n);
	}
}

void init()
{
	n=read();inv[1]=jc[0]=jc[1]=jcn[0]=jcn[1]=1;
	REP(i,2,n+1)jc[i]=(ll)i*jc[i-1]%mod,inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod,jcn[i]=(ll)inv[i]*jcn[i-1]%mod;
	f[0]=g[0]=1;g[1]=n+1;
	REP(i,1,n)f[i]=(ll)(mod-inv[i])*f[i-1]%mod;
	REP(i,2,n)g[i]=(ll)(power(i,n+1)-1)*inv[i-1]%mod*jcn[i]%mod;
	int N=1,l=0;while(N<=2*n)N<<=1,l++;
	REP(i,1,N-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	NTT(f,N,1);
	NTT(g,N,1);
	REP(i,0,N-1)f[i]=(ll)f[i]*g[i]%mod;
	NTT(f,N,-1);
}

void doing()
{
	int ans=0,bin=1;
	REP(i,0,n)ans=(ans+(ll)f[i]*jc[i]%mod*bin)%mod,bin=(bin<<1)%mod;
	printf("%d
",ans);
}

int main()
{
	init();
	doing();
	return 0;
}