Stirling Number of the Second Kind
(math/stirling_second.hpp)

View this file on GitHub
Last update: 2024-01-08 00:27:17+09:00
Include: #include "math/stirling_second.hpp"

Description

第2種 Stirling 数 ${n \brace k}$ は，以下の恒等式で定義される数である．

\[{n \brace k} = \frac{1}{k!} \sum_{i=0}^n (-1)^{k-i} \binom{k}{i} i^n\]

${n \brace k}$ は，$n$ 個の区別できるボールを， $k$ 個の区別できない箱に，すべての箱に1つ以上のボールが入るように分配する方法の数である．

Operations

vector<T> stirling_second_table(int n)
- ${n \brace k}$ を各 $k=0,1,\dots,n$ について計算する
- 時間計算量: $O(n\log n)$
T stirling_second(int n, int k)
- ${n \brace k}$ を計算する
- 時間計算量: $O(n+k\log n)$

Notes

第2種 Stirling 数について以下の式が成り立つ．

\[x^n = \sum_{k=0}^n {n\brace k} x(x-1)\cdots(x-(k-1))\] \[{n\brace k} = {n-1\brace k-1} + k{n-1 \brace k}\]

Reference

「写像12相」を総整理！〜数え上げ問題の学びの宝庫〜 - Qiita

Depends on

Verified with

test/yosupo/stirling_number_of_the_second_kind.test.cpp

Code

#pragma once
#include <vector>
#include "../convolution/ntt.hpp"
#include "combination.cpp"

template <typename T>
std::vector<T> stirling_second_table(int n) {
    T f = 1;
    for (int i = 1; i <= n; ++i) f *= i;
    f = T(1) / f;
    std::vector<T> a(n + 1), b(n + 1);
    for (int i = n; i >= 0; --i) {
        a[i] = f * (i % 2 ? -1 : 1);
        b[i] = f * T(i).pow(n);
        f *= i;
    }
    auto c = convolution(a, b);
    return std::vector(c.begin(), c.begin() + n + 1);
}

template <typename T>
T stirling_second(int n, int k) {
    Combination<T> comb(n);
    T res = 0;
    for (int i = 0; i <= k; ++i) {
        T tmp = comb.comb(k, i) * T(i).pow(n);
        if ((k - i) & 1) res -= tmp;
        else res += tmp;
    }
    res /= comb.fact(k);
    return res;
}

#line 2 "math/stirling_second.hpp"
#include <vector>
#line 2 "convolution/ntt.hpp"
#include <bit>
#line 4 "convolution/ntt.hpp"

constexpr int get_primitive_root(int mod) {
    if (mod == 167772161) return 3;
    if (mod == 469762049) return 3;
    if (mod == 754974721) return 11;
    if (mod == 998244353) return 3;
    if (mod == 1224736769) return 3;
}

template <typename mint>
void ntt(std::vector<mint>& a) {
    constexpr int mod = mint::mod();
    constexpr mint primitive_root = get_primitive_root(mod);

    const int n = a.size();
    for (int m = n; m > 1; m >>= 1) {
        mint omega = primitive_root.pow((mod - 1) / m);
        for (int s = 0; s < n / m; ++s) {
            mint w = 1;
            for (int i = 0; i < m / 2; ++i) {
                mint l = a[s * m + i];
                mint r = a[s * m + i + m / 2];
                a[s * m + i] = l + r;
                a[s * m + i + m / 2] = (l - r) * w;
                w *= omega;
            }
        }
    }
}

template <typename mint>
void intt(std::vector<mint>& a) {
    constexpr int mod = mint::mod();
    constexpr mint primitive_root = get_primitive_root(mod);

    const int n = a.size();
    for (int m = 2; m <= n; m <<= 1) {
        mint omega = primitive_root.pow((mod - 1) / m).inv();
        for (int s = 0; s < n / m; ++s) {
            mint w = 1;
            for (int i = 0; i < m / 2; ++i) {
                mint l = a[s * m + i];
                mint r = a[s * m + i + m / 2] * w;
                a[s * m + i] = l + r;
                a[s * m + i + m / 2] = l - r;
                w *= omega;
            }
        }
    }
}

template <typename mint>
std::vector<mint> convolution(std::vector<mint> a, std::vector<mint> b) {
    const int size = a.size() + b.size() - 1;
    const int n = std::bit_ceil((unsigned int)size);
    a.resize(n);
    b.resize(n);
    ntt(a);
    ntt(b);
    for (int i = 0; i < n; ++i) a[i] *= b[i];
    intt(a);
    a.resize(size);
    mint n_inv = mint(n).inv();
    for (int i = 0; i < size; ++i) a[i] *= n_inv;
    return a;
}
#line 3 "math/combination.cpp"

template <typename mint>
class Combination {
   public:
    Combination() = default;
    Combination(int n) : fact_(n + 1), fact_inv_(n + 1) {
        fact_[0] = 1;
        for (int i = 1; i <= n; ++i) fact_[i] = fact_[i - 1] * i;
        fact_inv_[n]=fact_[n].inv();
        for (int i = n; i > 0; --i) fact_inv_[i - 1] = fact_inv_[i] * i;
    }

    mint perm(int n, int k) const {
        if (k < 0 || n < k) return 0;
        return fact_[n] * fact_inv_[n - k];
    }

    mint comb(int n, int k) const {
        if (k < 0 || n < k) return 0;
        return fact_[n] * fact_inv_[k] * fact_inv_[n - k];
    }

    mint fact(int n) const { return fact_[n]; }
    mint fact_inv(int n) const { return fact_inv_[n]; }

   private:
    std::vector<mint> fact_, fact_inv_;
};

template <typename T>
T comb(long long n, int k) {
    if (k < 0 || n < k) return 0;
    T num = 1, den = 1;
    for (int i = 1; i <= k; ++i) {
        num = num * (n - i + 1);
        den = den * i;
    }
    return num / den;
}
#line 5 "math/stirling_second.hpp"

template <typename T>
std::vector<T> stirling_second_table(int n) {
    T f = 1;
    for (int i = 1; i <= n; ++i) f *= i;
    f = T(1) / f;
    std::vector<T> a(n + 1), b(n + 1);
    for (int i = n; i >= 0; --i) {
        a[i] = f * (i % 2 ? -1 : 1);
        b[i] = f * T(i).pow(n);
        f *= i;
    }
    auto c = convolution(a, b);
    return std::vector(c.begin(), c.begin() + n + 1);
}

template <typename T>
T stirling_second(int n, int k) {
    Combination<T> comb(n);
    T res = 0;
    for (int i = 0; i <= k; ++i) {
        T tmp = comb.comb(k, i) * T(i).pow(n);
        if ((k - i) & 1) res -= tmp;
        else res += tmp;
    }
    res /= comb.fact(k);
    return res;
}

Stirling Number of the Second Kind (math/stirling_second.hpp)