This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub suisen-cp/cp-library-cpp
#include "library/sequence/stirling_number1.hpp"
シグネチャ
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$ 個の要素を選んで積を取ったものの総和となります.形式的には,以下が成り立ちます.
組合せ的には,$\mathrm{S1}(n,i)$ はラベル付けされた $n$ 個の玉を $i$ 個の円環 (巡回列) に分割する方法の個数と一致します.
テンプレート引数
mint
返り値
$\{A_i\} _ {i=0} ^ k=\{\mathrm{S1}(n,n-i)\} _ {i=0} ^ {k}$
Note. $k>n$ の場合は $A_{n+1}=\cdots=A_{k}=0$ が保証される.
制約
時間計算量
参考
template <typename mint> std::vector<mint> stirling_number1(int n)
stirling number1 reversed (1) の列を逆順にしたもの,つまり (符号なし) 第一種スターリング数の列 $\{\mathrm{S1}(n,i)\} _ {i=0} ^ n$ を計算します.
#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]; } 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]; } 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]; } 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