BZOJ4555 TJOI2016&HEOI2016 求和

Problem

BZOJ

Solution

\[\sum_{i=0}^n\sum_{j=0}^i S(i,j)\times 2^j\times j!\tag{1} \]

我们令

\[F(i)=\sum_{j=0}^i S(i,j)\times 2^j\times j!\tag{2} \]

思考一下\(F(i)\)的含义，就是把 \(i\) 个球放入若干个盒子的方案数，其中盒子带标号，且有两种状态。

那么我们枚举最后一个盒子放多少个球即可得到递推式

\[F(i)=2\sum_{j=1}^i \binom i j F(i-j)=2i!\sum_{j=1}^i \frac {1} {j!}\times \frac {F(i-j)} {(i-j)!}\tag{3} \]

\[\frac {F(i)} {i!}=\sum_{j=1}^i \frac 2 {j!}\times \frac {F(i-j)} {(i-j)!}\tag{4} \]

我们令\(A(i)=\frac {F(i)}{i!},B(i)=\frac 2 {i!}\)

\[A(x)=\sum_{j=1}^i B(j)A(i-j)\tag{5} \]

不难发现这是一个分治FFT的模板，我们考虑转多项式求逆解决

---

分治FFT转多项式求逆的推导：

由于\(B(0)=0\)

\[\sum_{j=1}^iB(j)A(i-j)=\sum_{j=1}^iB(j)A(i-j)+B(0)A(i)\tag{6} \]

\[=\sum_{j=0}^i B(j)A(i-j)=(A*B)(i)\tag{7} \]

\(A(i)=(A*B)(i)\)

然而这样是不对的，注意\(A(0)\)是有初始值的。\(A(0)\not =0\)，而\((A*B)(0)=0\)

因此\(A(i)=(A*B)(i)+A(0)\)

\[A(x)=\frac {1} {A(0)-B(x)}\tag{8} \]

多项式求逆即可。

时间复杂度优化至\(O(n\log n)\)

Code

#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=300010,mod=998244353,G=3;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
    x=0;int f=0;char ch=getchar();
    while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
    if(ch=='-') f=1,ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    if(f) x=-x;
}
int n,ans,fac[maxn],f[maxn],g[maxn],h[maxn],rev[maxn];
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
int power(int x,int y)
{
	int res=1;
	for(;y;y>>=1,x=(ll)x*x%mod)
	  if(y&1)
	    res=(ll)res*x%mod;
	return res;
}
void NTT(int *a,int N,int f)
{
	for(int i=1;i<N;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=1;i<N;i<<=1)
	{
		int gn=power(G,(mod-1)/(i<<1));
		for(int j=0;j<N;j+=(i<<1))
		{
			int g=1,x,y;
			for(int k=0;k<i;++k,g=(ll)g*gn%mod)
			{
				x=a[j+k];y=(ll)g*a[j+k+i]%mod;
				a[j+k]=pls(x,y);a[j+k+i]=dec(x,y);
			}
		}
	}
	if(f==-1)
	{
		int inv=power(N,mod-2);reverse(a+1,a+N);
		for(int i=0;i<N;i++) a[i]=(ll)a[i]*inv%mod;
	}
}
void work(int n)
{
	if(n==1){f[0]=1;return ;}
	work((n+1)>>1);
	int N,l;
	for(N=1,l=0;N<n+n;N<<=1) ++l;
	for(int i=1;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0;i<n;i++) h[i]=g[i];
	for(int i=n;i<N;i++) h[i]=0;
	NTT(f,N,1);NTT(h,N,1);
	for(int i=0;i<N;i++) f[i]=dec(pls(f[i],f[i]),(ll)f[i]*f[i]%mod*h[i]%mod);
	NTT(f,N,-1);
	for(int i=n;i<N;i++) f[i]=0;
}
int main()
{
	read(n);
	fac[0]=1;g[0]=1;
	for(int i=1;i<=n;i++)
	{
		fac[i]=(ll)fac[i-1]*i%mod;
		g[i]=2*power(fac[i],mod-2);
		if(g[i]>=mod) g[i]-=mod;
		g[i]=dec(0,g[i]);
	}
	work(n+1);
	for(int i=0;i<=n;i++) ans=pls(ans,(ll)f[i]*fac[i]%mod);
	printf("%d\n",ans);
	return 0;
}