第一種スターリング数の末尾 k 項を計算する

こどふぉブログ https://codeforces.com/blog/entry/47764 の最後で紹介されているテクニックです．第一種スターリング数 $\begin{bmatrix}N\\ N-i\end{bmatrix}\bmod P$ を全ての $i=0,\ldots,K$ に対して $O(K ^ 2\log N +\log P)$ で計算できます．

解説

以下， $a _ {N, i}:=\begin{bmatrix}N\\ N-i\end{bmatrix}$ とおきます．

$\displaystyle \mathcal{S}_i(l,r)=\{S\subset\{l,\ldots,r-1\}\mid |S|=i\}$

とすると，任意の非負整数 $n$ に対して

$\displaystyle \begin{align} a_{2n,i} &=\sum_{S\in\mathcal{S}_i(0,2n)}\prod_{v\in S}v\\ &=\sum_{j=0}^i\left(\sum_{S\in\mathcal{S}_{i-j}(0,n)}\prod_{v\in S}v\right)\left(\sum_{S\in\mathcal{S}_j(n,2n)}\prod_{v\in S}v\right)\\ &=\sum_{j=0}^ia_{n,i-j}\left(\sum_{S\in\mathcal{S}_j(0,n)}\prod_{v\in S}(n+v)\right)\\ &=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}\left(\sum_{S\in\mathcal{S}_k(0,n)}\prod_{v\in S}v\right)\\ &=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}a_{n,k}\\ a_{2n+1,i}&=\begin{cases} 1&(i=0)\\ a_{2n,i}+2n\cdot a_{2n,i-1}&(i\gt 0)\ \end{cases} \end{align}$

が成り立つのでダブリング的に $a _ n$ を求めていくことが出来ます．

このままでは $a _ n$ から $a _ {2n}$ を求めるのに $\Theta(K ^ 3)$ 掛かってしまいそうですが， $\displaystyle b _ j := \sum _ {k=0} ^ j\begin{pmatrix}n-k \\ j - k \end{pmatrix} n ^ {j - k} a _ {n, k}$ は $i$ に依存しないことに注目します．すると，これは $O(K ^ 2)$ で前計算できるので，全体としても $O(K ^ 2)$ で $a _ {2n}$ を求めることができます．

ここで， $n$ が大きい場合は階乗の前計算が出来ないため二項係数の計算が少し難しくなる点に注意です． $b _ j$ に集める形 ( $j\rightarrow k$ 順のループ) で書こうとすると $n-k$ の値が動いて厄介なので， $b _ j$ に配る形 ( $k\rightarrow j$ 順のループ) で書くと楽だと思います．

以上を実装したのが次になります． $i=1,\ldots,K$ の $\mathrm{mod}\ P$ における乗法逆元の計算をサボっているので計算量は $\Theta(K ^ 2\log N + K\log P)$ となっていますが，実用上は誤差レベルだと思います．

実装例 (C++)

折り畳み

注意: $K\lt P$ を仮定しています．これが満たされない場合，逆元が存在しないことがあるので二項係数の計算がうまくいきません．

/**
 * return:
 *   vector<modint> v s.t. v[i] = S1[n,n-i] for i=0,...,k, where S1 is the stirling number of the first kind (unsigned).
 * constraints:
 * - 0 <= n <= 10^18
 * - 0 <= k <= 5000
 * - k < mod
 */
template <typename modint>
std::vector<modint> last_terms_of_stirling_numbers_of_the_first_kind(unsigned long long n, size_t k) {
    std::vector<modint> invs(k + 1);
    for (size_t i = 1; i <= k; ++i) invs[i] = modint(1) / i;
    std::vector<modint> a(k + 1, 0);
    a[0] = 1;
    size_t bit_length = 0;
    while (n >> bit_length != 0) ++bit_length;
    unsigned long long m = 0;
    for (size_t bit = bit_length; bit --> 0;) {
        std::vector<modint> b(k + 1, 0);
        for (size_t j = 0; j <= k; ++j) {
            modint tmp = 1;
            for (size_t i = j; i <= k; ++i) {
                b[i] += a[j] * tmp;
                tmp *= (m - i) * invs[i - j + 1] * m;
            }
        }
        for (size_t i = k + 1; i --> 0;) {
            modint sum = 0;
            for (size_t j = 0; j <= i; ++j) sum += a[j] * b[i - j];
            a[i] = sum;
        }
        m <<= 1;
        if ((n >> bit) & 1) {
            for (size_t i = k; i > 0; --i) a[i] += m * a[i - 1];
            ++m;
        }
    }
    return a;
}

verify: https://codeforces.com/contest/1516/submission/114464936

以上が元の記事に書かれていたことを私なりに再構築して説明した部分になります．

おまけ

紹介したテクニックを用いれば， $\begin{bmatrix}N\\ N-i\end{bmatrix}\bmod P$ を全ての $i=0,\ldots,N$ に対して $O(N\log N)$ で計算することができます．

$a _ {2n, i}$ を求める式をさらに変形します．

$\displaystyle \begin{align} a_{2n,i} &=\sum_{j=0}^ia_{n,i-j}\sum_{k=0}^j\begin{pmatrix}n-k\\ j-k\end{pmatrix}n^{j-k}a_{n,k}\\ &=\sum_{j=0}^ia_{n,i-j}\cdot\dfrac{1}{(n-j)!}\sum_{k=0}^j\dfrac{n^{j-k}}{(j-k)!}\cdot a_{n,k}\cdot(n-k)!\\ &=\sum_{j=0}^ia_{n,i-j}\cdot\dfrac{[x^j]f_n g_n}{(n-j)!}\quad\left(\begin{aligned} f_n:=&\sum_{k=0}^\infty \dfrac{n^k}{k!}x^k\\ g_n:=&\sum_{k=0}^n a_{n,k}\cdot(n-k)!\cdot x^k \end{aligned}\right)\\ &=[x^i]a_nh_n\quad\left(h_n:=\sum_{j=0}^n\dfrac{[x^j]f_ng_n}{(n-j)!}x^j\right) \end{align}$

上式より $a _ {2n} = a _ n h _ n$ が成り立ち，2 回の畳み込みで $a _ {2n}$ が求まることが分かります． $f_n$ は $n$ 次までの項で打ち切ってよいことから， $a_n$ から $a_{2n}$ を求める計算は $O(n\log n)$ でできることが分かりました．

$i$ 回目の反復で列の長さは $O(2 ^ i)$ であり，反復の回数は $O(\log N)$ 回であることから，全体の計算量は $O(N\log N)$ です．

実装例 (C++)

折り畳み

注意: 一部独自のクラスを用いているため，これをそのまま貼っても動きません

template <typename mint>
std::vector<mint> stirling_number1_reversed(int n) {
    std::vector<mint> fac(n + 1), fac_inv(n + 1);
    fac[0] = 1;
    for (int i = 1; i <= n; ++i) fac[i] = fac[i - 1] * i;
    fac_inv[n] = fac[n].inv();
    for (int i = n; i > 0; --i) fac_inv[i - 1] = fac_inv[i] * i;
    int l = 0;
    while ((n >> l) != 0) ++l;
    // 形式的冪級数のクラス
    FPS<mint> a {1};
    int m = 0;
    while (l --> 0) {
        // f, g は m 次の多項式
        FPS<mint> f(m + 1), g(m + 1);
        mint powm = 1;
        for (int i = 0; i <= m; ++i, powm *= m) {
            f[i] = powm * fac_inv[i];
            g[i] = a[i] * fac[m - i];
        }
        // f := f * g とし，(定数倍改善のために) f の m+1 次以降の係数を捨てる
        f *= g, f.pre_inplace(m);
        for (int i = 0; i <= m; ++i) f[i] *= fac_inv[m - i];
        // a = a * f とすることで a_{2n} が求められる
        // 次数 m を 2 倍し，m+1 次以降の係数を捨てる (必須)
        a *= f, m *= 2, a.pre_inplace(m);
        // a_{m} -> a_{m+1} の更新
        if ((n >> l) & 1) {
            a.push_back(0);
            for (int i = m; i > 0; --i) a[i] += m * a[i - 1];
            ++m;
        }
    }
    return a;
}

おまけ2

元の問題も同様に高速化することが出来ます． $K+1$ 次以上の項を無視してよいので，次のように $a_{2n}$ を求めることが出来ます．

$\displaystyle \begin{align} a_{2n} &=a_nh_n\quad\left(\begin{aligned} f_n:=&\sum_{k=0}^{\min(n,K)} \dfrac{n^k}{k!}x^k\\ g_n:=&\sum_{k=0}^{\min(n,K)} a_{n,k}\cdot(n-k)!\cdot x^k\\ h_n:=&\sum_{j=0}^{\min(n,K)} \dfrac{[x^j]f_ng_n}{(n-j)!}x^j \end{aligned}\right) \end{align}$

ここで， $g _ n$ や $h _ n$ の係数を求める際に最大で $N/2$ 程度の階乗を計算する必要があります． $N$ が小さい場合は前計算が出来るので問題とはなりませんでしたが， $N$ が大きくなってくるとそうもいきません．

これは簡単には求まりませんが，1 つの階乗を $O(\sqrt{P}(\log P) ^ 2)$ あるいは $O(\sqrt{P}\log P)$ で求めるアルゴリズムが存在します ¹．従って， $0!,\ldots,K!$ を前計算しておき，反復の度に都度 $(n-\min(n,K))!$ を高速階乗計算により求めることで全体 $O((K\log K + \sqrt{P}\log P)\log N)$ で問題を解くことが出来ます．