bzoj4555: 求和sum 快速傅立叶变换

题目大意

给定\(S(n,m)\)表示第二类斯特林数,定义函数\(f(n)\)
\[f(n) = \sum_{i=0}^n\sum_{j=0}^iS(i,j)*2^j*(j!)\]
给定正整数\(n,(n\leq 10^5)\),求\(f(n)\)

题解

我们都知道第二类斯特林数的递推公式为
\[S(i,j) = S(i-1,j-1) + j*S(i-1,j),(1 \leq j \leq i-1)\]
且有边界\(S(i,i) = 1(0 \leq i),S(i,0) = 0(1 \leq i)\)
第二类斯特林数\(S(i,j)\)的含义是把\(i\)个元素划分成\(j\)个无序的集合的方案
假设允许空集合的存在的话,方案即为\(m^n\)
我们应用容斥原理,枚举至少有多少空集合空集合,那么有
\[S(n,m) = \frac{1}{m!}\sum_{k=1}^{m}C_m^k(m-k)^n(-1)^k\]
设\(g(n) = \sum_{i=0}^nS(n,i)2^i(i!)\)
那么我们将\(S(n,m)\)代入\(g(n)\)化简得
\[g(n) = \sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]
那么将\(g(n)\)带入答案表达式中,有
\[ans = \sum_{n=0}^x\sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]
这时我们发现每次最外层的\(n -> (n+1)\)时,都相当于在内部的\(\frac{(m-k)^n}{(m-k)!}\)一项上又加上了一个\(\frac{(m-k)^{n+1}}{(m-k)!}\)
所以我们把这一项做等比数列求和.
设\(g(x) = \frac{x^{n+1} - x}{(x-1)(x!)}\)
那么上式变成了
\[ans = \sum_{m=0}^n2^m(m!)\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k)\]
于是我们在\(\sum_{k=0}^m\frac{(-1)^k}{k!}g(m-k)\)进行FFT计算卷积
这样就只剩下了一个sigma式,for循环一边即可.
复杂度\(O(nlogn)\)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
template<typename T>inline void read(T &x){
    x=0;char ch;bool flag = false;
    while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
    while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
const int maxn = 600010;
const int mod = 998244353;
const int pri_rt = 3;
int w[maxn];
inline int qpow(int x,int p){
    int ret = 1;
    for(;p;p>>=1,x=1LL*x*x%mod) if(p&1) ret=1LL*ret*x % mod;
    return ret;
}
inline void FNT(int *x,int n,int p){
    for(int i=0,t=0;i<n;++i){
        if(i > t) swap(x[i],x[t]);
        for(int j=n>>1;(t^=j)<j;j>>=1);
    }
    for(int m=2;m<=n;m<<=1){
        int k = m>>1;
        int wn = qpow(pri_rt,p == 1 ? (mod-1)/m : (mod-1) - (mod-1)/m);
        for(int i=1;i<k;++i) w[i] = 1LL*w[i-1]*wn % mod;
        w[0] = 1;
        for(int i=0;i<n;i+=m){
            for(int j=0;j<k;++j){
                int u = 1LL*x[i+j+k]*w[j] % mod;
                x[i+j+k] = x[i+j] - u;
                if(x[i+j+k] < 0) x[i+j+k] += mod;
                x[i+j] += u;
                if(x[i+j] >= mod) x[i+j] -= mod;
            }
        }
    }
    if(p == -1){
        int inv = qpow(n,mod-2);
        for(int i=0;i<n;++i) x[i] = 1LL*x[i]*inv % mod;
    }
}
int fac[maxn],inv[maxn];
inline void init(int n){
    fac[0] = 1;
    for(int i=1;i<=n;++i) fac[i] = 1LL*fac[i-1]*i % mod;
    inv[n] = qpow(fac[n],mod-2);
    for(int i = n-1;i>=0;--i) inv[i] = 1LL*inv[i+1]*(i+1) % mod;
}
int A[maxn],B[maxn];
int main(){
    int n;read(n);
    int len;for(len=1;len <= (n+1);len<<=1);len<<=1;
    init(n);
    for(int i=0;i<=n;++i){
        if(i&1) A[i] = -inv[i] + mod;
        else A[i] = inv[i];
    }
    for(int i=2;i<=n;++i){
        B[i] = qpow(i,n+1) - i + mod;
        if(B[i] < 0) B[i] += mod;
        B[i] = (1LL*B[i]*qpow(i-1,mod-2)%mod*inv[i]) % mod;
    }B[1] = n;
    FNT(A,len,1);FNT(B,len,1);
    for(int i=0;i<len;++i) A[i] = 1LL*A[i]*B[i] % mod;
    FNT(A,len,-1);
    int ans = 1;
    for(int i=1,f2=2;i<=n;++i){
        ans = (ans + 1LL*A[i]*f2%mod*fac[i]) % mod;
        f2 = (f2<<1) % mod;
    }printf("%d\n",ans);
    getchar();getchar();
    return 0;
}

转载于:https://www.cnblogs.com/Skyminer/p/6402254.html