算法学习FFT系列(3):多项式求逆详解——NTT+分治&&bzoj4555: [Tjoi2016&Heoi2016]求和例题详解

bzoj4555: [Tjoi2016&Heoi2016]

Description

在2016年，佳媛姐姐刚刚学习了第二类斯特林数，非常开心。
现在他想计算这样一个函数的值:

S(i, j)表示第二类斯特林数，递推公式为:
S(i, j) = j ∗ S(i − 1, j) + S(i − 1, j − 1), 1 <= j <= i − 1。
边界条件为：S(i, i) = 1(0 <= i), S(i, 0) = 0(1 <= i)
你能帮帮他吗?

Input

输入只有一个正整数

Output

输出f(n)。由于结果会很大，输出f(n)对998244353(7 × 17 × 223 + 1)取模的结果即可。1 ≤ n ≤ 100000

Sample Input

3

Sample Output

87

知识点：多项式求逆

基本概念

对于一个多项式$A(x)$，定义$degA$为多项式的度。
对于多项式$A(x),B(x)$存在唯一$Q(x),R(x)$使得

$$A(x)=B(x)Q(x)+R(x)(degR<degB)$$
我们就称 $Q(x)$为$A(x)$除以$B(x)$的商，$R(x)$为其余数，记作
$A(x) \equiv R(x) (\mod B(x))$

多项式的逆元

对于一个多项式$A(x)$如果存在唯一的$B(x)$且$degB\le degA$，有

$$A(x)B(x) \equiv 1 (\mod x^n)$$
我们称 $B(x)$为 $A(x) \mod x^n$意义下的逆元，记作 $A^{-1}(x)$

求解过程

考虑如何求解$A^{-1}(x)$，如果$n=1$，$A^{-1}(x)\equiv c^{-1} (\mod x)$，其中c为常数。这样$A^{-1} = c^{-1}$
对于$n> 1$的情况，类比FFT的分治思想，我们先考虑已经做完$\mod x^{\frac{n}{2}}$意义下的情况拓展到$\mod x^n$意义下的答案。
由定义可知

$$(1) A(x)B(x) \equiv 1 (\mod x^n)$$
由当前条件有
$$(2) A(x)B'(x) \equiv 1 (\mod x^{\frac{n}{2}})$$
将(1)式放到 $\mod \frac{n}{2}$意义下可得
$$(3) A(x)B(x) \equiv 1 (\mod x^{\frac{n}{2}})$$
我们由(2)-(3)可知
$$B(x)-B'(x) \equiv 0 (\mod x^{frac{n}{2}})$$
两边平方可得
$$B^2(x)-2*B'(x)B(x)+B'^2(x) \equiv 0 (\mod x^n)$$
为什么模数会平方？因为左边多项式0~n-1次方每一项系数都为0，平方之后考虑0~2n-1项，由多项式乘法可得其必定为零。
同乘 $A(x)$可以得到
$B(x) \equiv 2*B'(x)-A(x)B'^2(x) (\mod \frac{n}{2})$
于是我们就可以用FFT或NTT加速优化这个过程就可以结束这个算法啦！

复杂度分析

$T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)$
我有一个绝妙的证明方法，可惜我太弱，不会写

分析

可惜这道题，不是模板题。
这道题其实有两种算法，第一种是暴力拆斯特林数，拆完搞搞式子可以直接卷积。这里不介绍。
第二种算法就是观察式子的意义。
$ans=\sum\limits_i^n \sum\limits_j^i S(i,j)*j!*2^j$
设$g(n)= \sum\limits_j^n S(n,j)*j!*2^j$
我们考虑这个式子的意义，就是从n个集合里面挑j个数出来，考虑顺序，每种集合有两种状态。
于是我们尝试写出这个式子的递推式。
$g(n)= \sum\limits_j^n 2*C_n^jg(n-j)$
也就是枚举第一个集合的数的个数j，先从挑j个数出来，乘上状态即可。
然后神奇的事情发生了。
$g(n)=\sum\limits_j^n2*\frac{n!}{j!(n-j)!}g(n-j)$
$\frac{g(n)}{n!}=\sum\limits_j^n \frac{2}{j!} \cdot \frac{g(n-j)}{(n-j)!}$
我们把C拆开，发现了指数型生成函数的形式。
于是，构造地
令$G(x)=\sum\limits_i^ {\infty}\frac{g(i)}{i!} \cdot x^i, H(x)=\sum\limits_i^{\infty}\frac{2}{i!}\cdot x^i$
于是有$G=G*H+1$
然后$G=\frac{1}{1-H}$
多项式求逆即可。

代码

/**************************************************************
    Problem: 4555
    User: 2014lvzelong
    Language: C++
    Result: Accepted
    Time:4372 ms
    Memory:15628 kb
****************************************************************/

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 524288, P = 998244353, K = 18;
int read() {
    char ch = getchar(); int x = 0;
    while(ch < '0' || ch > '9') ch = getchar();
    for(;ch >= '0' && ch <= '9'; ch = getchar()) x = (x << 1) + (x << 3) - '0' + ch;
    return x;
}
int F[N], R[N], inv[N + 1], g[K + 1], ng[K + 1], tp[N], b[N], a[N], f[N];
int pow(int a, int x, int P) {
    int b = 1;
    for(; x;x >>= 1,  a = 1LL * a * a % P) if(x & 1) b = 1LL * b * a % P;
    return b; 
}
int Inv(int x) {return pow(x, P - 2, P);}
void NTT(int *F, int n, int f) {
    for(int i = 0;i < n; ++i) if(i < R[i]) swap(F[i], F[R[i]]);
    for(int d = 0;(1 << d) < n; ++d) {
        int wn = ~f ? g[d] : ng[d], m = 1 << d, m2 = m << 1;
        for(int j = 0;j < n; j += m2)
            for(int w = 1, l = 0; l < m; ++l, w = 1LL * w * wn % P)  {
                int &x = F[j + l], &y = F[m + j + l], t = 1LL * w * y % P;
                y = x - t; if(y < 0) y += P;
                x = x + t; if(x >= P) x -= P;
            }
    }
    if(!(~f)) for(int i = 0;i < n; ++i) F[i] = 1LL * F[i] * inv[n] % P;
}

void Work(int n) {
    if(n == 1) return void(b[0] = Inv(a[0]));
    Work(n >> 1);
    for(int i = 0;i < n; ++i) tp[i] = a[i], tp[n + i] = 0;
    int L = 0; while(!(n >> L & 1)) ++L;
    for(int i = 1;i < (n << 1); ++i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << L);
    NTT(tp, n << 1, 1); NTT(b, n << 1, 1);
    for(int i = 0;i < (n << 1); ++i) tp[i] = 1LL * b[i] * (2 + P - 1LL * tp[i] * b[i] % P) % P;
    NTT(tp, n << 1, -1);
    for(int i = 0;i < n; ++i) b[i] = tp[i], b[n + i] = 0;
}

int main() {
    int i, n = read(), m;
    for(g[K] = pow(3, (P - 1) / N, P), ng[K] = Inv(g[K]), i = K - 1; ~i; --i)
        g[i] = 1LL * g[i + 1] * g[i + 1] % P, ng[i] = 1LL * ng[i + 1] * ng[i + 1] % P;
    for(inv[1] = 1, i = 2;i <= N; ++i) inv[i] = 1LL * (P - inv[P % i]) * (P / i) % P;
    f[0] = 1; for(i = 1;i <= n; ++i) f[i] = 1LL * f[i - 1] * inv[i] % P;
    a[0] = 1; for(i = 1;i <= n; ++i) a[i] = ((P - f[i]) << 1) % P; 
    for(m = 1; m <= n; m <<= 1) ;
    Work(m);
    int ans = b[n];
    for(int i = n; i; --i)   ans = (1LL * ans * i + b[i - 1]) % P;
    printf("%d\n", ans);
    return 0;
}