[BZOJ4555][Tjoi2016&Heoi2016]求和(FFT)

=== ===

这里放传送门

=== ===

题解

首先值得吐槽的就是这个题如果不知道第二类斯特林数的通项公式的话是很难做的。。想想万一考试的时候遇到这个题然后不知道通项公式两眼抓瞎不是很GG。。而且题目还偏偏给了你递推公式。。明显就是在说【有本事你自己推嘛】之类的嘛。。。= =

然后ATP看到这个题第一反应就是去查了通项公式。。

S(n,m)=1m!k=1m(1)kCkm(mk)n

没错就决定是它了。
题目要求的式子是
i=1nj=1iS(i,j)×2j×j!

第二重枚举的上限跟 i 有关让人看起来很不爽,我们可以把枚举上界强行扩大到n。因为斯特林数一个直观意义就是把n个不同的小球放到m个相同的盒子里并且不允许有空的方案数,当 n 小于m的时候当然就是0。所以多枚举的那些对最后的答案不会有影响。

然后把通项公式代进去:

i=1nj=1n1j!k=1j(1)kCjk(jk)i×2j×j!

显然 j! 1j! 可以约掉了,那我们就把它约掉。然后怎么办呢?似乎看起来有什么 (1)k (jk)i 让我们能够联想到卷积之类的东西,但那个组合数好像很碍事,那我们把它化成阶乘的形式。
然后式子变成了:
i=1nj=1nk=1j(1)kj!k!(jk)!×(jk)i×2j

那个 2j j! 好像跟 k 无关诶那我们把它扔到第三层枚举外面去吧,然后把和k有关的量和跟 jk 有关的量写到一起去,把式子变成:

i=1nj=1n2j×j!k=1j(1)kk!×(jk)i(jk)!

这是卷积吗?好像不大是。。最大的不和谐因素就是后面那个 (jk)i ,这个东西跟i有关,但我们又不能做n遍卷积。。。可是当我们固定一个 j 的时候,i=1..n的所有值都被算了一遍。那么我们可以把 i 那个东西拿到里面去让它先把所有的都加起来。然后式子就变成了:
j=1n2j×j!k=1j(1)kk!×ni=1(jk)i(jk)!

那么我们设 f(i)=(1)ii! g(i)=k=1niki! h(x)=f(x)×g(x) ,那么式子就变成了

j=1n2j×j!×h(j)

一个卷积就可以搞出来啦!

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const long long Mod=998244353;
const long long G=3;
int n,L,R[300010],N,M;
long long a[300010],b[300010],mul[300010],ans,inv;
long long powww(long long a,int t){
    long long ans=1;a%=Mod;
    while (t!=0){
        if (t&1) ans=(ans*a)%Mod;
        a=(a*a)%Mod;t>>=1;
    }
    return ans;
}
void NTT(long long *a,int N,int opt){
    long long wn,w,x,y;
    for (int i=0;i<N;i++)
      if (i<R[i]) swap(a[i],a[R[i]]);
    for (int k=1;k<N;k<<=1){
        wn=powww(G,(Mod-1)/(k<<1));
        for (int p=(k<<1),i=0;i<N;i+=p){
            w=1;
            for (int j=0;j<k;j++){
                x=a[i+j];y=a[i+j+k]*w%Mod;
                a[i+j]=(x+y)%Mod;
                a[i+j+k]=(x-y+Mod)%Mod;
                w=(w*wn)%Mod;
            }
        }
    }
    if (opt==-1) reverse(a+1,a+N);
}
int main()
{
    scanf("%d",&n);mul[0]=1;
    for (int i=1;i<=n;i++)
      mul[i]=mul[i-1]*i%Mod;
    M=2*n;
    for (N=1;N<=M;N<<=1) L++;
    for (int i=0;i<=N;i++)
      R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=0;i<=n;i++){
        a[i]=powww(mul[i],Mod-2);
        if  (i&1) a[i]=Mod-a[i];
        b[i]=(1-powww(i,n+1))%Mod;
        b[i]=b[i]*powww(1-i,Mod-2)%Mod;
        b[i]=b[i]*powww(mul[i],Mod-2)%Mod;
        b[i]=(b[i]+Mod)%Mod;
    }
    b[1]=n+1;
    NTT(a,N,1);NTT(b,N,1);
    for (int i=0;i<=N;i++) a[i]=(a[i]*b[i])%Mod;
    NTT(a,N,-1);
    inv=powww(N,Mod-2)%Mod;
    for (int i=0;i<=n;i++){
        long long tmp=powww(2,i);
        tmp=tmp*mul[i]%Mod;
        tmp=tmp*a[i]%Mod*inv%Mod;
        ans=(ans+tmp)%Mod;
    }
    ans=(ans+Mod)%Mod;
    printf("%I64d\n",ans);
    return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值