斯特林数
定义
斯特林数分为第一类斯特林数和第二类斯特林数两种。
-
第一类斯特林数：将 \(n\) 个元素划分成 \(m\) 个圆排列的方案数。
-
第二类斯特林数：将 \(n\) 个元素划分成 \(m\) 个无序集合的方案数。
这两类组合意义有时候会经常用到。
第一类斯特林数考虑的是划分有序集合，但是等价于同一圆排列的
第二类斯特林数类似于逆无标号计数（……？）考虑将全集划分成若干子集时可以使用。
递推
直接根据组合意义。
第一类斯特林数：考虑第 \(n\) 个元素所处的圆排列。
-
自成一个圆排列
-
插入已有的圆排列。
所以 \(S_1(n, k) = S_1(n - 1, k - 1) + (n - 1) S_1(n - 1, k)\).
第二类斯特林数：类似。
有 \(S_2(n, k) = S_2(n - 1, k - 1) + k S_2(n - 1, k)\).
快速求值
朴素递推两类斯特林数的复杂度是 \(O(nk)\)，如果 \(n\) 或 \(m\) 是固定的，可以用 poly 优化复杂度到 \(O(n \log n)\).
虽然真要出斯特林数都可以暴力递推就是了
茴字的四种写法
-
求一整行第一类斯特林数
-
求一整列第一类斯特林数
-
求一整行第二类斯特林数
-
求一整列第二类斯特林数
前置知识：上升幂、下降幂。
下降幂
考虑求出 \(x^{\overline{n}}\).
注意到 \(x^{\overline{2n}} = x^{\overline{n}} (x + n)^{\overline{n}}\).
于是套上多项式移位，类似快速幂地倍增求就行。
void times(ll *F, ll *G, int len, int lim)
{
int k = 1;
while (k < (len << 1)) k <<= 1;
NTT(F, k), NTT(G, k);
for (int i = 0; i < k; i++) F[i] = F[i] * G[i] % mod;
INTT(F, k), INTT(G, k);
}
void fplus(ll *F, ll *G, int len, int c)
{
for (int i = 0; i < len; i++) Ft[len - i - 1] = G[i] * fac[i] % mod;
ll coe = 1;
for (int i = 0; i < len; i++, coe = coe * c % mod) F[i] = coe * invf[i] % mod;
times(Ft, F, len, len);
for (int i = 0; i < len; i++) F[len - i - 1] = Ft[i] * invf[len - i - 1] % mod;
for (int i = len; i < (len << 1); i++) Ft[i] = 0;
}
void solve(ll *F, ll *G, int n)
{
if (n == 1) return F[0] = 0, F[1] = 1, void();
if (n & 1)
{
solve(F, G, n - 1);
F[n] = 0;
for (int i = n; i > 0; i--) F[i] = (F[i - 1] + (n - 1) * F[i]) % mod;
F[0] = F[0] * (n - 1) % mod;
return;
}
solve(F, G, n >> 1);
fplus(G, F, (n >> 1) + 1, n >> 1);
times(F, G, (n >> 1) + 1, n + 1);
}
第二类斯特林数·行
第二类斯特林数有性质：\(m^n = \sum\limits_{i = 0}^m {m \choose i} {n \brace i} i!\)，组合意义是 选出非空集合 - 划分到钦定集合 - 钦定顺序。
根据二项式反演得 \({n \brace m} m! = \sum\limits_{i = 0}^m (-1)^{m - i} {m \choose i} i^n\).
拆开组合数整理得到 \({n \brace m} = \sum\limits_{i = 0}^m \frac{(-1)^{m - i}}{(m - i)!} \frac{i^n}{i!}\).
令 \(F(x) = \sum\limits_{i = 0}^m \frac{(-1)^i}{i!}, G(x) = \sum\limits_{i = 0}^m \frac{i^n}{i!}\)，上式等价于 \(F * G\).
时间复杂度 \(O(n \log n)\).
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 1e6 + 5;
const int mod = 167772161;
const int g = 3;
int n;
int rev[maxn], invf[maxn];
ll F[maxn], G[maxn], wp[maxn];
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod;
power >>= 1;
}
return res;
}
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len, mod);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}
int main()
{
scanf("%d", &n);
invf[0] = invf[1] = 1;
for (int i = 2; i <= n; i++) invf[i] = 1ll * (mod - mod / i) * invf[mod % i] % mod;
for (int i = 1; i <= n; i++) invf[i] = 1ll * invf[i] * invf[i - 1] % mod;
for (int i = 0; i <= n; i++) F[i] = qpow(i, n, mod) * invf[i] % mod, G[i] = (i & 1 ? mod - invf[i] : invf[i]);
int k = 1;
while (k <= (n << 1)) k <<= 1;
NTT(F, k), NTT(G, k);
for (int i = 0; i < k; i++) F[i] = F[i] * G[i] % mod;
INTT(F, k);
for (int i = 0; i <= n; i++) printf("%lld ", F[i]);
putchar('\n');
return 0;
}
第一类斯特林数·行
第一类斯特林数的性质：\(x^{\overline{n}} = \sum\limits_{i = 0}^n {n \brack i} x^i\)，问题转化成求 \(x^{\overline{n}}\).
这个可以倍增求：\(x^{\overline{2n}} = x^{\overline{n}} (x + n)^{\overline{n}}\).
也就是说现在有一个多项式 \(F(x)\)，我们希望求出 \(F(x + c)\)，也就是多项式移位。
表示出 \(F(x + c)\)：
\(F(x + c) = \sum\limits_{i = 0}^{n - 1} f_i (x + c)^i\).
考虑二项式定理展开：
\(F(x + c) = \sum\limits_{i = 0}^{n - 1} f_i \sum\limits_{j = 0}^{i} {i \choose j} x^i c^{i - j}\).
组合数拆开，交换一下求和顺序：
\(F(x + c) = \sum\limits_{j = 0}^{n - 1} x^j \frac{1}{j!} \sum\limits_{i = j}^{n - 1} i! f_i \frac{c^{i - j}}{(i - j)!}\).
发现是差卷积的形式，于是反转再卷积一下就可以倍增求了。
时间复杂度 \(T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)\).
#include <cstdio>
#include <algorithm>
using namespace std;
#define swap(x, y) x ^= y ^= x ^= y
typedef long long ll;
const int maxn = (1 << 19);
const int mod = 167772161;
const int g = 3;
int n;
int rev[maxn], fac[maxn], invf[maxn];
ll wp[maxn], F[maxn], G[maxn], Ft[maxn];
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod, power >>= 1;
}
return res;
}
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len, mod);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}
void times(ll *F, ll *G, int len, int lim)
{
int k = 1;
while (k < (len << 1)) k <<= 1;
NTT(F, k), NTT(G, k);
for (int i = 0; i < k; i++) F[i] = F[i] * G[i] % mod;
INTT(F, k), INTT(G, k);
}
void fplus(ll *F, ll *G, int len, int c)
{
for (int i = 0; i < len; i++) Ft[len - i - 1] = G[i] * fac[i] % mod;
ll coe = 1;
for (int i = 0; i < len; i++, coe = coe * c % mod) F[i] = coe * invf[i] % mod;
times(Ft, F, len, len);
for (int i = 0; i < len; i++) F[len - i - 1] = Ft[i] * invf[len - i - 1] % mod;
for (int i = len; i < (len << 1); i++) Ft[i] = 0;
}
void solve(ll *F, ll *G, int n)
{
if (n == 1) return F[0] = 0, F[1] = 1, void();
if (n & 1)
{
solve(F, G, n - 1);
F[n] = 0;
for (int i = n; i > 0; i--) F[i] = (F[i - 1] + (n - 1) * F[i]) % mod;
F[0] = F[0] * (n - 1) % mod;
return;
}
solve(F, G, n >> 1);
fplus(G, F, (n >> 1) + 1, n >> 1);
times(F, G, (n >> 1) + 1, n + 1);
}
int main()
{
scanf("%d", &n);
fac[0] = invf[0] = fac[1] = invf[1] = 1;
for (int i = 2; i <= n; i++) fac[i] = 1ll * fac[i - 1] * i % mod, invf[i] = 1ll * (mod - mod / i) * invf[mod % i] % mod;
for (int i = 2; i <= n; i++) invf[i] = 1ll * invf[i] * invf[i - 1] % mod;
solve(F, G, n);
for (int i = 0; i <= n; i++) printf("%lld ", F[i]);
putchar('\n');
return 0;
}
第二类斯特林数·列
观察递推公式：
\({n \brace m} = {n - 1 \brace m - 1} + m {n - 1 \brace m}\).
考虑生成函数 \(H_k(x) = \sum\limits_{i = 0} {i \brace k} x^i\).
\(H_k(x) = \sum\limits_{i = 0} ({i - 1 \brace k - 1} + k {i - 1 \brace k}) x^i\).
简单整理再移位可得：
\(H_k(x) = x \sum\limits_{i = 0} {i \brace k - 1} x^i + kx \sum\limits_{i = 0} {i \brace k} x^i = xH_{k - 1}(x) + kx H_k(x)\).
于是可以得到 \(H_k(x)\) 的封闭形式：
\(H_k(x) = \frac{x}{1 - kx} H_{k - 1}(x), H_0(x) = 1\).
现在的问题是：计算 \(H_k(x) = \prod\limits_{i = 1}^k (1 - ix) = \prod\limits_{i = 1}^k x(\frac{1}{x} - i) = x^k \prod\limits_{i = 1}^k (\frac{1}{x} - i)\).
这里有结论：\(\prod\limits_{i = 1}^k x(\frac{1}{x} - i)\) 反转过来就是 \(x^k \prod\limits_{i = 1}^k (\frac{1}{x} - i)\)，证明之后补。
注意到 \(\prod\limits_{i = 1}^k x(\frac{1}{x} - i) = \frac{x^{\overline{k + 1}}}{x}\)，考虑求 \(x\) 的 \(k + 1\) 次上升幂。
直接倍增 \(O(n \log n)\) 做就行。
#include <cstdio>
#include <algorithm>
using namespace std;
#define swap(x, y) x ^= y ^= x ^= y
typedef long long ll;
const int maxn = (1 << 19);
const int mod = 167772161;
const int g = 3;
int n;
int rev[maxn], fac[maxn], invf[maxn];
ll wp[maxn], F[maxn], G[maxn], Ft[maxn];
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod, power >>= 1;
}
return res;
}
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len, mod);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i++) A[i] = 1ll * A[i] * inv % mod;
}
void times(ll *F, ll *G, int len, int lim)
{
int k = 1;
while (k < (len << 1)) k <<= 1;
NTT(F, k), NTT(G, k);
for (int i = 0; i < k; i++) F[i] = F[i] * G[i] % mod;
INTT(F, k), INTT(G, k);
}
void fplus(ll *F, ll *G, int len, int c)
{
for (int i = 0; i < len; i++) Ft[len - i - 1] = G[i] * fac[i] % mod;
ll coe = 1;
for (int i = 0; i < len; i++, coe = coe * c % mod) F[i] = coe * invf[i] % mod;
times(Ft, F, len, len);
for (int i = 0; i < len; i++) F[len - i - 1] = Ft[i] * invf[len - i - 1] % mod;
for (int i = len; i < (len << 1); i++) Ft[i] = 0;
}
void solve(ll *F, ll *G, int n)
{
if (n == 1) return F[0] = 0, F[1] = 1, void();
if (n & 1)
{
solve(F, G, n - 1);
F[n] = 0;
for (int i = n; i > 0; i--) F[i] = (F[i - 1] + (n - 1) * F[i]) % mod;
F[0] = F[0] * (n - 1) % mod;
return;
}
solve(F, G, n >> 1);
fplus(G, F, (n >> 1) + 1, n >> 1);
times(F, G, (n >> 1) + 1, n + 1);
}
int main()
{
scanf("%d", &n);
fac[0] = invf[0] = fac[1] = invf[1] = 1;
for (int i = 2; i <= n; i++) fac[i] = 1ll * fac[i - 1] * i % mod, invf[i] = 1ll * (mod - mod / i) * invf[mod % i] % mod;
for (int i = 2; i <= n; i++) invf[i] = 1ll * invf[i] * invf[i - 1] % mod;
solve(F, G, n);
for (int i = 0; i <= n; i++) printf("%lld ", F[i]);
putchar('\n');
return 0;
}
第一类斯特林数·列
第一类斯特林数的 EGF：\(\frac{(-\ln(1 - x)^k)}{k!}\)
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = (1 << 19);
const int mod = 167772161;
const int g = 3;
int n, k;
int rev[maxn], fac[maxn], inv[maxn];
ll F[maxn], G[maxn];
ll wp[maxn], Ft[maxn], Rt[maxn], ft[maxn], rt[maxn], Fn[maxn], Rn[maxn];
ll qpow(ll base, ll power, ll mod)
{
ll res = 1;
while (power)
{
if (power & 1) res = res * base % mod;
base = base * base % mod;
power >>= 1;
}
return res;
}
void calc_rev(int k) { for (int i = 1; i < k; i++) rev[i] = (rev[i >> 1] >> 1 | (i & 1 ? k >> 1 : 0)); }
void NTT(ll *A, int n)
{
calc_rev(n);
for (int i = 1; i < n; i++)
if (rev[i] > i) swap(A[i], A[rev[i]]);
for (int len = 2, m = 1; len <= n; m = len, len <<= 1)
{
ll wn = qpow(g, (mod - 1) / len, mod);
wp[0] = 1;
for (int i = 1; i <= len; i++) wp[i] = wp[i - 1] * wn % mod;
for (int l = 0, r = len - 1; r <= n; l += len, r += len)
{
int w = 0;
for (int p = l; p < l + m; p++, w++)
{
ll x = A[p], y = wp[w] * A[p + m] % mod;
A[p] = (x + y) % mod, A[p + m] = (x - y + mod) % mod;
}
}
}
}
void INTT(ll *A, int n)
{
NTT(A, n);
reverse(A + 1, A + n);
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i++) A[i] = A[i] * inv % mod;
}
void diffp(ll *f, ll *der, int n)
{
for (int i = 1; i < n; i++) der[i - 1] = f[i] * i % mod;
der[n - 1] = 0;
}
void intep(ll *f, ll *inte, int n)
{
for (int i = 1; i < n; i++)
if (!inv[i]) inv[i] = (i == 1 ? 1 : 1ll * inv[mod % i] * (mod - mod / i) % mod);
for (int i = 1; i < n; i++) inte[i] = f[i - 1] * inv[i] % mod;
inte[0] = 0;
}
void invp(ll *f, ll *r, int n)
{
int k = 1;
while (k < n) k <<= 1;
r[0] = qpow(f[0], mod - 2, mod);
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Rt[i] = r[i], Ft[i] = f[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = 0; i < m; i++) Rt[i] = 0; Rt[0] = 1;
for (int i = 0; i < len; i++) Ft[i] = r[i];
NTT(Ft, len), NTT(Rt, len);
for (int i = 0; i < len; i++) Rt[i] = Rt[i] * Ft[i] % mod;
INTT(Rt, len);
for (int i = m; i < len; i++) r[i] = (r[i] * 2ll - Rt[i] + mod) % mod;
}
memset(Ft, 0, k * sizeof(ll));
memset(Rt, 0, k * sizeof(ll));
for (int i = n; i < k; i++) r[i] = 0;
}
void lnp(ll *f, ll *ln, int n)
{
diffp(f, ft, n), invp(f, rt, n);
int k = 1;
while (k < (n << 1)) k <<= 1;
NTT(ft, k), NTT(rt, k);
for (int i = 0; i < k; i++) ft[i] = ft[i] * rt[i] % mod;
INTT(ft, k);
intep(ft, ln, n);
for (int i = 0; i < k; i++) ft[i] = rt[i] = 0;
}
void expp(ll *f, ll *exp, int n)
{
int k = 1;
while (k < n) k <<= 1;
exp[0] = 1;
for (int len = 2, m = 1; len <= k; m = len, len <<= 1)
{
for (int i = 0; i < len; i++) Fn[i] = exp[i];
lnp(Fn, Rn, len);
for (int i = 0; i < len; i++) Rn[i] = (f[i] - Rn[i] + mod) % mod;
Rn[0] = (Rn[0] + 1) % mod;
NTT(Fn, len << 1), NTT(Rn, len << 1);
for (int i = 0; i < (len << 1); i++) Fn[i] = Fn[i] * Rn[i] % mod;
INTT(Fn, len << 1);
for (int i = 0; i < len; i++) exp[i] = Fn[i];
}
for (int i = 0; i < (k << 1); i++) Fn[i] = Rn[i] = 0;
for (int i = n; i < k; i++) exp[i] = 0;
}
int main()
{
scanf("%d%d", &n, &k);
if (k > n)
{
for (int i = 0; i <= n; i++) printf("%d ", 0);
return 0;
}
n++;
fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = 1ll * fac[i - 1] * i % mod;
for (int i = 0; i < n; i++) F[i] = qpow(i + 1, mod - 2, mod);
lnp(F, G, n);
for (int i = 0; i < n; i++) G[i] = G[i] * k % mod;
memset(F, 0, sizeof(F));
expp(G, F, n);
for (int i = n - k - 1; i >= 0; i--) F[i + k] = F[i];
for (int i = 0; i < k; i++) F[i] = 0;
int invk = 1;
for (int i = 1; i <= k; i++) invk = 1ll * invk * i % mod;
invk = qpow(invk, mod - 2, mod);
for (int i = 0; i < n; i++) F[i] = F[i] * fac[i] % mod * invk % mod;
for (int i = 0; i < n; i++) printf("%lld ", F[i]);
return 0;
}
