Stirling Number1
(library/sequence/stirling_number1.hpp)

View this file on GitHub
Last update: 2026-06-28 22:26:07+09:00
Include: #include "library/sequence/stirling_number1.hpp"

Stirling Number1

stirling_number1_reversed

シグネチャ

template <typename mint>
std::vector<mint> stirling_number1_reversed(int n) // (1)
template <typename mint>
std::vector<mint> stirling_number1_reversed(const long long n, const int k) // (2)

概要

(符号なし) 第一種スターリング数 $\mathrm{S1}(n,\cdot)$ を逆順に並べた列 $A=(\mathrm{S1}(n,n),\mathrm{S1}(n,n-1),\ldots)$ を計算します．つまり，$A_i$ は集合 $\{0,\ldots,n-1\}$ から $i$ 個の要素を選んで積を取ったものの総和となります．形式的には，以下が成り立ちます．
\[A_i=\sum_{\overset{\scriptstyle S\subset\{0,\ldots,n-1\},}{|S|=i}}\prod_{j\in S} j\]
組合せ的には，$\mathrm{S1}(n,i)$ はラベル付けされた $n$ 個の玉を $i$ 個の円環 (巡回列) に分割する方法の個数と一致します．
テンプレート引数
- mint: modint 型を想定
返り値
1. $\{A_i\} _ {i=0} ^ n=\{\mathrm{S1}(n,n-i)\} _ {i=0} ^ {n}$
2. $\{A_i\} _ {i=0} ^ k=\{\mathrm{S1}(n,n-i)\} _ {i=0} ^ {k}$
  
  Note. $k>n$ の場合は $A_{n+1}=\cdots=A_{k}=0$ が保証される．
制約
1. - $0\leq n\leq 10 ^ 6$
2. - $0\leq n\leq 10 ^ {18}$
  - $0\leq k\leq 5000$
  - $k\lt \mathrm{mod}$
時間計算量
1. $O(n\log n)$
2. $O(k ^ 2\log n)$
参考
- スターリング数 - Wikipedia
- 第一種スターリング数の末尾項を計算する

stirling_number1

シグネチャ

template <typename mint>
std::vector<mint> stirling_number1(int n)

概要

stirling number1 reversed (1) の列を逆順にしたもの，つまり (符号なし) 第一種スターリング数の列 $\{\mathrm{S1}(n,i)\} _ {i=0} ^ n$ を計算します．

Depends on

Required by

Stirling Number of the First Kind (Small Prime Mod) (library/sequence/stirling_number1_small_prime_mod.hpp)

Verified with

Code

#ifndef SUISEN_STIRLING_NUMBER_1
#define SUISEN_STIRLING_NUMBER_1

#include <algorithm>
#include "library/math/inv_mods.hpp"
#include "library/math/factorial.hpp"

namespace suisen {
    /**
     * return:
     *   vector<mint> v s.t. v[i] = S1[n,n-i] for i=0,...,k (unsigned)
     * constraints:
     *   0 <= n <= 10^6
     */
    template <typename FPSType>
    std::vector<typename FPSType::value_type> stirling_number1_reversed(int n) {
        using mint = typename FPSType::value_type;
        factorial<mint> fac(n);
        int l = 0;
        while ((n >> l) != 0) ++l;
        FPSType a{ 1 };
        int m = 0;
        while (l-- > 0) {
            FPSType f(m + 1), g(m + 1);
            mint powm = 1;
            for (int i = 0; i <= m; ++i, powm *= m) {
                f[i] = powm * fac.fac_inv(i);
                g[i] = a[i] * fac.fac(m - i);
            }
            f *= g, f.cut(m + 1);
            for (int i = 0; i <= m; ++i) f[i] *= fac.fac_inv(m - i);
            a *= f, m *= 2, a.cut(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;
    }
    template <typename FPSType>
    std::vector<typename FPSType::value_type> stirling_number1(int n) {
        std::vector<typename FPSType::value_type> a(stirling_number1_reversed<FPSType>(n));
        std::reverse(a.begin(), a.end());
        return a;
    }
    /**
     * return:
     *   vector<mint> 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 mint>
    std::vector<mint> stirling_number1_reversed(const long long n, const int k) {
        inv_mods<mint> invs(k + 1);
        std::vector<mint> a(k + 1, 0);
        a[0] = 1;
        int l = 0;
        while (n >> l) ++l;
        mint m = 0;
        while (l-- > 0) {
            std::vector<mint> b(k + 1, 0);
            for (int j = 0; j <= k; ++j) {
                mint tmp = 1;
                for (int i = j; i <= k; ++i) {
                    b[i] += a[j] * tmp;
                    tmp *= (m - i) * invs[i - j + 1] * m;
                }
            }
            for (int i = k + 1; i-- > 0;) {
                mint sum = 0;
                for (int j = 0; j <= i; ++j) sum += a[j] * b[i - j];
                a[i] = sum;
            }
            m *= 2;
            if ((n >> l) & 1) {
                for (int i = k; i > 0; --i) a[i] += m * a[i - 1];
                ++m;
            }
        }
        return a;
    }
    template <typename mint>
    std::vector<std::vector<mint>> stirling_number1_table(int n) {
        std::vector dp(n + 1, std::vector<mint>{});
        for (int i = 0; i <= n; ++i) {
            dp[i].resize(i + 1);
            dp[i][0] = 0, dp[i][i] = 1;
            for (int j = 1; j < i; ++j) dp[i][j] = dp[i - 1][j - 1] + (i - 1) * dp[i - 1][j];
        }
        return dp;
    }
} // namespace suisen

#endif // SUISEN_STIRLING_NUMBER_1

#line 1 "library/sequence/stirling_number1.hpp"



#include <algorithm>
#line 1 "library/math/inv_mods.hpp"



#include <vector>

namespace suisen {
    template <typename mint>
    class inv_mods {
    public:
        inv_mods() = default;
        inv_mods(int n) { ensure(n); }
        const mint& operator[](int i) const {
            ensure(i);
            return invs[i];
        }
        static void ensure(int n) {
            int sz = invs.size();
            if (sz < 2) invs = { 0, 1 }, sz = 2;
            if (sz < n + 1) {
                invs.resize(n + 1);
                for (int i = sz; i <= n; ++i) invs[i] = mint(mod - mod / i) * invs[mod % i];
            }
        }
    private:
        static std::vector<mint> invs;
        static constexpr int mod = mint::mod();
    };
    template <typename mint>
    std::vector<mint> inv_mods<mint>::invs{};

    template <typename mint>
    std::vector<mint> get_invs(const std::vector<mint>& vs) {
        const int n = vs.size();

        mint p = 1;
        for (auto& e : vs) {
            p *= e;
            assert(e != 0);
        }
        mint ip = p.inv();

        std::vector<mint> rp(n + 1);
        rp[n] = 1;
        for (int i = n - 1; i >= 0; --i) {
            rp[i] = rp[i + 1] * vs[i];
        }
        std::vector<mint> res(n);
        for (int i = 0; i < n; ++i) {
            res[i] = ip * rp[i + 1];
            ip *= vs[i];
        }
        return res;
    }
}


#line 1 "library/math/factorial.hpp"



#include <cassert>
#line 6 "library/math/factorial.hpp"

namespace suisen {
    // 引数として与える値に対して、法が十分大きいことを仮定する
    template <typename T, typename U = T>
    struct factorial {
        factorial() = default;
        factorial(int n) { ensure(n); }

        static void ensure(const int n) {
            int sz = _fac.size();
            if (n + 1 <= sz) return;
            int new_size = std::max(n + 1, sz * 2);
            _fac.resize(new_size), _fac_inv.resize(new_size);
            for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
            _fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
            for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
        }

        T fac(const int i) {
            ensure(i);
            return _fac[i];
        }
        T operator()(int i) {
            return fac(i);
        }
        U fac_inv(const int i) {
            ensure(i);
            return _fac_inv[i];
        }
        // i の逆数
        // i = 0 の場合は assert 違反となる
        U inv(const int i) {
            assert(i > 0);
            ensure(i);
            return _fac_inv[i] * _fac[i - 1];
        }
        U binom(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
        }
        // binom(n, r) の逆数
        // binom(n, r) = 0 の場合は assert 違反となる
        U binom_inv(const int n, const int r) {
            assert(r >= 0 and n >= r);
            ensure(n);
            return _fac_inv[n] * _fac[r] * _fac[n - r];
        }
        // n 種類から重複を許して r 個選ぶ場合の数
        // x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数でもある
        // multichoose(n, r) = binom(n + r - 1, r)
        U multichoose(const int n, const int r) {
            if (n < 0 or r < 0) return 0;
            return r > 0 ? binom(n + r - 1, r) : U(1);
        }
        // n 種類から重複を許して r 個選ぶ場合の数 multichoose(n, r) の逆数
        // x_1+x_2+...+x_n=r（x_i は非負整数）となる x の個数の逆数でもある
        // multichoose(n, r) = binom(n + r - 1, r)
        // multichoose(n, r) = 0 の場合は assert 違反となる
        U multichoose_inv(const int n, const int r) {
            assert(n >= 0 and r >= 0);
            return r > 0 ? binom_inv(n + r - 1, r) : U(1);
        }
        template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>
        U polynom(const int n, const Ds& ...ds) {
            if (n < 0) return 0;
            ensure(n);
            int sumd = 0;
            U res = _fac[n];
            for (int d : { ds... }) {
                if (d < 0 or d > n) return 0;
                sumd += d;
                res *= _fac_inv[d];
            }
            if (sumd > n) return 0;
            res *= _fac_inv[n - sumd];
            return res;
        }
        U perm(const int n, const int r) {
            if (n < 0 or r < 0 or n < r) return 0;
            ensure(n);
            return _fac[n] * _fac_inv[n - r];
        }
        // perm(n, r) の逆数
        // perm(n, r) = 0 の場合は assert 違反となる
        U perm_inv(const int n, const int r) {
            assert(r >= 0 and n >= r);
            ensure(n);
            return _fac_inv[n] * _fac[n - r];
        }
    private:
        static std::vector<T> _fac;
        static std::vector<U> _fac_inv;
    };
    template <typename T, typename U>
    std::vector<T> factorial<T, U>::_fac{ 1 };
    template <typename T, typename U>
    std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen


#line 7 "library/sequence/stirling_number1.hpp"

namespace suisen {
    /**
     * return:
     *   vector<mint> v s.t. v[i] = S1[n,n-i] for i=0,...,k (unsigned)
     * constraints:
     *   0 <= n <= 10^6
     */
    template <typename FPSType>
    std::vector<typename FPSType::value_type> stirling_number1_reversed(int n) {
        using mint = typename FPSType::value_type;
        factorial<mint> fac(n);
        int l = 0;
        while ((n >> l) != 0) ++l;
        FPSType a{ 1 };
        int m = 0;
        while (l-- > 0) {
            FPSType f(m + 1), g(m + 1);
            mint powm = 1;
            for (int i = 0; i <= m; ++i, powm *= m) {
                f[i] = powm * fac.fac_inv(i);
                g[i] = a[i] * fac.fac(m - i);
            }
            f *= g, f.cut(m + 1);
            for (int i = 0; i <= m; ++i) f[i] *= fac.fac_inv(m - i);
            a *= f, m *= 2, a.cut(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;
    }
    template <typename FPSType>
    std::vector<typename FPSType::value_type> stirling_number1(int n) {
        std::vector<typename FPSType::value_type> a(stirling_number1_reversed<FPSType>(n));
        std::reverse(a.begin(), a.end());
        return a;
    }
    /**
     * return:
     *   vector<mint> 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 mint>
    std::vector<mint> stirling_number1_reversed(const long long n, const int k) {
        inv_mods<mint> invs(k + 1);
        std::vector<mint> a(k + 1, 0);
        a[0] = 1;
        int l = 0;
        while (n >> l) ++l;
        mint m = 0;
        while (l-- > 0) {
            std::vector<mint> b(k + 1, 0);
            for (int j = 0; j <= k; ++j) {
                mint tmp = 1;
                for (int i = j; i <= k; ++i) {
                    b[i] += a[j] * tmp;
                    tmp *= (m - i) * invs[i - j + 1] * m;
                }
            }
            for (int i = k + 1; i-- > 0;) {
                mint sum = 0;
                for (int j = 0; j <= i; ++j) sum += a[j] * b[i - j];
                a[i] = sum;
            }
            m *= 2;
            if ((n >> l) & 1) {
                for (int i = k; i > 0; --i) a[i] += m * a[i - 1];
                ++m;
            }
        }
        return a;
    }
    template <typename mint>
    std::vector<std::vector<mint>> stirling_number1_table(int n) {
        std::vector dp(n + 1, std::vector<mint>{});
        for (int i = 0; i <= n; ++i) {
            dp[i].resize(i + 1);
            dp[i][0] = 0, dp[i][i] = 1;
            for (int j = 1; j < i; ++j) dp[i][j] = dp[i - 1][j - 1] + (i - 1) * dp[i - 1][j];
        }
        return dp;
    }
} // namespace suisen

Stirling Number1 (library/sequence/stirling_number1.hpp)

Stirling Number1

stirling_number1_reversed

stirling_number1

Depends on

Required by

Verified with

Code

Stirling Number1
(library/sequence/stirling_number1.hpp)