斯特林的试炼
Reference
第一类Strling数
定义
将 \(n\) 个互异元素放入 \(m\) 个轮换（这 \(m\) 个轮换不分顺序）的方案数称为第一类Strling数，记做 \(\begin{bmatrix}n\newline m\end{bmatrix}\)。
一个轮换实际上就是一个圆排列，例如 \([1,2,3,4]\) 和 \([4,3,2,1]\)、\([3,4,1,2]\)、\([2,3,4,1]\) 都本质相同（可以看作是将它们放在一个可以随意转动的圆上，第一个数和最后一个数相邻）。
例如将 \(4\) 个数放入 \(2\) 个轮换的方案数为 \(11\)，方案为
\[\begin{aligned}
&[1,2,3][4]，[1,2,4][3]，[1,3,4][2]，[2,3,4][1]，\newline
&[1,3,2][4]，[1,4,2][3]，[1,4,3][2]，[2,4,3][1]，\newline
&[1,2][3,4]，[1,3][2,4]，[1,4][2,3]
\end{aligned}
\nonumber
\]
于是 \(\begin{bmatrix}4\newline 2\end{bmatrix}=11\)。
递推式也不难得知，考虑第 \(n\) 个元素是单独作为一个轮换还是和其它元素在同一个轮换，就有
\[\begin{bmatrix}n\newline m\end{bmatrix}=\begin{bmatrix}n-1\newline m-1\end{bmatrix}+(n-1)\begin{bmatrix}n-1\newline m\end{bmatrix}
\nonumber
\]
特别地，我们有 \(\begin{bmatrix}n\newline 0\end{bmatrix} = [n=0]\)，一个特殊值是 \(\begin{bmatrix}n\newline 1\end{bmatrix}=(n-1)!\)。
性质
第一类Strling数在OI中的应用远不如第二类Strling数广泛，但是它仍有一些值得学习的性质。
\[\sum_{k=0}^{n}\begin{bmatrix}n \newline k\end{bmatrix}=n!
\]
考虑用组合意义解释，考虑 \(1 \sim n\) 的一个任意排列 \(\pi\)，连边 \((i,p_i)\) 即可得到一张由一堆环构成的图，而由这张图也可以唯一还原出排列 \(\pi\)，我们得到了一个双射！
枚举环的个数，设有 \(k\) 个，方案数正是第一类Strling数 \(\begin{bmatrix}n \newline k\end{bmatrix}\)。
一个更重要的应用是上升幂与普通幂的转化，记上升幂 \(x^{\overline{n}}=\prod_{k=0}^{n-1}(x+k)\)，下降幂 \(x^{\underline{n}}=\prod_{k=0}^{n-1}(x-k)\) 那么有以下恒等式
\[x^{\overline{n}}=\sum_{k=0}^{n}\begin{bmatrix}n\newline k\end{bmatrix}x^k
\]
\[x^{\underline{n}}=\sum_{k=0}^{n}(-1)^{n-k}\begin{bmatrix}n\newline k\end{bmatrix}x^k
\]
证明考虑归纳。
对于负号的添加，具体数学中给出了一个便于记忆的方法：注意到当 \(x > n > 1\) 时有 \(x^{\overline{n}} > x^n > x^{\underline{n}}\)，而Strling数是非负的，因此将“大的”幂变成“小的”幂时需要添加负号。
计算
第一类Strling数没有实用的通项公式，但我们可以借助多项式科技做到 \(\mathcal{O}(n \log n)\) 求一行或一列的值。所谓“一行”就是上标固定，也就是 \(\begin{bmatrix}n\newline i\end{bmatrix}\)，其中 \(n\) 是定值；“一列”就是指 \(\begin{bmatrix}i\newline n\end{bmatrix}\)。之所以这么叫是因为将它排成标后，上标相同的在同一行，下标相同的在同一列。
第一类Strling数-行
根据之前有的上升幂和普通幂之间的转化，我们有
\[x^{\overline{n}}=\prod_{k=0}^{n-1}(x+k)=\sum_{k=0}^{n}\begin{bmatrix}n\newline k\end{bmatrix}x^k
\nonumber
\]
注意到最后一个等号右边正是 \(\begin{bmatrix}n\newline k\end{bmatrix}\) 的生成函数，于是我们只需要想办法快速求出 \(x^{\overline{n}}=\prod_{k=0}^{n-1}(x+k)\)。显然可以直接暴力分治做到 \(\mathcal{O}(n \log^2 n)\)，但是我们有更优雅的 \(\mathcal{O}(n \log n)\) 做法。
套路地考虑倍增，显然 \(x^{\overline{2n}}=x^{\overline{n}}(x+n)^{\overline{n}}\)，于是设 \(f(x)=x^{\overline{n}}\)，那么 \(x^{\overline{2n}}=f(x)f(x+n)\)。
我们展开 \(f(x+n)\)，做一些常见的变换得到
\[\begin{aligned}
f(x+n)&=\sum_{i=0}^{n}(x+n)^i[x^i]f(x)\newline
&=\sum_{i=0}^{n}\sum_{j=0}^{i}\binom{i}{j}x^j n^{i-j} [x^i]f(x)\newline
&=\sum_{j=0}^{n}\sum_{i=j}^{n}\binom{i}{j}x^j n^{i-j} [x^i]f(x)\newline
&= \sum_{j=0}^{n}x^j\sum_{i=j}^{n}\binom{i}{j}n^{i-j}[x^i]f(x)\newline
&=\sum_{j=0}^{n}x^j\sum_{i=0}^{n-j}\binom{i+j}{j}n^{i}[x^{i+j}]f(x)
\end{aligned}
\nonumber
\]
最后，整理为我们喜欢的形式，得到
\[f(x+n)=\sum_{i=0}^{n}\frac{x^i}{i!}\sum_{j=0}^{n-i}\frac{(i+j)!}{j!}n^j[x^{i+j}]f(x)
\nonumber
\]
记 \(P(x)=\sum_{i}\big(i![x^i]f(x)\big)x^i\)，\(Q(x)=\sum_i \frac{n^i}{i!}x^i\)，就有
\[f(x+n)=\sum_{i=0}^{n}\frac{x^i}{i!}\sum_{j=0}^{n-i}\Big([x^{i+j}]P(x) \Big)\Big([x^j]Q(x)\Big)
\nonumber
\]
记 \(\mathscr{R}P\) 表示反转 \(P\) 的各项系数得到的新多项式，即
\[\mathscr{R}P(x)=\sum_{i=0}^{n}\Big([x^{n-i}]P(x)\Big)x^i
\nonumber
\]
于是
\[f(x+n)=\sum_{i=0}^{n}\frac{x^i}{i!}\sum_{j=0}^{n-i}\Big([x^{n-i-j}]\mathscr{R}P(x)\Big)\Big([x^j]Q(x)\Big)
\nonumber
\]
后面是个卷积的形式，即
\[f(x+n)=\sum_{i=0}^{n}\frac{x^i}{i!}\Big([x^{n-i}]\big(\mathscr{R}P(x) \times Q(x)\big)\Big)
\nonumber
\]
只需做一个多项式乘法。当 \(n\) 是奇数的时候，线性将 \(2\lfloor \frac{n}{2} \rfloor\) 扩展为 \(n\) 即可。
时间复杂度 \(T(n)=T(\frac{n}{2})+\mathcal{O}(n \log n)=\mathcal{O}(n \log n)\)。
code for Luogu P5408
#include <bits/stdc++.h>
using namespace std;
typedef vector<int> poly;
const int N = 5e6 + 10, MOD = 998244353;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace My_poly {
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
int rev[N], iv[N];
inline void make_iv() {
// 处理 [1, N - 1] 的逆元并存到 iv[] 中
iv[1] = 1;
for(int i = 2; i < N; i ++)
iv[i] = Minus(0, 1ll * iv[MOD % i] * (MOD / i) % MOD);
}
inline void make_rev(int n) {
for(int i = 1; i < n; i ++)
rev[i] = ((rev[i >> 1] >> 1) | (i & 1) * (n >> 1));
}
inline void NTT(poly &A, int type) {
// 调用前请确保 rev[] 数组的正确
// type == 1:= DFT
// type == -1:= IDFT
static const int g = 3; // 原根
int n = A.size();
for(int i = 0; i < n; i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int h = 2; h <= n; h <<= 1) {
long long step = ksm(g, (MOD - 1) / h);
if(type == -1) step = ksm(step, MOD - 2);
for(int i = 0; i < n; i += h)
for(int j = i, mul = 1; j < i + (h >> 1); j ++, mul = mul * step % MOD) {
int u = A[j], v = 1ll * A[j + (h >> 1)] * mul % MOD;
A[j] = Plus(u, v), A[j + (h >> 1)] = Minus(u, v);
}
}
if(type == -1) {
long long mul = ksm(n, MOD - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * mul % MOD;
}
}
poly operator + (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Plus(res[i], rhs[i]);
return res;
}
poly operator - (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Minus(res[i], rhs[i]);
return res;
}
poly operator * (const int lhs, const poly &rhs) {
poly res = rhs;
for(int i = 0; i < rhs.size(); i ++)
res[i] = 1ll * res[i] * lhs % MOD;
return res;
}
poly operator * (const poly &lhs, const int rhs) {
poly res = lhs;
for(int i = 0; i < lhs.size(); i ++)
res[i] = 1ll * res[i] * rhs % MOD;
return res;
}
poly operator * (poly A, poly B) {
int h = 1;
while(h <= A.size() + B.size()) h <<= 1;
A.resize(h, 0), B.resize(h, 0);
make_rev(h);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < h; i ++) A[i] = 1ll * A[i] * B[i] % MOD;
NTT(A, -1); return A;
}
inline poly derivative(poly A) {
// 多项式求导
for(int i = 1; i < A.size(); i ++)
A[i - 1] = 1ll * i * A[i] % MOD;
if(!A.empty()) A.pop_back();
return A;
}
inline poly integrate(poly A) {
// 多项式积分
// 使用前请保证调用过make_iv函数或init函数
A.emplace_back(0);
for(int i = A.size() - 1; i >= 1; i --)
A[i] = 1ll * A[i - 1] * iv[i] % MOD;
A[0] = 0; // 不定积分的常数 C
return A;
}
inline poly inv(poly A, int n) {
// 多项式求逆
// 返回模 x^n 意义下 A 的逆元
int h = 1; while(h < n) h <<= 1; A.resize(h, 0);
if(A.empty()) return poly();
poly res(1, ksm(A[0], MOD - 2));
for(int i = 2; i <= h; i <<= 1) {
poly q(A.begin(), A.begin() + i);
res.resize(2 * i, 0), q.resize(2 * i, 0);
make_rev(2 * i), NTT(res, 1), NTT(q, 1);
for(int j = 0; j < 2 * i; j ++) res[j] = 1ll * res[j] * Minus(2, 1ll * res[j] * q[j] % MOD) % MOD;
NTT(res, -1); res.resize(i);
}
res.resize(n); return res;
}
inline poly ln(poly A, int n) {
// 多项式 ln，返回 ln A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 1
A.resize(n, 0);
A = derivative(A) * inv(A, n); A.resize(n);
return integrate(A);
}
inline poly exp(poly A, int n) {
// 多项式 exp，返回exp A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 0
if(n == 1) return poly(1, 1);
int t = (n + 1) >> 1;
poly res = exp(poly(A.begin(), A.begin() + t), t);
res = res * (poly(1, 1) - ln(res, n) + A);
res.resize(n); return res;
}
inline poly power(poly A, int k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
long long val = ksm(A[p], k), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(k * ln(A, n), n);
p = min(1ll * p * k, 1ll * n);
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly power(poly A, string k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
long long mod1 = 0, mod2 = 0;
bool zero = false; // 答案是否为 0
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
for(int i = 0; i < k.size(); i ++) {
if((mod1 * 10 + (k[i] - '0')) * p >= n) {zero = true; break; }
mod1 = (mod1 * 10 + (k[i] - '0')) % MOD;
mod2 = (mod2 * 10 + (k[i] - '0')) % (MOD - 1); // 指数对 \varphi(MOD) 取模
}
if(zero) return poly(n, 0);
long long val = ksm(A[p], mod2), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(mod1 * ln(A, n), n);
p = p * mod1;
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly sqrt1(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时请保证 A[0] = 1，可以稍加修改得到 A[0] 为完全平方数时的算法
int h = 1; while(h < n) h <<= 1;
A.resize(h, 0); poly res(1, 1);
for(int i = 2; i <= h; i <<= 1) {
res = (res * res + poly(A.begin(), A.begin() + i)) * inv(2 * res, i);
res.resize(i);
}
res.resize(n); return res;
}
// 二次剩余相关
long long w;
struct Complex {
int x, y;
Complex(int _x = 0, int _y = 0): x(_x), y(_y) {}
};
Complex operator * (const Complex &lhs, const Complex &rhs)
{return Complex(Plus(1ll * lhs.x * rhs.x % MOD, 1ll * lhs.y * rhs.y % MOD * w % MOD),
Plus(1ll * lhs.x * rhs.y % MOD, 1ll * lhs.y * rhs.x % MOD)); }
inline Complex complex_ksm(Complex A, int b) {
Complex r(1, 0);
for(; b; b >>= 1, A = A * A)
if(b & 1) r = r * A;
return r;
}
long long Cipolla(long long x) {
if (ksm(x,(MOD - 1) >> 1) == MOD - 1) return -1;
while(true) {
long long a = Rand() % MOD;
w = (a * a % MOD + MOD - x) % MOD;
if(ksm(w,(MOD - 1) >> 1) == MOD - 1) {
int res = complex_ksm(Complex(a, 1), (MOD + 1) >> 1).x;
return std::min(res, MOD - res); // 选用首项较小的解
}
}
}
inline poly sqrt(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时不需要保证 A[0] = 1
if(A.empty()) return A;
A.resize(n, 0);
long long val = A[0], mul = ksm(A[0], MOD - 2);
for(int i = 0; i < n; i ++) A[i] = A[i] * mul % MOD;
A = sqrt1(A, n);
val = Cipolla(val);
for(int i = 0; i < n; i ++) A[i] = A[i] * val % MOD;
return A;
}
inline void Rev(poly &A) {
// 反转系数
reverse(A.begin(), A.end());
}
inline pair<poly, poly> Div(poly A, poly B) {
// 多项式除法，返回的 std::pair 中第一维是商，第二维是余数
while(!A.empty() && A.back() == 0) A.pop_back();
while(!B.empty() && B.back() == 0) B.pop_back();
if(A.size() < B.size()) return {poly(1, 0), A};
int n = (int)A.size() - 1, m = (int)B.size() - 1;
Rev(A), Rev(B); poly p = A * inv(B, n - m + 1); p.resize(n - m + 1, 0);
Rev(A), Rev(B), Rev(p); poly r = A - B * p; r.resize(m, 0);
return {p, r};
}
inline void init() {make_iv(); }
} // namespace My_poly
using namespace My_poly;
int main() {
ios::sync_with_stdio(false), cin.tie(0);
init(); // poly 初始化
return 0;
}
第一类Strling数-列
从组合意义上来讲，单个轮换的EGF就是
\[\sum_{i=1}^{+\infty}\frac{(i-1)!x^i}{i!}=\sum_{i=1}^{+\infty}\frac{x^i}{i}
\nonumber
\]
其 \(k\) 次幂就是 \(k\) 个轮换的EGF，也就是
\[\sum_{i=0}^{+\infty}\begin{bmatrix}i\newline k\end{bmatrix}\frac{x^i}{i!}
\nonumber
\]
多项式快速幂碾过去，时间复杂度 \(\mathcal{O}(n \log n)\)。
从代数意义上讲，考虑将 \((x+1)^k\) 二项式定理展开后的下降幂用第一类Strling数展开，立即得到：
\[(x+1)^k=\sum_{i=0}^{+\infty}\frac{k^{\underline{i}}}{i!}x^i=\sum_{i=0}^{+\infty}\frac{x^i}{i!}\sum_{j=0}^{i}(-1)^{i-j}\begin{bmatrix}i\newline j\end{bmatrix}k^j
\nonumber
\]
整理得到
\[(x+1)^k=\sum_{i=0}^{+\infty}k^i\sum_{j=i}^{+\infty}(-1)^{j-i}\begin{bmatrix}j\newline i\end{bmatrix}\frac{x^j}{j!}
\nonumber
\]
另一方面，将 \((x+1)^k\) 先 \(\ln\) 再 \(\exp\)，得到
\[(x+1)^k=\exp\big(k\ln(x+1)\big)=\sum_{i=0}^{+\infty}k^i\frac{\big(\ln(x+1)\big)^i}{i!}
\nonumber
\]
将两个式合起来，得到
\[\sum_{i=0}^{+\infty}k^i\sum_{j=i}^{+\infty}(-1)^{j-i}\begin{bmatrix}j\newline i\end{bmatrix}\frac{x^j}{j!}=\sum_{i=0}^{+\infty}k^i\frac{\big(\ln(x+1)\big)^i}{i!}
\nonumber
\nonumber
\]
于是就有
\[\sum_{j=i}^{+\infty}(-1)^{j-i}\begin{bmatrix}j\newline i\end{bmatrix}\frac{x^j}{j!}=\frac{\big(\ln(x+1)\big)^i}{i!}
\]
做一个多项式 \(\ln\) 和多项式快速幂即可。
实际上可以进一步化简，得到与组合意义相同的结果：
\[\sum_{j=0}^{+\infty}(-1)^j\begin{bmatrix}j\newline i\end{bmatrix}\frac{x^j}{j!}=\frac{\Big(\ln\left(\frac{1}{1+x}\right)\Big)^i}{i!}
\nonumber
\]
用 \(-x\) 替换 \(x\) 就有
\[\sum_{j=0}^{+\infty}\begin{bmatrix}j\newline i\end{bmatrix}\frac{x^j}{j!}=\frac{\Big(\ln\left(\frac{1}{1-x}\right)\Big)^i}{i!}
\nonumber
\]
而注意到 \(\ln\left(\frac{1}{1-x}\right)\) 在 \(0\) 处的Taylor展开正是
\[\sum_{j=1}^{+\infty}\frac{x^j}{j}
\nonumber
\]
得到了相同的结果。
code for Luogu P5409
#include <bits/stdc++.h>
using namespace std;
typedef vector<int> poly;
const int N = 5e6 + 10, MOD = 167772161;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace My_poly {
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
int rev[N], iv[N];
inline void make_iv() {
// 处理 [1, N - 1] 的逆元并存到 iv[] 中
iv[1] = 1;
for(int i = 2; i < N; i ++)
iv[i] = Minus(0, 1ll * iv[MOD % i] * (MOD / i) % MOD);
}
inline void make_rev(int n) {
for(int i = 1; i < n; i ++)
rev[i] = ((rev[i >> 1] >> 1) | (i & 1) * (n >> 1));
}
inline void NTT(poly &A, int type) {
// 调用前请确保 rev[] 数组的正确
// type == 1:= DFT
// type == -1:= IDFT
static const int g = 3; // 原根
int n = A.size();
for(int i = 0; i < n; i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int h = 2; h <= n; h <<= 1) {
long long step = ksm(g, (MOD - 1) / h);
if(type == -1) step = ksm(step, MOD - 2);
for(int i = 0; i < n; i += h)
for(int j = i, mul = 1; j < i + (h >> 1); j ++, mul = mul * step % MOD) {
int u = A[j], v = 1ll * A[j + (h >> 1)] * mul % MOD;
A[j] = Plus(u, v), A[j + (h >> 1)] = Minus(u, v);
}
}
if(type == -1) {
long long mul = ksm(n, MOD - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * mul % MOD;
}
}
poly operator + (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Plus(res[i], rhs[i]);
return res;
}
poly operator - (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Minus(res[i], rhs[i]);
return res;
}
poly operator * (const int lhs, const poly &rhs) {
poly res = rhs;
for(int i = 0; i < rhs.size(); i ++)
res[i] = 1ll * res[i] * lhs % MOD;
return res;
}
poly operator * (const poly &lhs, const int rhs) {
poly res = lhs;
for(int i = 0; i < lhs.size(); i ++)
res[i] = 1ll * res[i] * rhs % MOD;
return res;
}
poly operator * (poly A, poly B) {
int h = 1;
while(h <= A.size() + B.size()) h <<= 1;
A.resize(h, 0), B.resize(h, 0);
make_rev(h);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < h; i ++) A[i] = 1ll * A[i] * B[i] % MOD;
NTT(A, -1); return A;
}
inline poly derivative(poly A) {
// 多项式求导
for(int i = 1; i < A.size(); i ++)
A[i - 1] = 1ll * i * A[i] % MOD;
if(!A.empty()) A.pop_back();
return A;
}
inline poly integrate(poly A) {
// 多项式积分
// 使用前请保证调用过make_iv函数或init函数
A.emplace_back(0);
for(int i = A.size() - 1; i >= 1; i --)
A[i] = 1ll * A[i - 1] * iv[i] % MOD;
A[0] = 0; // 不定积分的常数 C
return A;
}
inline poly inv(poly A, int n) {
// 多项式求逆
// 返回模 x^n 意义下 A 的逆元
int h = 1; while(h < n) h <<= 1; A.resize(h, 0);
if(A.empty()) return poly();
poly res(1, ksm(A[0], MOD - 2));
for(int i = 2; i <= h; i <<= 1) {
poly q(A.begin(), A.begin() + i);
res.resize(2 * i, 0), q.resize(2 * i, 0);
make_rev(2 * i), NTT(res, 1), NTT(q, 1);
for(int j = 0; j < 2 * i; j ++) res[j] = 1ll * res[j] * Minus(2, 1ll * res[j] * q[j] % MOD) % MOD;
NTT(res, -1); res.resize(i);
}
res.resize(n); return res;
}
inline poly ln(poly A, int n) {
// 多项式 ln，返回 ln A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 1
A.resize(n, 0);
A = derivative(A) * inv(A, n); A.resize(n);
return integrate(A);
}
inline poly exp(poly A, int n) {
// 多项式 exp，返回exp A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 0
if(n == 1) return poly(1, 1);
int t = (n + 1) >> 1;
poly res = exp(poly(A.begin(), A.begin() + t), t);
res = res * (poly(1, 1) - ln(res, n) + A);
res.resize(n); return res;
}
inline poly power(poly A, int k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
long long val = ksm(A[p], k), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(k * ln(A, n), n);
p = min(1ll * p * k, 1ll * n);
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly power(poly A, string k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
long long mod1 = 0, mod2 = 0;
bool zero = false; // 答案是否为 0
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
for(int i = 0; i < k.size(); i ++) {
if((mod1 * 10 + (k[i] - '0')) * p >= n) {zero = true; break; }
mod1 = (mod1 * 10 + (k[i] - '0')) % MOD;
mod2 = (mod2 * 10 + (k[i] - '0')) % (MOD - 1); // 指数对 \varphi(MOD) 取模
}
if(zero) return poly(n, 0);
long long val = ksm(A[p], mod2), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(mod1 * ln(A, n), n);
p = p * mod1;
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly sqrt1(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时请保证 A[0] = 1，可以稍加修改得到 A[0] 为完全平方数时的算法
int h = 1; while(h < n) h <<= 1;
A.resize(h, 0); poly res(1, 1);
for(int i = 2; i <= h; i <<= 1) {
res = (res * res + poly(A.begin(), A.begin() + i)) * inv(2 * res, i);
res.resize(i);
}
res.resize(n); return res;
}
// 二次剩余相关
long long w;
struct Complex {
int x, y;
Complex(int _x = 0, int _y = 0): x(_x), y(_y) {}
};
Complex operator * (const Complex &lhs, const Complex &rhs)
{return Complex(Plus(1ll * lhs.x * rhs.x % MOD, 1ll * lhs.y * rhs.y % MOD * w % MOD),
Plus(1ll * lhs.x * rhs.y % MOD, 1ll * lhs.y * rhs.x % MOD)); }
inline Complex complex_ksm(Complex A, int b) {
Complex r(1, 0);
for(; b; b >>= 1, A = A * A)
if(b & 1) r = r * A;
return r;
}
long long Cipolla(long long x) {
if (ksm(x,(MOD - 1) >> 1) == MOD - 1) return -1;
while(true) {
long long a = Rand() % MOD;
w = (a * a % MOD + MOD - x) % MOD;
if(ksm(w,(MOD - 1) >> 1) == MOD - 1) {
int res = complex_ksm(Complex(a, 1), (MOD + 1) >> 1).x;
return std::min(res, MOD - res); // 选用首项较小的解
}
}
}
inline poly sqrt(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时不需要保证 A[0] = 1
if(A.empty()) return A;
A.resize(n, 0);
long long val = A[0], mul = ksm(A[0], MOD - 2);
for(int i = 0; i < n; i ++) A[i] = A[i] * mul % MOD;
A = sqrt1(A, n);
val = Cipolla(val);
for(int i = 0; i < n; i ++) A[i] = A[i] * val % MOD;
return A;
}
inline void Rev(poly &A) {
// 反转系数
reverse(A.begin(), A.end());
}
inline pair<poly, poly> Div(poly A, poly B) {
// 多项式除法，返回的 std::pair 中第一维是商，第二维是余数
while(!A.empty() && A.back() == 0) A.pop_back();
while(!B.empty() && B.back() == 0) B.pop_back();
if(A.size() < B.size()) return {poly(1, 0), A};
int n = (int)A.size() - 1, m = (int)B.size() - 1;
Rev(A), Rev(B); poly p = A * inv(B, n - m + 1); p.resize(n - m + 1, 0);
Rev(A), Rev(B), Rev(p); poly r = A - B * p; r.resize(m, 0);
return {p, r};
}
inline void init() {make_iv(); }
} // namespace My_poly
using namespace My_poly;
int fac[N], ifac[N];
int main() {
ios::sync_with_stdio(false), cin.tie(0);
init(); // poly 初始化
fac[0] = ifac[0] = 1;
for(int i = 1; i < N; i ++) {
fac[i] = 1ll * fac[i - 1] * i % MOD;
ifac[i] = 1ll * ifac[i - 1] * iv[i] % MOD;
}
int n, k; cin >> n >> k;
poly f(n + 1, 0);
for(int i = 1; i <= n; i ++) f[i] = iv[i];
f = power(f, k, n + 1);
for(int i = 0; i <= n; i ++)
cout << 1ll * f[i] * ifac[k] % MOD * fac[i] % MOD << " \n"[i == n];
return 0;
}
第二类Strling数
定义
将 \(n\) 个两两不同的元素放入 \(m\) 个互不区分的非空集合中的方案数称为第二类Strling数，记作 \(\begin{Bmatrix}n\newline m\end{Bmatrix}\)。
例如将 \(4\) 个数放入 \(2\) 个集合中有 \(7\) 种方法：
\[\begin{aligned}
&\{1,2,3\}\{4\}，\{1,2,4\}\{3\}，\{1,3,4\}\{2\}，\{2,3,4\}\{1\}\newline
&\{1,2\}\{3,4\}，\{1,3\}\{2,4\}，\{1,4\}\{2,3\}
\end{aligned}
\nonumber
\]
显然我们应当有
\[\begin{Bmatrix}n\newline m\end{Bmatrix} \le \begin{bmatrix}n\newline m\end{bmatrix}
\]
因为一个集合总可以对应多个轮换。
同样用组合意义求递推公式，考虑最后一个元素是单独作为一个集合还是被加入之前的集合中，得到
\[\begin{Bmatrix}n\newline m\end{Bmatrix}=\begin{Bmatrix}n-1\newline m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\newline m\end{Bmatrix}
\]
特别地，有 \(\begin{Bmatrix}n\newline 0\end{Bmatrix}=[n=0]\)。
性质
和第一类Strling数一样，第二类Strling数也可以将普通幂和上升幂、下降幂互相转换：
\[x^n=\sum_{k=0}^{n}\begin{Bmatrix}n\newline k\end{Bmatrix}x^{\underline{k}}
\]
\[x^n=\sum_{k=0}^{n}(-1)^{n-k}\begin{Bmatrix}n\newline k\end{Bmatrix}x^{\overline{k}}
\]
计算
通项公式
第二类Strling数有一个好用的通项公式：
\[\begin{Bmatrix}n\newline m\end{Bmatrix}=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}
\]
证明考虑二项式反演，设将 \(n\) 个两两不同的元素划分到 \(i\) 个可以是空集的、两两不同的集合中的方案数为 \(G_i\)，划分到 \(i\) 个非空的、两两不同的集合中的方案数为 \(F_i\)，那么就有
\[G_i=\sum_{j=0}^{i}\binom{i}{j}F_j=i^n
\nonumber
\]
立即得到
\[F_i=\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}G_j=\sum_{j=0}^{i}(-1)^{i-j}\binom{i}{j}j^n
\nonumber
\]
考虑第二类Strling数和 \(F\) 的区别只有集合之间是否区分，于是 \(\begin{Bmatrix}n\newline i\end{Bmatrix}=\dfrac{F_i}{i!}\)，整理得
\[\begin{Bmatrix}n\newline m\end{Bmatrix}=\sum_{i=0}^{m}\frac{(-1)^{m-i}i^n}{i!(m-i)!}
\nonumber
\]
证毕。
第二类Strling数-行
根据通项公式做卷积即可，时间复杂度 \(\mathcal{O}(n \log n)\)。
code for Luogu P5395
#include <bits/stdc++.h>
using namespace std;
typedef vector<int> poly;
const int N = 5e6 + 10, MOD = 167772161;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace My_poly {
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
int rev[N], iv[N];
inline void make_iv() {
// 处理 [1, N - 1] 的逆元并存到 iv[] 中
iv[1] = 1;
for(int i = 2; i < N; i ++)
iv[i] = Minus(0, 1ll * iv[MOD % i] * (MOD / i) % MOD);
}
inline void make_rev(int n) {
for(int i = 1; i < n; i ++)
rev[i] = ((rev[i >> 1] >> 1) | (i & 1) * (n >> 1));
}
inline void NTT(poly &A, int type) {
// 调用前请确保 rev[] 数组的正确
// type == 1:= DFT
// type == -1:= IDFT
static const int g = 3; // 原根
int n = A.size();
for(int i = 0; i < n; i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int h = 2; h <= n; h <<= 1) {
long long step = ksm(g, (MOD - 1) / h);
if(type == -1) step = ksm(step, MOD - 2);
for(int i = 0; i < n; i += h)
for(int j = i, mul = 1; j < i + (h >> 1); j ++, mul = mul * step % MOD) {
int u = A[j], v = 1ll * A[j + (h >> 1)] * mul % MOD;
A[j] = Plus(u, v), A[j + (h >> 1)] = Minus(u, v);
}
}
if(type == -1) {
long long mul = ksm(n, MOD - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * mul % MOD;
}
}
poly operator + (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Plus(res[i], rhs[i]);
return res;
}
poly operator - (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Minus(res[i], rhs[i]);
return res;
}
poly operator * (const int lhs, const poly &rhs) {
poly res = rhs;
for(int i = 0; i < rhs.size(); i ++)
res[i] = 1ll * res[i] * lhs % MOD;
return res;
}
poly operator * (const poly &lhs, const int rhs) {
poly res = lhs;
for(int i = 0; i < lhs.size(); i ++)
res[i] = 1ll * res[i] * rhs % MOD;
return res;
}
poly operator * (poly A, poly B) {
int h = 1;
while(h <= A.size() + B.size()) h <<= 1;
A.resize(h, 0), B.resize(h, 0);
make_rev(h);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < h; i ++) A[i] = 1ll * A[i] * B[i] % MOD;
NTT(A, -1); return A;
}
inline poly derivative(poly A) {
// 多项式求导
for(int i = 1; i < A.size(); i ++)
A[i - 1] = 1ll * i * A[i] % MOD;
if(!A.empty()) A.pop_back();
return A;
}
inline poly integrate(poly A) {
// 多项式积分
// 使用前请保证调用过make_iv函数或init函数
A.emplace_back(0);
for(int i = A.size() - 1; i >= 1; i --)
A[i] = 1ll * A[i - 1] * iv[i] % MOD;
A[0] = 0; // 不定积分的常数 C
return A;
}
inline poly inv(poly A, int n) {
// 多项式求逆
// 返回模 x^n 意义下 A 的逆元
int h = 1; while(h < n) h <<= 1; A.resize(h, 0);
if(A.empty()) return poly();
poly res(1, ksm(A[0], MOD - 2));
for(int i = 2; i <= h; i <<= 1) {
poly q(A.begin(), A.begin() + i);
res.resize(2 * i, 0), q.resize(2 * i, 0);
make_rev(2 * i), NTT(res, 1), NTT(q, 1);
for(int j = 0; j < 2 * i; j ++) res[j] = 1ll * res[j] * Minus(2, 1ll * res[j] * q[j] % MOD) % MOD;
NTT(res, -1); res.resize(i);
}
res.resize(n); return res;
}
inline poly ln(poly A, int n) {
// 多项式 ln，返回 ln A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 1
A.resize(n, 0);
A = derivative(A) * inv(A, n); A.resize(n);
return integrate(A);
}
inline poly exp(poly A, int n) {
// 多项式 exp，返回exp A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 0
if(n == 1) return poly(1, 1);
int t = (n + 1) >> 1;
poly res = exp(poly(A.begin(), A.begin() + t), t);
res = res * (poly(1, 1) - ln(res, n) + A);
res.resize(n); return res;
}
inline poly power(poly A, int k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
long long val = ksm(A[p], k), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(k * ln(A, n), n);
p = min(1ll * p * k, 1ll * n);
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly power(poly A, string k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
long long mod1 = 0, mod2 = 0;
bool zero = false; // 答案是否为 0
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
for(int i = 0; i < k.size(); i ++) {
if((mod1 * 10 + (k[i] - '0')) * p >= n) {zero = true; break; }
mod1 = (mod1 * 10 + (k[i] - '0')) % MOD;
mod2 = (mod2 * 10 + (k[i] - '0')) % (MOD - 1); // 指数对 \varphi(MOD) 取模
}
if(zero) return poly(n, 0);
long long val = ksm(A[p], mod2), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(mod1 * ln(A, n), n);
p = p * mod1;
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly sqrt1(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时请保证 A[0] = 1，可以稍加修改得到 A[0] 为完全平方数时的算法
int h = 1; while(h < n) h <<= 1;
A.resize(h, 0); poly res(1, 1);
for(int i = 2; i <= h; i <<= 1) {
res = (res * res + poly(A.begin(), A.begin() + i)) * inv(2 * res, i);
res.resize(i);
}
res.resize(n); return res;
}
// 二次剩余相关
long long w;
struct Complex {
int x, y;
Complex(int _x = 0, int _y = 0): x(_x), y(_y) {}
};
Complex operator * (const Complex &lhs, const Complex &rhs)
{return Complex(Plus(1ll * lhs.x * rhs.x % MOD, 1ll * lhs.y * rhs.y % MOD * w % MOD),
Plus(1ll * lhs.x * rhs.y % MOD, 1ll * lhs.y * rhs.x % MOD)); }
inline Complex complex_ksm(Complex A, int b) {
Complex r(1, 0);
for(; b; b >>= 1, A = A * A)
if(b & 1) r = r * A;
return r;
}
long long Cipolla(long long x) {
if (ksm(x,(MOD - 1) >> 1) == MOD - 1) return -1;
while(true) {
long long a = Rand() % MOD;
w = (a * a % MOD + MOD - x) % MOD;
if(ksm(w,(MOD - 1) >> 1) == MOD - 1) {
int res = complex_ksm(Complex(a, 1), (MOD + 1) >> 1).x;
return std::min(res, MOD - res); // 选用首项较小的解
}
}
}
inline poly sqrt(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时不需要保证 A[0] = 1
if(A.empty()) return A;
A.resize(n, 0);
long long val = A[0], mul = ksm(A[0], MOD - 2);
for(int i = 0; i < n; i ++) A[i] = A[i] * mul % MOD;
A = sqrt1(A, n);
val = Cipolla(val);
for(int i = 0; i < n; i ++) A[i] = A[i] * val % MOD;
return A;
}
inline void Rev(poly &A) {
// 反转系数
reverse(A.begin(), A.end());
}
inline pair<poly, poly> Div(poly A, poly B) {
// 多项式除法，返回的 std::pair 中第一维是商，第二维是余数
while(!A.empty() && A.back() == 0) A.pop_back();
while(!B.empty() && B.back() == 0) B.pop_back();
if(A.size() < B.size()) return {poly(1, 0), A};
int n = (int)A.size() - 1, m = (int)B.size() - 1;
Rev(A), Rev(B); poly p = A * inv(B, n - m + 1); p.resize(n - m + 1, 0);
Rev(A), Rev(B), Rev(p); poly r = A - B * p; r.resize(m, 0);
return {p, r};
}
inline void init() {make_iv(); }
} // namespace My_poly
using namespace My_poly;
int fac[N], ifac[N];
int main() {
ios::sync_with_stdio(false), cin.tie(0);
init(); // poly 初始化
fac[0] = ifac[0] = 1;
for(int i = 1; i < N; i ++) {
fac[i] = 1ll * fac[i - 1] * i % MOD;
ifac[i] = 1ll * ifac[i - 1] * iv[i] % MOD;
}
int n; cin >> n;
poly f(n + 1, 0), g(n + 1, 0);
for(int i = 0; i <= n; i ++) {
f[i] = 1ll * ksm(i, n) * ifac[i] % MOD;
g[i] = ifac[i]; if(i & 1) g[i] = Minus(0, g[i]);
}
f = f * g;
for(int i = 0; i <= n; i ++)
cout << f[i] << " \n"[i == n];
return 0;
}
第二类Strling数-列
组合意义，一个非空集合的EGF为
\[\sum_{i=1}^{+\infty}\frac{x^i}{i!}=e^x-1
\nonumber
\]
将这个EGF \(k\) 次幂就可以得到将若干个不同元素放入 \(k\) 个有标号集合的EGF，再除以 \(k!\) 就可以得到同一列第二类 Strling数的EGF了。
时间复杂度 \(\mathcal{O}(n \log n)\)。
code for Luogu P5396
#include <bits/stdc++.h>
using namespace std;
typedef vector<int> poly;
const int N = 5e6 + 10, MOD = 167772161;
inline int Plus(int a, int b) {return a + b >= MOD ? a + b - MOD : a + b; }
inline int Minus(int a, int b) {return a - b < 0 ? a - b + MOD : a - b; }
inline int ksm(long long a, int b) {
long long r = 1;
for(; b; b >>= 1, a = a * a % MOD)
if(b & 1) r = r * a % MOD;
return r;
}
namespace My_poly {
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());
int rev[N], iv[N];
inline void make_iv() {
// 处理 [1, N - 1] 的逆元并存到 iv[] 中
iv[1] = 1;
for(int i = 2; i < N; i ++)
iv[i] = Minus(0, 1ll * iv[MOD % i] * (MOD / i) % MOD);
}
inline void make_rev(int n) {
for(int i = 1; i < n; i ++)
rev[i] = ((rev[i >> 1] >> 1) | (i & 1) * (n >> 1));
}
inline void NTT(poly &A, int type) {
// 调用前请确保 rev[] 数组的正确
// type == 1:= DFT
// type == -1:= IDFT
static const int g = 3; // 原根
int n = A.size();
for(int i = 0; i < n; i ++)
if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int h = 2; h <= n; h <<= 1) {
long long step = ksm(g, (MOD - 1) / h);
if(type == -1) step = ksm(step, MOD - 2);
for(int i = 0; i < n; i += h)
for(int j = i, mul = 1; j < i + (h >> 1); j ++, mul = mul * step % MOD) {
int u = A[j], v = 1ll * A[j + (h >> 1)] * mul % MOD;
A[j] = Plus(u, v), A[j + (h >> 1)] = Minus(u, v);
}
}
if(type == -1) {
long long mul = ksm(n, MOD - 2);
for(int i = 0; i < n; i ++)
A[i] = A[i] * mul % MOD;
}
}
poly operator + (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Plus(res[i], rhs[i]);
return res;
}
poly operator - (const poly &lhs, const poly &rhs) {
poly res = lhs; res.resize(max(lhs.size(), rhs.size()), 0);
for(int i = 0; i < rhs.size(); i ++) res[i] = Minus(res[i], rhs[i]);
return res;
}
poly operator * (const int lhs, const poly &rhs) {
poly res = rhs;
for(int i = 0; i < rhs.size(); i ++)
res[i] = 1ll * res[i] * lhs % MOD;
return res;
}
poly operator * (const poly &lhs, const int rhs) {
poly res = lhs;
for(int i = 0; i < lhs.size(); i ++)
res[i] = 1ll * res[i] * rhs % MOD;
return res;
}
poly operator * (poly A, poly B) {
int h = 1;
while(h <= A.size() + B.size()) h <<= 1;
A.resize(h, 0), B.resize(h, 0);
make_rev(h);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < h; i ++) A[i] = 1ll * A[i] * B[i] % MOD;
NTT(A, -1); return A;
}
inline poly derivative(poly A) {
// 多项式求导
for(int i = 1; i < A.size(); i ++)
A[i - 1] = 1ll * i * A[i] % MOD;
if(!A.empty()) A.pop_back();
return A;
}
inline poly integrate(poly A) {
// 多项式积分
// 使用前请保证调用过make_iv函数或init函数
A.emplace_back(0);
for(int i = A.size() - 1; i >= 1; i --)
A[i] = 1ll * A[i - 1] * iv[i] % MOD;
A[0] = 0; // 不定积分的常数 C
return A;
}
inline poly inv(poly A, int n) {
// 多项式求逆
// 返回模 x^n 意义下 A 的逆元
int h = 1; while(h < n) h <<= 1; A.resize(h, 0);
if(A.empty()) return poly();
poly res(1, ksm(A[0], MOD - 2));
for(int i = 2; i <= h; i <<= 1) {
poly q(A.begin(), A.begin() + i);
res.resize(2 * i, 0), q.resize(2 * i, 0);
make_rev(2 * i), NTT(res, 1), NTT(q, 1);
for(int j = 0; j < 2 * i; j ++) res[j] = 1ll * res[j] * Minus(2, 1ll * res[j] * q[j] % MOD) % MOD;
NTT(res, -1); res.resize(i);
}
res.resize(n); return res;
}
inline poly ln(poly A, int n) {
// 多项式 ln，返回 ln A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 1
A.resize(n, 0);
A = derivative(A) * inv(A, n); A.resize(n);
return integrate(A);
}
inline poly exp(poly A, int n) {
// 多项式 exp，返回exp A 在模 x^n 意义下的结果
// 调用前请保证 A[0] = 0
if(n == 1) return poly(1, 1);
int t = (n + 1) >> 1;
poly res = exp(poly(A.begin(), A.begin() + t), t);
res = res * (poly(1, 1) - ln(res, n) + A);
res.resize(n); return res;
}
inline poly power(poly A, int k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
long long val = ksm(A[p], k), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(k * ln(A, n), n);
p = min(1ll * p * k, 1ll * n);
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly power(poly A, string k, int n) {
// 多项式快速幂，返回 A^k 在模 x^n 意义下的结果
// 不要求 A[0] = 1
long long mod1 = 0, mod2 = 0;
bool zero = false; // 答案是否为 0
int p = -1; A.resize(n, 0);
for(int i = 0; i < n; i ++)
if(A[i] != 0) {p = i; break; }
if(p == -1) return A; // A = 0
for(int i = 0; i < k.size(); i ++) {
if((mod1 * 10 + (k[i] - '0')) * p >= n) {zero = true; break; }
mod1 = (mod1 * 10 + (k[i] - '0')) % MOD;
mod2 = (mod2 * 10 + (k[i] - '0')) % (MOD - 1); // 指数对 \varphi(MOD) 取模
}
if(zero) return poly(n, 0);
long long val = ksm(A[p], mod2), mul = ksm(A[p], MOD - 2);
for(int i = p; i < n; i ++) A[i - p] = A[i] * mul % MOD;
for(int i = n - p; i < n; i ++) A[i] = 0;
A = exp(mod1 * ln(A, n), n);
p = p * mod1;
for(int i = n - 1; i >= p; i --) A[i] = A[i - p] * val % MOD;
for(int i = 0; i < p; i ++) A[i] = 0;
return A;
}
inline poly sqrt1(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时请保证 A[0] = 1，可以稍加修改得到 A[0] 为完全平方数时的算法
int h = 1; while(h < n) h <<= 1;
A.resize(h, 0); poly res(1, 1);
for(int i = 2; i <= h; i <<= 1) {
res = (res * res + poly(A.begin(), A.begin() + i)) * inv(2 * res, i);
res.resize(i);
}
res.resize(n); return res;
}
// 二次剩余相关
long long w;
struct Complex {
int x, y;
Complex(int _x = 0, int _y = 0): x(_x), y(_y) {}
};
Complex operator * (const Complex &lhs, const Complex &rhs)
{return Complex(Plus(1ll * lhs.x * rhs.x % MOD, 1ll * lhs.y * rhs.y % MOD * w % MOD),
Plus(1ll * lhs.x * rhs.y % MOD, 1ll * lhs.y * rhs.x % MOD)); }
inline Complex complex_ksm(Complex A, int b) {
Complex r(1, 0);
for(; b; b >>= 1, A = A * A)
if(b & 1) r = r * A;
return r;
}
long long Cipolla(long long x) {
if (ksm(x,(MOD - 1) >> 1) == MOD - 1) return -1;
while(true) {
long long a = Rand() % MOD;
w = (a * a % MOD + MOD - x) % MOD;
if(ksm(w,(MOD - 1) >> 1) == MOD - 1) {
int res = complex_ksm(Complex(a, 1), (MOD + 1) >> 1).x;
return std::min(res, MOD - res); // 选用首项较小的解
}
}
}
inline poly sqrt(poly A, int n) {
// 多项式开根，返回 \sqrt{A} 在模 x^n 意义下的结果
// 调用时不需要保证 A[0] = 1
if(A.empty()) return A;
A.resize(n, 0);
long long val = A[0], mul = ksm(A[0], MOD - 2);
for(int i = 0; i < n; i ++) A[i] = A[i] * mul % MOD;
A = sqrt1(A, n);
val = Cipolla(val);
for(int i = 0; i < n; i ++) A[i] = A[i] * val % MOD;
return A;
}
inline void Rev(poly &A) {
// 反转系数
reverse(A.begin(), A.end());
}
inline pair<poly, poly> Div(poly A, poly B) {
// 多项式除法，返回的 std::pair 中第一维是商，第二维是余数
while(!A.empty() && A.back() == 0) A.pop_back();
while(!B.empty() && B.back() == 0) B.pop_back();
if(A.size() < B.size()) return {poly(1, 0), A};
int n = (int)A.size() - 1, m = (int)B.size() - 1;
Rev(A), Rev(B); poly p = A * inv(B, n - m + 1); p.resize(n - m + 1, 0);
Rev(A), Rev(B), Rev(p); poly r = A - B * p; r.resize(m, 0);
return {p, r};
}
inline void init() {make_iv(); }
} // namespace My_poly
using namespace My_poly;
int fac[N], ifac[N];
int main() {
ios::sync_with_stdio(false), cin.tie(0);
init(); // poly 初始化
fac[0] = ifac[0] = 1;
for(int i = 1; i < N; i ++) {
fac[i] = 1ll * fac[i - 1] * i % MOD;
ifac[i] = 1ll * ifac[i - 1] * iv[i] % MOD;
}
int n, k; cin >> n >> k;
poly f(n + 1, 0);
for(int i = 1; i <= n; i ++)
f[i] = ifac[i];
f = power(f, k, n + 1);
for(int i = 0; i <= n; i ++)
cout << 1ll * f[i] * ifac[k] % MOD * fac[i] % MOD << " \n"[i == n];
return 0;
}
浙公网安备 33010602011771号