P5408 第一类斯特林数·行 --- 题解

Fast_Gai_Tle · 2021-10-28 13:46:28 · 题解

第一类斯特林数的性质

定义 x^{\underline{n}} 为 x 的下降幂，即 \prod_{i=0}^{n-1}(x-i)，x^{\overline{n}} 为 x 的上升幂，即 \prod_{i=0}^{n-1}(x+i)。

x^{\underline{n}}=\sum_{i=1}^{n}(-1)^{n-i}\begin{bmatrix}n\\i\end{bmatrix}x^{i}

证明如下：当 n=1 时，x^{\underline{n}}=(-1)^{0}x=x，显然成立。

若 n=k 时成立，则当 n=k+1 时：
x^{\underline{k+1}}=(x-k)x^{\underline{k}} =(x-k)\sum_{i=1}^{k}(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}x^{i} =\sum_{i=1}^{k}(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}x^{i+1}-k\sum_{i=1}^{k}(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}x^{i} =\sum_{i=0}^{k}(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}x^{i+1}-k\sum_{i=1}^{k}(-1)^{k-i}\begin{bmatrix}k\\i\end{bmatrix}x^{i} =\sum_{i=1}^{k+1}(-1)^{k-i+1}\begin{bmatrix}k\\i-1\end{bmatrix}x^{i}+k\sum_{i=1}^{k}(-1)^{k-i+1}\begin{bmatrix}k\\i\end{bmatrix}x^{i} =x^{k+1}+\sum_{i=1}^{k}(-1)^{k-i+1}\begin{bmatrix}k\\i-1\end{bmatrix}x^{i}+k\sum_{i=1}^{k}(-1)^{k-i+1}\begin{bmatrix}k\\i\end{bmatrix}x^{i} =x^{k+1}+\sum_{i=1}^{k}(-1)^{k-i+1}\begin{bmatrix}k+1\\i\end{bmatrix}x^{i} =\sum_{i=1}^{k+1}(-1)^{k+1-i}\begin{bmatrix}k+1\\i\end{bmatrix}x^{i}

求第一类斯特林数

相信 O(n^2) 的递推大家都知道。

下面介绍一种 O(nlogn) 求第一类斯特林数某一行的方法：

不难发现 x^{\overline{n}}，即为 \begin{bmatrix}n\\i\end{bmatrix} 的生成函数，又因为有 x^{\overline{2n}}=x^{\overline{n}}(x+n)^{\overline{n}}，考虑用分治的方式来求 x^{\overline{n}}。

令 F_n(x)=x^{\overline{n}}=\sum_{i=0}^{n}\begin{bmatrix}n\\i\end{bmatrix}x^{i}，有：

F_n(x+n)=\sum_{i=0}^{n}\begin{bmatrix}n\\i\end{bmatrix}(x+n)^{i}

=\sum_{i=0}^{n}\begin{bmatrix}n\\i\end{bmatrix}\sum_{j=0}^{i}\begin{pmatrix}i\\j\end{pmatrix}x^{j}n^{i-j}

=\sum_{j=0}^{n}\frac{x^{j}}{j!}\sum_{i=j}^{n}\begin{bmatrix}n\\i\end{bmatrix}i!\frac{n^{i-j}}{(i-j)!}

~~这个柿子是不是看起来很好求？~~

令 a_i=\begin{bmatrix}n\\i\end{bmatrix}i!，b_i=\dfrac{n^{i}}{i!}，则原式可化为：

=\sum_{j=0}^{n}\frac{x^{j}}{j!}\sum_{i=j}^{n}a_ib_{i-j}

=\sum_{j=0}^{n}\frac{x^{j}}{j!}\sum_{i-k=j}a_ib_k

可以发现这个柿子后半部分很像一个卷积的形式。只可惜不是 i+k=j，这时候只需把 b 数组翻转过来求一次FFT，再把所求答案往左移 n 位即可算出 F_n(x+n)。至此 F_n(x+n) 就可以用 F_n(x) 来计算，而 F_{2n}(x) 就是 F_n(x) 与 F_n(x+n) 的卷积，分治即可。如果感觉有点糊，可以看看代码：

#include <cstdio>
#include <algorithm>
#define re register
using namespace std;
typedef long long LL;
const LL mod = 167772161LL;
const int N = 4e5 + 5;
LL g, gi, a[N], b[N], f[N], rev[N], fac[N], inv[N];
int n, l, li;
LL ksm(LL a, LL n) {
    LL res = 1LL, b = a;
    for(; n; n >>= 1LL, b = b * b % mod)
        if(n & 1LL) res = res * b % mod;
    return res;
}
void NTT(LL *A, bool op, int lim) {
    for(re int i = 0; i < lim; ++i)
        if(i < rev[i]) swap(A[i], A[rev[i]]);
    for(re int m = 1; m < lim; m <<= 1) {
        LL wn = ksm(op ? g : gi, (mod - 1LL) / (m << 1LL));
        for(re int L = m << 1, i = 0; i < lim; i += L) {
            LL w = 1LL;
            for(re int j = 0; j < m; ++j, w = w * wn % mod) {
                LL x = A[i + j], y = w * A[i + m + j] % mod;
                A[i + m + j] = (x - y + mod) % mod;
                A[i + j] = (x + y) % mod;
            } 
        }
    }
}
void solve(int x) {
    if(x == 1) { f[1] = 1; return ; }
    solve(x >> 1); int n = x >> 1; LL ba = 1LL;

    // 求出a数组和b数组
    for(re int i = 0; i <= n; ++i) a[i] = f[i] * fac[i] % mod;
    for(re int i = 0; i <= n; ++i, ba = ba * n % mod) b[i] = ba * inv[i] % mod;

    // 把b数组翻转，求a，b的卷积
    for(re int i = 0; i <= n / 2; ++i) swap(b[i], b[n - i]);
    li = 1; l = 0; while(li <= 2 * n) li <<= 1, ++l;
    for(re int i = 0; i < li; ++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
    for(re int i = n + 1; i < li; ++i) a[i] = b[i] = 0;
    NTT(a, 1, li); NTT(b, 1, li);
    for(re int i = 0; i < li; ++i) a[i] = a[i] * b[i];
    NTT(a, 0, li); LL invli = ksm(li, mod - 2LL);

    // 把求出的Fn(x+n)向左移n位
    for(re int i = 0; i <= n; ++i) a[i] = a[i + n] * invli % mod;

    // 求Fn(x)和Fn(x+n)的卷积
    for(re int i = n + 1; i <= 2 * n; ++i) a[i] = 0;
    for(re int i = 0; i <= n; ++i) a[i] = a[i] * inv[i] % mod;
    NTT(f, 1, li); NTT(a, 1, li);
    for(re int i = 0; i < li; ++i) f[i] = f[i] * a[i] % mod;
    NTT(f, 0, li);
    for(re int i = 0; i <= 2 * n; ++i) f[i] = f[i] * invli % mod;

    // 如果n不能被二整除，则需再乘一个(x+n-1)，也就是把多项式左移一位，再加上每一位乘(n-1)
    // 注意上述的x与递归的x不同，这里指多项式的自变量
    if(x & 1) {
        for(re int i = 2 * n + 1; i >= 1; --i) f[i] = (f[i - 1] + (x - 1) * f[i] % mod) % mod;
        f[0] = (x - 1) * f[0] % mod;
    }
}
int main() {
    scanf("%d", &n); g = 3LL; gi = ksm(g, mod - 2LL);
    fac[0] = fac[1] = inv[0] = inv[1] = 1LL;
    for(re int i = 2; i <= n; ++i) {
        fac[i] = fac[i - 1] * i % mod;
        inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    }
    for(re int i = 2; i <= n; ++i)
        inv[i] = inv[i - 1] * inv[i] % mod;
    solve(n);
    for(re int i = 0; i <= n; ++i)
        printf("%lld ", f[i]);
    puts("");
    return 0;
}

向学习更多斯特林数的芝士可以看这个