Sequences (4): Stirling numbers of the second kind

Stirling Number of the Second Kind

Stirling numbers of the second kind(number of Sterling subsets) $\begin{Bmatrix}n\\ k\end{Bmatrix}$ , can also be recorded as $S(n,k)$ , indicating that the $n$ Pairwise different elements are divided into $k$ The number of alternatives that are indistinguishable from each other and are non-empty subsets.

Recursive

\begin{Bmatrix}n\\ k\end{Bmatrix}=\begin{Bmatrix}n-1\\ k-1\end{Bmatrix}+k\begin{Bmatrix}n-1\\ k\end{Bmatrix}

The boundary is $\begin{Bmatrix}n\\ 0\end{Bmatrix}=[n=0]$ 。

Consider proof in combinatorial sense.

When we insert a new element, we have two options:

Put new elements into a subset separately, with $\begin{Bmatrix}n-1\\ k-1\end{Bmatrix}$ kind of plan;
Put the new element into an existing non-empty subset, with $k\begin{Bmatrix}n-1\\ k\end{Bmatrix}$ kind of plan.

According to the principle of addition, the recursive formula can be obtained by adding the two formulas.

general formula

\begin{Bmatrix}n\\m\end{Bmatrix}=\sum\limits_{i=0}^m\dfrac{(-1)^{m-i}i^n}{i!(m-i)!}

Prove this formula using the inclusion-exclusion principle. Set general $n$ two different elements, divided into $k$ The number of solutions for pairwise different sets (empty sets are allowed) is $G_i$ , will $n$ two different elements, divided into $k$ The number of solutions for pairwise different non-empty sets (empty sets are not allowed) is $F_i$ 。

obviously

G_i=k^n\\ G_i=\sum\limits_{j=0}^i\binom{i}{j}F_j

According to the binomial inversion

\begin{aligned} F_i&=\sum\limits_{j=0}^{i}(-1)^{i-j}\binom{i}{j}G_j\\ &=\sum\limits_{j=0}^{i}(-1)^{i-j}\binom{i}{j}j^n\\ &=\sum\limits_{j=0}^{i}\dfrac{i!(-1)^{i-j}j^n}{j!(i-j)!} \end{aligned}

consider $F_i$ with $\begin{Bmatrix}n\\i\end{Bmatrix}$ relationship. Stirling numbers of the second kind require sets to be indistinguishable from each other, so $F_i$ Exactly $\begin{Bmatrix}n\\i\end{Bmatrix}$ of $i!$ times. So

\begin{Bmatrix}n\\m\end{Bmatrix}=\dfrac{F_m}{m!}=\sum\limits_{i=0}^m\dfrac{(-1)^{m-i}i^n}{i!(m-i)!}

Calculation of Stirling numbers of the second kind in the same row

The second type of Stirling numbers in "the same row" refers to those with different $i$ , the same $n$ a series of $\begin{Bmatrix}n\\i\end{Bmatrix}$ . Find all Stirling numbers of the second kind in the same row, which is $i=0..n$ Find out the general $n$ different elements are divided into $i$ The number of options for a non-empty set.

Method 1. Directly use the general formula

According to the general formula given above, convolution calculation is enough. The time complexity of this approach is $O(n \log n)$ 。

The code below uses a file called poly The polynomial class is for reference only.

 #ifndef _FEISTDLIB_POLY_ #define _FEISTDLIB_POLY_
/* * This file is part of the fstdlib project. * Version: Build v0.0.2 * You can check for details at https://github.com/FNatsuka/fstdlib */ #include <algorithm>
#include <cmath>
#include <cstdio>
#include <vector> namespace fstdlib { typedef long long ll;
int mod = 998244353, grt = 3; class poly { private: std::vector<int> data; void out(void) { for (int i = 0; i < (int)data.size(); ++i) printf("%d ", data[i]); puts(""); } public: poly(std::size_t len = std::size_t(0)) { data = std::vector<int>(len); } poly(const std::vector<int> &b) { data = b; } poly(const poly &b) { data = b.data; } void resize(std::size_t len, int val = 0) { data.resize(len, val); } std::size_t size(void) const { return data.size(); } void clear(void) { data.clear(); }
#if __cplusplus >= 201103L void shrink_to_fit(void) { data.shrink_to_fit(); }
#endif int &operator[](std::size_t b) { return data[b]; } const int &operator[](std::size_t b) const { return data[b]; } poly operator*(const poly &h) const; poly operator*=(const poly &h); poly operator*(const int &h) const; poly operator*=(const int &h); poly operator+(const poly &h) const; poly operator+=(const poly &h); poly operator-(const poly &h) const; poly operator-=(const poly &h); poly operator<<(const std::size_t &b) const; poly operator<<=(const std::size_t &b); poly operator>>(const std::size_t &b) const; poly operator>>=(const std::size_t &b); poly operator/(const int &h) const; poly operator/=(const int &h); poly operator==(const poly &h) const; poly operator!=(const poly &h) const; poly operator+(const int &h) const; poly operator+=(const int &h); poly inv(void) const; poly inv(const int &h) const; friend poly sqrt(const poly &h); friend poly log(const poly &h); friend poly exp(const poly &h);
}; int qpow(int a, int b, int p = mod) { int res = 1; while (b) { if (b & 1) res = (ll)res * a % p; a = (ll)a * a % p, b >>= 1; } return res;
} std::vector<int> rev; void dft_for_module(std::vector<int> &f, int n, int b) { static std::vector<int> w; w.resize(n); for (int i = 0; i < n; ++i) if (i < rev[i]) std::swap(f[i], f[rev[i]]); for (int i = 2; i <= n; i <<= 1) { w[0] = 1, w[1] = qpow(grt, (mod - 1) / i); if (b == -1) w[1] = qpow(w[1], mod - 2); for (int j = 2; j < i / 2; ++j) w[j] = (ll)w[j - 1] * w[1] % mod; for (int j = 0; j < n; j += i) for (int k = 0; k < i / 2; ++k) { int p = f[j + k], q = (ll)f[j + k + i / 2] * w[k] % mod; f[j + k] = (p + q) % mod, f[j + k + i / 2] = (p - q + mod) % mod; } }
} poly poly::operator*(const poly &h) const { int N = 1; while (N < (int)(size() + h.size() - 1)) N <<= 1; std::vector<int> f(this->data), g(h.data); f.resize(N), g.resize(N); rev.resize(N); for (int i = 0; i < N; ++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? N >> 1 : 0); dft_for_module(f, N, 1), dft_for_module(g, N, 1); for (int i = 0; i < N; ++i) f[i] = (ll)f[i] * g[i] % mod; dft_for_module(f, N, -1), f.resize(size() + h.size() - 1); for (int i = 0, inv = qpow(N, mod - 2); i < (int)f.size(); ++i) f[i] = (ll)f[i] * inv % mod; return f;
} poly poly::operator*=(const poly &h) { return *this = *this * h; } poly poly::operator*(const int &h) const { std::vector<int> f(this->data); for (int i = 0; i < (int)f.size(); ++i) f[i] = (ll)f[i] * h % mod; return f;
} poly poly::operator*=(const int &h) { for (int i = 0; i < (int)size(); ++i) data[i] = (ll)data[i] * h % mod; return *this;
} poly poly::operator+(const poly &h) const { std::vector<int> f(this->data); if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] + h[i]) % mod; return f;
} poly poly::operator+=(const poly &h) { std::vector<int> &f = this->data; if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] + h[i]) % mod; return f;
} poly poly::operator-(const poly &h) const { std::vector<int> f(this->data); if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] - h[i] + mod) % mod; return f;
} poly poly::operator-=(const poly &h) { std::vector<int> &f = this->data; if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] - h[i] + mod) % mod; return f;
} poly poly::operator<<(const std::size_t &b) const { std::vector<int> f(size() + b); for (int i = 0; i < (int)size(); ++i) f[i + b] = data[i]; return f;
} poly poly::operator<<=(const std::size_t &b) { return *this = (*this) << b; } poly poly::operator>>(const std::size_t &b) const { std::vector<int> f(size() - b); for (int i = 0; i < (int)f.size(); ++i) f[i] = data[i + b]; return f;
} poly poly::operator>>=(const std::size_t &b) { return *this = (*this) >> b; } poly poly::operator/(const int &h) const { std::vector<int> f(this->data); int inv = qpow(h, mod - 2); for (int i = 0; i < (int)f.size(); ++i) f[i] = (ll)f[i] * inv % mod; return f;
} poly poly::operator/=(const int &h) { int inv = qpow(h, mod - 2); for (int i = 0; i < (int)data.size(); ++i) data[i] = (ll)data[i] * inv % mod; return *this;
} poly poly::inv(void) const { int N = 1; while (N < (int)(size() + size() - 1)) N <<= 1; std::vector<int> f(N), g(N), d(this->data); d.resize(N), f[0] = qpow(d[0], mod - 2); for (int w = 2; w < N; w <<= 1) { for (int i = 0; i < w; ++i) g[i] = d[i]; rev.resize(w << 1); for (int i = 0; i < w * 2; ++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? w : 0); dft_for_module(f, w << 1, 1), dft_for_module(g, w << 1, 1); for (int i = 0; i < w * 2; ++i) f[i] = (ll)f[i] * (2 + mod - (ll)f[i] * g[i] % mod) % mod; dft_for_module(f, w << 1, -1); for (int i = 0, inv = qpow(w << 1, mod - 2); i < w; ++i) f[i] = (ll)f[i] * inv % mod; for (int i = w; i < w * 2; ++i) f[i] = 0; } f.resize(size()); return f;
} poly poly::operator==(const poly &h) const { if (size() != h.size()) return 0; for (int i = 0; i < (int)size(); ++i) if (data[i] != h[i]) return 0; return 1;
} poly poly::operator!=(const poly &h) const { if (size() != h.size()) return 1; for (int i = 0; i < (int)size(); ++i) if (data[i] != h[i]) return 1; return 0;
} poly poly::operator+(const int &h) const { poly f(this->data); f[0] = (f[0] + h) % mod; return f;
} poly poly::operator+=(const int &h) { return *this = (*this) + h; } poly poly::inv(const int &h) const { poly f(*this); f.resize(h); return f.inv();
} int modsqrt(int h, int p = mod) { return 1; } poly sqrt(const poly &h) { int N = 1; while (N < (int)(h.size() + h.size() - 1)) N <<= 1; poly f(N), g(N), d(h); d.resize(N), f[0] = modsqrt(d[0]); for (int w = 2; w < N; w <<= 1) { g.resize(w); for (int i = 0; i < w; ++i) g[i] = d[i]; f = (f + f.inv(w) * g) / 2; f.resize(w); } f.resize(h.size()); return f;
} poly log(const poly &h) { poly f(h); for (int i = 1; i < (int)f.size(); ++i) f[i - 1] = (ll)f[i] * i % mod; f[f.size() - 1] = 0, f = f * h.inv(), f.resize(h.size()); for (int i = (int)f.size() - 1; i > 0; --i) f[i] = (ll)f[i - 1] * qpow(i, mod - 2) % mod; f[0] = 0; return f;
} poly exp(const poly &h) { int N = 1; while (N < (int)(h.size() + h.size() - 1)) N <<= 1; poly f(N), g(N), d(h); f[0] = 1, d.resize(N); for (int w = 2; w < N; w <<= 1) { f.resize(w), g.resize(w); for (int i = 0; i < w; ++i) g[i] = d[i]; f = f * (g + 1 - log(f)); f.resize(w); } f.resize(h.size()); return f;
} struct comp { long double x, y; comp(long double _x = 0, long double _y = 0) : x(_x), y(_y) {} comp operator*(const comp &b) const { return comp(x * b.x - y * b.y, x * b.y + y * b.x); } comp operator+(const comp &b) const { return comp(x + b.x, y + b.y); } comp operator-(const comp &b) const { return comp(x - b.x, y - b.y); } comp conj(void) { return comp(x, -y); }
}; const int EPS = 1e-9; template <typename FLOAT_T>
FLOAT_T fabs(const FLOAT_T &x) { return x > 0 ? x : -x;
} template <typename FLOAT_T>
FLOAT_T sin(const FLOAT_T &x, const long double &EPS = fstdlib::EPS) { FLOAT_T res = 0, delt = x; int d = 0; while (fabs(delt) > EPS) { res += delt, ++d; delt *= -x * x / ((2 * d) * (2 * d + 1)); } return res;
} template <typename FLOAT_T>
FLOAT_T cos(const FLOAT_T &x, const long double &EPS = fstdlib::EPS) { FLOAT_T res = 0, delt = 1; int d = 0; while (fabs(delt) > EPS) { res += delt, ++d; delt *= -x * x / ((2 * d) * (2 * d - 1)); } return res;
} const long double PI = std::acos((long double)(-1)); void dft_for_complex(std::vector<comp> &f, int n, int b) { static std::vector<comp> w; w.resize(n); for (int i = 0; i < n; ++i) if (i < rev[i]) std::swap(f[i], f[rev[i]]); for (int i = 2; i <= n; i <<= 1) { w[0] = comp(1, 0), w[1] = comp(cos(2 * PI / i), b * sin(2 * PI / i)); for (int j = 2; j < i / 2; ++j) w[j] = w[j - 1] * w[1]; for (int j = 0; j < n; j += i) for (int k = 0; k < i / 2; ++k) { comp p = f[j + k], q = f[j + k + i / 2] * w[k]; f[j + k] = p + q, f[j + k + i / 2] = p - q; } }
} class arbitrary_module_poly { private: std::vector<int> data; int construct_element(int D, ll x, ll y, ll z) const { x %= mod, y %= mod, z %= mod; return ((ll)D * D * x % mod + (ll)D * y % mod + z) % mod; } public: int mod; arbitrary_module_poly(std::size_t len = std::size_t(0), int module_value = 1e9 + 7) { mod = module_value; data = std::vector<int>(len); } arbitrary_module_poly(const std::vector<int> &b, int module_value = 1e9 + 7) { mod = module_value; data = b; } arbitrary_module_poly(const arbitrary_module_poly &b) { mod = b.mod; data = b.data; } void resize(std::size_t len, const int &val = 0) { data.resize(len, val); } std::size_t size(void) const { return data.size(); } void clear(void) { data.clear(); }
#if __cplusplus >= 201103L void shrink_to_fit(void) { data.shrink_to_fit(); }
#endif int &operator[](std::size_t b) { return data[b]; } const int &operator[](std::size_t b) const { return data[b]; } arbitrary_module_poly operator*(const arbitrary_module_poly &h) const; arbitrary_module_poly operator*=(const arbitrary_module_poly &h); arbitrary_module_poly operator*(const int &h) const; arbitrary_module_poly operator*=(const int &h); arbitrary_module_poly operator+(const arbitrary_module_poly &h) const; arbitrary_module_poly operator+=(const arbitrary_module_poly &h); arbitrary_module_poly operator-(const arbitrary_module_poly &h) const; arbitrary_module_poly operator-=(const arbitrary_module_poly &h); arbitrary_module_poly operator<<(const std::size_t &b) const; arbitrary_module_poly operator<<=(const std::size_t &b); arbitrary_module_poly operator>>(const std::size_t &b) const; arbitrary_module_poly operator>>=(const std::size_t &b); arbitrary_module_poly operator/(const int &h) const; arbitrary_module_poly operator/=(const int &h); arbitrary_module_poly operator==(const arbitrary_module_poly &h) const; arbitrary_module_poly operator!=(const arbitrary_module_poly &h) const; arbitrary_module_poly inv(void) const; arbitrary_module_poly inv(const int &h) const; friend arbitrary_module_poly sqrt(const arbitrary_module_poly &h); friend arbitrary_module_poly log(const arbitrary_module_poly &h);
}; arbitrary_module_poly arbitrary_module_poly::operator*( const arbitrary_module_poly &h) const { int N = 1; while (N < (int)(size() + h.size() - 1)) N <<= 1; std::vector<comp> f(N), g(N), p(N), q(N); const int D = std::sqrt(mod); for (int i = 0; i < (int)size(); ++i) f[i].x = data[i] / D, f[i].y = data[i] % D; for (int i = 0; i < (int)h.size(); ++i) g[i].x = h[i] / D, g[i].y = h[i] % D; rev.resize(N); for (int i = 0; i < N; ++i) rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? N >> 1 : 0); dft_for_complex(f, N, 1), dft_for_complex(g, N, 1); for (int i = 0; i < N; ++i) { p[i] = (f[i] + f[(N - i) % N].conj()) * comp(0.50, 0) * g[i]; q[i] = (f[i] - f[(N - i) % N].conj()) * comp(0, -0.5) * g[i]; } dft_for_complex(p, N, -1), dft_for_complex(q, N, -1); std::vector<int> r(size() + h.size() - 1); for (int i = 0; i < (int)r.size(); ++i) r[i] = construct_element(D, p[i].x / N + 0.5, (p[i].y + q[i].x) / N + 0.5, q[i].y / N + 0.5); return arbitrary_module_poly(r, mod);
} arbitrary_module_poly arbitrary_module_poly::operator*=( const arbitrary_module_poly &h) { return *this = *this * h;
} arbitrary_module_poly arbitrary_module_poly::operator*(const int &h) const { std::vector<int> f(this->data); for (int i = 0; i < (int)f.size(); ++i) f[i] = (ll)f[i] * h % mod; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator*=(const int &h) { for (int i = 0; i < (int)size(); ++i) data[i] = (ll)data[i] * h % mod; return *this;
} arbitrary_module_poly arbitrary_module_poly::operator+( const arbitrary_module_poly &h) const { std::vector<int> f(this->data); if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] + h[i]) % mod; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator+=( const arbitrary_module_poly &h) { if (size() < h.size()) resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) data[i] = (data[i] + h[i]) % mod; return *this;
} arbitrary_module_poly arbitrary_module_poly::operator-( const arbitrary_module_poly &h) const { std::vector<int> f(this->data); if (f.size() < h.size()) f.resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) f[i] = (f[i] + mod - h[i]) % mod; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator-=( const arbitrary_module_poly &h) { if (size() < h.size()) resize(h.size()); for (int i = 0; i < (int)h.size(); ++i) data[i] = (data[i] + mod - h[i]) % mod; return *this;
} arbitrary_module_poly arbitrary_module_poly::operator<<( const std::size_t &b) const { std::vector<int> f(size() + b); for (int i = 0; i < (int)size(); ++i) f[i + b] = data[i]; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator<<=(const std::size_t &b) { return *this = (*this) << b;
} arbitrary_module_poly arbitrary_module_poly::operator>>( const std::size_t &b) const { std::vector<int> f(size() - b); for (int i = 0; i < (int)f.size(); ++i) f[i] = data[i + b]; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator>>=(const std::size_t &b) { return *this = (*this) >> b;
} arbitrary_module_poly arbitrary_module_poly::inv(void) const { int N = 1; while (N < (int)(size() + size() - 1)) N <<= 1; arbitrary_module_poly f(1, mod), g(N, mod), h(*this), f2(1, mod); f[0] = qpow(data[0], mod - 2, mod), h.resize(N), f2[0] = 2; for (int w = 2; w < N; w <<= 1) { g.resize(w); for (int i = 0; i < w; ++i) g[i] = h[i]; f = f * (f * g - f2) * (mod - 1); f.resize(w); } f.resize(size()); return f;
} arbitrary_module_poly arbitrary_module_poly::inv(const int &h) const { arbitrary_module_poly f(*this); f.resize(h); return f.inv();
} arbitrary_module_poly arbitrary_module_poly::operator/(const int &h) const { int inv = qpow(h, mod - 2, mod); std::vector<int> f(this->data); for (int i = 0; i < (int)f.size(); ++i) f[i] = (ll)f[i] * inv % mod; return arbitrary_module_poly(f, mod);
} arbitrary_module_poly arbitrary_module_poly::operator/=(const int &h) { int inv = qpow(h, mod - 2, mod); for (int i = 0; i < (int)size(); ++i) data[i] = (ll)data[i] * inv % mod; return *this;
} arbitrary_module_poly arbitrary_module_poly::operator==( const arbitrary_module_poly &h) const { if (size() != h.size() || mod != h.mod) return 0; for (int i = 0; i < (int)size(); ++i) if (data[i] != h[i]) return 0; return 1;
} arbitrary_module_poly arbitrary_module_poly::operator!=( const arbitrary_module_poly &h) const { if (size() != h.size() || mod != h.mod) return 1; for (int i = 0; i < (int)size(); ++i) if (data[i] != h[i]) return 1; return 0;
} arbitrary_module_poly sqrt(const arbitrary_module_poly &h) { int N = 1; while (N < (int)(h.size() + h.size() - 1)) N <<= 1; arbitrary_module_poly f(1, mod), g(N, mod), d(h); f[0] = modsqrt(h[0], mod), d.resize(N); for (int w = 2; w < N; w <<= 1) { g.resize(w); for (int i = 0; i < w; ++i) g[i] = d[i]; f = (f + f.inv(w) * g) / 2; f.resize(w); } f.resize(h.size()); return f;
} arbitrary_module_poly log(const arbitrary_module_poly &h) { arbitrary_module_poly f(h); for (int i = 1; i < (int)f.size(); ++i) f[i - 1] = (ll)f[i] * i % f.mod; f[f.size() - 1] = 0, f = f * h.inv(), f.resize(h.size()); for (int i = (int)f.size() - 1; i > 0; --i) f[i] = (ll)f[i - 1] * qpow(i, f.mod - 2, f.mod) % f.mod; f[0] = 0; return f;
} typedef arbitrary_module_poly m_poly;
} // namespace fstdlib #endif

int main() { scanf("%d", &n); fact[0] = 1; for (int i = 1; i <= n; ++i) fact[i] = (ll)fact[i - 1] * i % mod; exgcd(fact[n], mod, ifact[n], ifact[0]), ifact[n] = (ifact[n] % mod + mod) % mod; for (int i = n - 1; i >= 0; --i) ifact[i] = (ll)ifact[i + 1] * (i + 1) % mod; poly f(n + 1), g(n + 1); for (int i = 0; i <= n; ++i) g[i] = (i & 1 ? mod - 1ll : 1ll) * ifact[i] % mod, f[i] = (ll)qpow(i, n) * ifact[i] % mod; f *= g, f.resize(n + 1); for (int i = 0; i <= n; ++i) printf("%d ", f[i]); return 0;
}

Method 2. Using exponential generating functions

in a box $i$ The number of options for items and the box is not empty is $[i>0]$ . We can write its exponential generating function as $F(x)=\sum\limits_{i=1}^{+\infty}\dfrac{x^i}{i!} = e^x-1$ . After previous study, we understand $F^k(x)$ That’s it $i$ put a numbered item into $k$ an exponential generating function in a labeled box, $\exp F(x)=\sum\limits_{i=0}^{+\infty}\dfrac{F^i(x)}{i!}$ That’s it $i$ An exponential generating function for placing labeled items into any number of unlabeled boxes (EXP is divided by one $i!$ Box markings removed). There is a lot of "labeled" and "unlabeled" content involved here, so pay attention to the analysis.

Then $\begin{Bmatrix}i\\k\end{Bmatrix}=\dfrac{\left[\dfrac{x^i}{i!}\right]F^k(x)}{k!}$ ， $O(n\log n)$ Just calculate the polynomial powers. Actual usage time ratio $O(n\log n)$ Method 1 is slower.

int main() { scanf("%d%d", &n, &k); poly f(n + 1); fact[0] = 1; for (int i = 1; i <= n; ++i) fact[i] = (ll)fact[i - 1] * i % mod; for (int i = 1; i <= n; ++i) f[i] = qpow(fact[i], mod - 2); f = exp(log(f >> 1) * k) << k, f.resize(n + 1); int inv = qpow(fact[k], mod - 2); for (int i = 0; i <= n; ++i) printf("%lld ", (ll)f[i] * fact[i] % mod * inv % mod); return 0;
}

Sequences (4): Stirling numbers of the second kind

# Stirling Number of the Second Kind

# Recursive

# general formula

# Calculation of Stirling numbers of the second kind in the same row

# Method 1. Directly use the general formula