「洛谷P5408」第一类斯特林数·行 题解
题意
求 $\left[\begin{matrix}n\\i\end{matrix}\right],i=0,1,\dots,n$
其中 $\left[\begin{matrix}n\\i\end{matrix}\right]$ 为第一类斯特林数，组合意义为将 $n$ 个不同元素放入$i$个相同圆排列（非空）的方案数，有递推式 $\left[\begin{matrix}n\\m\end{matrix}\right]=\left[\begin{matrix}n-1\\m-1\end{matrix}\right]+\left[\begin{matrix}n-1\\m\end{matrix}\right]\times(n-1)$
推导
一个结论：
$$\boxed{\sum\limits^n_{i=0}\left[\begin{matrix}n\\i\end{matrix}\right]x^i=x^{\overline{n}}}$$
$x^{\overline{n}}$ 为 $x$ 的 $n$ 次上升幂，等于$\prod\limits^{n-1}_{i=0}(x+i)$
考虑用数学归纳法证明：
$n=1$ 时该结论显然成立
证明为 $n$ 时成立则为 $n+1$ 时也成立：
$$\begin{aligned}x^{\overline{n+1}}&=x^{\overline{n}}\times(x+n)\\&=(x+n)\times\sum\limits^n_{i=0}\left[\begin{matrix}n\\i\end{matrix}\right]x^i\\&=\sum\limits^{n+1}_{i=1}\left[\begin{matrix}n\\i-1\end{matrix}\right]x^i+\sum\limits^n_{i=0}\left(n\times\left[\begin{matrix}n\\i\end{matrix}\right]x^i\right)\\&=\sum\limits^{n+1}_{i=0}\left(\left[\begin{matrix}n\\i-1\end{matrix}\right]+n\times\left[\begin{matrix}n\\i\end{matrix}\right]\right)x^i\\&=\sum\limits^{n+1}_{i=0}\left[\begin{matrix}n+1\\i\end{matrix}\right]x^i\end{aligned}$$
于是原命题成立
直接求上升幂是 $O(n^2)$，考虑如何优化：
假设 $n$ 是偶数：
$$x^{\overline{n}}=x^{\overline{\frac{n}{2}}}\times(x+\frac{n}{2})^{\overline{\frac{n}{2}}}$$
如果我们已经得到了前半部分对应的多项式 $f(x)=x^{\overline{\frac{n}{2}}}=\sum\limits^{\frac{n}{2}}_{i=0}a_ix^i$
考虑如何通过前半部分的结果求出后半部分：
$$\begin{aligned}(x+\frac{n}{2})^{\overline{\frac{n}{2}}}&=f(x+\frac{n}{2})\\&=\sum\limits^{\frac{n}{2}}_{i=0}a_i(x+\frac{n}{2})^i\\&=\sum\limits^{\frac{n}{2}}_{i=0}\left(a_i\times\sum\limits^{i}_{j=0}\left(\binom{i}{j}\times x^j\times\left(\frac{n}{2}\right)^{i-j}\right)\right)\\&=\sum\limits^{\frac{n}{2}}_{j=0}\frac{1}{j!}\times\left(\sum\limits^{\frac{n}{2}}_{i=j}\left(i!\cdot a_i\right)\left(\frac{\left(\frac{n}{2}\right)^{i-j}}{\left(i-j\right)!}\right)\right)\times x^j\end{aligned}$$
到这里可能还不能发现它是卷积的形式，但是注意到 $i-(i-j)=j$ 所以将后面一部分的 $i$ 次项改为 $-i$ 次项即可，但是负次数多项式不好实现，可以先将后半部分对应的多项式反过来（即将负次数多项式乘上 $x^{\frac{n}{2}}$），即将$\sum\limits^{\frac{n}{2}}_{i=0}a_ii!x^i$ 和 $\sum\limits^{\frac{n}{2}}_{i=0}\frac{\left(\frac{n}{2}\right)^i}{i!}x^{\frac{n}{2}-i}$ 乘起来的多项式第 $i$ 项的系数为要求的多项式的第 $i-\frac{n}{2}$ 项的系数
于是就可以由 $x^{\overline{\frac{n}{2}}}$ 在 $O(n\log_2n)$ 时间内求出 $(x+\frac{n}{2})^{\overline{\frac{n}{2}}}$
而 $x^{\overline{\frac{n}{2}}}$ 可以递归求解
若 $n$ 是奇数，则 $x^{\overline{n}}=x^{\overline{n-1}}\times(x+n-1)$
时间复杂度 $T(n)=T(\frac{n}{2})+\Theta(nlog_2n)=\Theta(nlog_2n)$
代码
本题卡常，NTT的常数需要写得较小
const int N = 530000;
const int Mod = 167772161, g = 3;
class Polynomial {
public:
int deg, *a;
Polynomial() {
a = new int[N];
memset(a, 0, sizeof(int) * N);
}
void NTT(int, bool);
};
int r[N];
Polynomial solve(int);
int fac[N], ifac[N];
int omega[26];
int siz;
int main () {
int n;
read(n);
int lim = 1;
while (lim <= (n << 1)) {
lim <<= 1;
}
siz = __builtin_ctz(lim);
omega[25] = pow(g, (Mod - 1) >> 25, Mod);
for (int i = 24; i >= 0; --i) {
omega[i] = omega[i + 1] * 1ll * omega[i + 1] % Mod;
}
for (int i = 0; i < lim; ++i) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
}
fac[0] = 1;
for (int i = 1; i <= lim; ++i) {
fac[i] = fac[i - 1] * 1ll * i % Mod;
}
ifac[lim] = pow(fac[lim], Mod - 2, Mod);
for (int i = lim - 1; i >= 0; --i) {
ifac[i] = ifac[i + 1] * 1ll * (i + 1) % Mod;
}
Polynomial ans = solve(n);
for (int i = 0; i <= n; ++i) {
write(ans.a[i]), SP;
}
EL;
return 0;
}
void Polynomial::NTT(int lim, bool opt) {
if (opt) {
std::reverse(a + 1, a + lim);
}
int tp = siz - __builtin_ctz(lim);
for (int i = 0; i < lim; ++i) {
if (i < (r[i] >> tp)) {
std::swap(a[i], a[r[i] >> tp]);
}
}
for (int l = 1, t = 1; l < lim; l <<= 1, ++t) {
int gw = omega[t];
for (int i = 0; i < lim; i += (l << 1)) {
for (int j = 0, gn = 1; j < l; ++j, gn = gn * 1ll * gw % Mod) {
int x = a[i + j], y = a[i + j + l] * 1ll * gn % Mod;
a[i + j] = (x + y) % Mod;
a[i + j + l] = (x - y + Mod) % Mod;
}
}
}
if (opt) {
int iv = pow(lim, Mod - 2, Mod);
for (int i = 0; i < lim; ++i) {
a[i] = a[i] * 1ll * iv % Mod;
}
}
}
void calc(const Polynomial &f, int n, int lim, Polynomial &g) {
Polynomial h;
memset(g.a + n + 1, 0, sizeof(int) * (lim - n));
g.a[0] = 1, h.a[0] = f.a[0];
for (int i = 1; i <= n; ++i) {
g.a[i] = g.a[i - 1] * 1ll * n % Mod;
h.a[i] = f.a[i] * 1ll * fac[i] % Mod;
}
for (int i = 1; i <= n; ++i) {
g.a[i] = g.a[i] * 1ll * ifac[i] % Mod;
}
std::reverse(g.a, g.a + n + 1);
g.NTT(lim, false), h.NTT(lim, false);
for (int i = 0; i < lim; ++i) {
g.a[i] = g.a[i] * 1ll * h.a[i] % Mod;
}
g.NTT(lim, true);
for (int i = 0; i <= n; ++i) {
g.a[i] = g.a[i + n] * 1ll * ifac[i] % Mod;
}
memset(g.a + n + 1, 0, sizeof(int) * (lim - n));
}
Polynomial solve(int n) {
std::stack<bool> S;
while (n) {
S.push(n & 1);
n >>= 1;
}
int lim = 1;
Polynomial A, B;
n = A.a[1] = S.top();
S.pop();
while (!S.empty()) {
while (lim <= (n << 1)) {
lim <<= 1;
}
calc(A, n, lim, B);
A.NTT(lim, false), B.NTT(lim, false);
for (int i = 0; i < lim; ++i) {
A.a[i] = A.a[i] * 1ll * B.a[i] % Mod;
}
A.NTT(lim, true);
n <<= 1;
if (!S.top()) {
S.pop();
continue;
}
S.pop();
A.a[n + 1] = A.a[n];
for (int i = n; i > 0; --i) {
A.a[i] = (A.a[i - 1] + A.a[i] * 1ll * n % Mod) % Mod;
}
A.a[0] = A.a[0] * 1ll * n % Mod;
++n;
}
A.deg = n;
return A;
}
本博客采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可
本文链接：https://blog.seniorious.cc/2020/luogu-5408/