「Luogu4091」[HEOI2016/TJOI2016]求和

「Luogu4091」[HEOI2016/TJOI2016]求和

problem

Solution

我们知道当(i<j)时，(S(i,j)=0)，所以

[sum_{i=0}^nsum_{j=0}^iS(i,j) imes 2^j imes j! ]

[=sum_{j=0}^n2^j imes j!sum_{i=0}^nS(i,j) ]

将第二类斯特林数的展开式代入得

[sum_{j=0}^n2^j imes j!sum_{i=0}^nfrac{1}{j!}sum_{k=0}^j(-1)^k imesegin{pmatrix}j\kend{pmatrix}(j-k)^i ]

[=sum_{j=0}^n2^j imes j!sum_{i=0}^nsum_{k=0}^jfrac{(-1)^k}{k!} imesfrac{(j-k)^i}{(j-k)!} ]

[=sum_{j=0}^n2^j imes j!sum_{k=0}^jfrac{(-1)^k}{k!} imesfrac{sum_{i=0}^n(j-k)^i}{(j-k)!} ]

[=sum_{j=0}^n2^j imes j!sum_{k=0}^jfrac{(-1)^k}{k!} imesfrac{(j-k)^{n+1}-1}{(j-k-1)(j-k)!} ]

后面那玩意看起来很像两个函数的卷积

我们令(f(x)=frac{(-1)^x}{x!},g(x)=frac{x^{n+1}-1}{(x-1)x!})

代入得

[sum_{j=0}^n2^j imes j!sum_{k=0}^jf(k) imes g(j-k) ]

[=sum_{j=0}^n2^j imes j!(f*g)(j) ]

然后就可以用(NTT)来预处理出((f*g))了

Code

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <iostream>
#define inv(x) ((fastpow((x),mod-2)))
using namespace std;
typedef long long ll;

template <typename T>void read(T &t)
{
	t=0;int f=0;char c=getchar();
	while(!isdigit(c)){f|=c=='-';c=getchar();}
	while(isdigit(c)){t=t*10+c-'0';c=getchar();}
	if(f)t=-t;
}

const ll mod=998244353,gg=3,ig=332748118;
const int maxn=100000+5;
int n;
ll f[maxn<<3],g[maxn<<3];
int len;
int rev[maxn<<3];

ll fastpow(ll a,ll b)
{
	ll re=1,base=a;
	while(b)
	{
		if(b&1)
			re=re*base%mod;
		base=base*base%mod;
		b>>=1;
	}
	return re;
}

void NTT(ll *f,int type)
{
    for(register int i=0;i<len;++i)
        if(i<rev[i])
            swap(f[i],f[rev[i]]);
    for(register int p=2;p<=len;p<<=1)
    {
        int length=p>>1;
        ll unr=fastpow(type==1?gg:ig,(mod-1)/p);
        for(register int l=0;l<len;l+=p)
        {
            register ll w=1;
            for(register int i=l;i<l+length;++i,w=w*unr%mod)
            {
                ll tt=f[i+length]*w%mod;
                f[i+length]=(f[i]-tt+mod)%mod;
                f[i]=(f[i]+tt)%mod;
            }
        }
    }
	if(!type)
	{
		ll ilen=inv(len);
		for(register int i=0;i<len;++i)
			f[i]=f[i]*ilen%mod;
	}
}

int main()
{
	read(n);
	f[0]=g[0]=1;
	g[1]=n+1,f[1]=mod-1;
	for(register ll i=2,fac=2;i<=n;++i,fac=fac*i%mod)
	{
		f[i]=f[i-1]*(mod-1)%mod*inv(i)%mod;
		g[i]=(fastpow(i,n+1)-1+mod)%mod*inv(i-1)%mod*inv(fac)%mod;
	}
	for(len=1;len<=n+n;len<<=1);
	for(register int i=0;i<len;++i)
	    rev[i]=((rev[i>>1]>>1)|((i&1)?len>>1:0));
	NTT(f,1);
	NTT(g,1);
	for(register int i=0;i<len;++i)
		f[i]=f[i]*g[i]%mod;
	NTT(f,0);
	ll ans=0;
	for(register ll i=0,pt=1,fac=1;i<=n;++i,pt=pt*2%mod,fac=fac*i%mod)
		ans=(ans+pt*fac%mod*f[i]%mod)%mod;
	printf("%lld",ans);
	return 0;
}