BZOJ 4555: [Tjoi2016&Heoi2016]求和 [分治FFT 组合计数 | 多项式求逆]

4555: [Tjoi2016&Heoi2016]求和

题意：求\[ \sum_{i=0}^n \sum_{j=0}^i S(i,j)\cdot 2^j\cdot j! \\ S是第二类斯特林数 \]

---

首先你要把这个组合计数肝出来，于是我去翻了一波《组合数学》

用斯特林数容斥原理推导那个式子可以直接出卷积形式，见下一篇，本篇是分治fft做法

组合计数

斯特林数 \(S(n,i)\)表示将n个不同元素划分成i个相同集合非空的方案数
Bell数 \(B(n)=\sum\limits_{i=0}^n S(n,i)\)也就是分成任意个相同的集合
有一个递推式
\[ B(n) = \sum_{i=0}^{n-1} \binom{n-1}{n-1-i}B(i) = \sum_{i=0}^{n-1} \binom{n-1}{i}B(i) \]
\(推导\)：
考虑和第一个元素在同一个集合的元素有几个，剩下的是子问题
或者这么想，因为每个集合相同，我们规定第一个元素在一个集合后这个集合就不相同了，剩下的元素分入这个不相同的集合或者成为子问题

\(S'(n,i)=S(n,i)*i!\)就是集合不相同啦，
\[ B'(n) = \sum_{i=1}^n \binom{n}{i}B'(n-i) \]
\(推导\)：
因为每个集合都不相同，我们没有必要规定第一个元素在集合里了，直接枚举第一个集合元素个数就可以了

然后那个\(2^j\)可以认为每个集合有两种颜色，枚举当前集合顺便枚举颜色，乘上一个2就行了

令\(f(n)=2 \sum\limits_{i=1}^n \binom{n}{i} f(n-i)\)，答案就是\(\sum_{i=0}^n f(i)\)
快速求出f就可以做了，整理一下发现
\[ f(n) = 2 n! \sum_{i=1}^n \frac{1}{i!} \frac{f(n-i)}{(n-i)!} \]
是一个卷积的形式，但是两边都有f，所以要用CDQ分治套FFT

分治fft

和CDQ分治一样，每次用\([l,mid]\)更新\([mid+1,r]\)
可以发现，\(f(i)[mid+1,r]\)都是\(\frac{f(i)}{i!}[l,mid]\)卷上\(\frac{1}{i!}[1,r-l]\)的结果
我们把两个函数\([l,mid]\)和\([1,r-l]\)的部分拿出来做卷积，更新\([mid+1,r]\)
可以认为先把向量移动了l位

复杂度\(O(nlog^2n)\)，13560ms，貌似别人的分治fft更慢20000ms多

注意初始值\(f[0]=1\)因为\(S(0,0)=1\)，以及cdq(0,n)

多项式求逆

当然也可以，具体看51nod这篇讨论吧
把\(f_0 x^0\)单独拿出来的思想比较有意思,\(A(x)\)从0开始，\(B(x)\)从1开始（\(B_0 = 0\)）
\[A(x) = F_0 x^0 + A(x) B(x)\]
update 2017.5.3 : 发现tls貌似补了一句\(2x e^x\)，求这个东西来做也可以，但不如直接处理处B来方便

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=(1<<18)+5, INF=1e9;
const ll P=998244353, g=3;
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

ll Pow(ll a, ll b) {
    ll ans=1;
    for(; b; b>>=1, a=a*a%P)
        if(b&1) ans=ans*a%P;
    return ans;
}

int n, rev[N];
ll inv[N], fac[N], facInv[N];
ll f[N], a[N], b[N];
void dft(ll *a, int n, int flag) {
    for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
    for(int l=2; l<=n; l<<=1) {
        int m=l>>1;
        ll wn = Pow(g, flag==1 ? (P-1)/l : P-1-(P-1)/l);
        for(ll *p=a; p!=a+n; p+=l) {
            ll w=1;
            for(int k=0; k<m; k++) {
                ll t = w*p[k+m];
                p[k+m]=(p[k]-t+P)%P;
                p[k]=(p[k]+t)%P;
                w=w*wn%P;
            }
        }
    }
    if(flag==-1) {
        ll inv=Pow(n, P-2);
        for(int i=0; i<n; i++) a[i]=a[i]*inv%P;
    }
}
void cdq(int l, int r) {
    if(l==r) return;
    int mid=(l+r)>>1;
    cdq(l, mid);
    int lim=r-l+1, n=1, k=0;
    while(n<lim) n<<=1, k++;
    for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));

    for(int i=0; i<n; i++) a[i]=b[i]=0;
    for(int i=l; i<=mid; i++) a[i-l] = f[i];
    for(int i=0; i<=r-l; i++) b[i] = facInv[i];
    dft(a, n, 1); dft(b, n, 1);
    for(int i=0; i<n; i++) a[i]=a[i]*b[i]%P;
    dft(a, n, -1);

    for(int i=mid+1; i<=r; i++) f[i]=(f[i] + 2*a[i-l])%P;
    cdq(mid+1, r);
}
int main() {
    freopen("in","r",stdin);
    n=read();
    inv[1]=1; fac[0]=facInv[0]=1;
    for(int i=1; i<=n; i++) {
        if(i!=1) inv[i] = (P-P/i)*inv[P%i]%P;
        fac[i] = fac[i-1]*i%P;
        facInv[i] = facInv[i-1]*inv[i]%P;
    }
    f[0]=1;
    cdq(0, n);
    ll ans=0;
    for(int i=0; i<=n; i++) ans = (ans + f[i]*fac[i]%P)%P;
    if(ans<0) ans+=P;
    printf("%lld\n", ans);
}