1. 定義
第二種スタirling数 \(S(n, k)\)（または \(\genfrac{\{}{\}}{0pt}{}{n}{k}\) と表記）は、\(n\) 個のラベル付きの要素を \(k\) 個のラベルなしの空でない部分集合に分割する方法の数を表します。
組み合わせ論的意味から、\(n < k\) のとき \(\genfrac{\{}{\}}{0pt}{}{n}{k} = 0\)、\(n = k\) のとき \(\genfrac{\{}{\}}{0pt}{}{n}{k} = 1\)、\(k = 0\) のとき \(\genfrac{\{}{\}}{0pt}{}{n}{k} = [n = 0]\) となります。
2. 計算方法
2.1. 動的計画法による計算
以下の漸化式が成り立ちます。
この漸化式により、\(O(n^2)\) で前計算を行うことができ、法に対する特別な制約もありません。
この漸化式の組み合わせ論的証明を考えてみましょう。
すでに \(n-1\) 個の要素を配置し、\(n\) 番目の要素を追加するとします。追加後、合計で \(k\) 個の空でない部分集合が存在します。この新たな要素には2つの選択肢があります：新しい部分集合として単独で存在する、または既存の部分集合のいずれかに加わる。
まず、前者の場合を考えます。この場合、新たに追加された要素の状態は確定しており、それ以上の操作は不要です。残りの \(n-1\) 個の要素を \(k-1\) 個の空でない部分集合に分割する必要があります（残りの1つは新しく生成された部分集合です）。したがって、この場合の方法の数は \(\genfrac{\{}{\}}{0pt}{}{n-1}{k-1}\) となります。
次に、後者の場合を考えます。新たに追加された要素は既存の部分集合のいずれかに加わる必要があり、その選択肢は \(k\) 通りあります。また、既存の \(n-1\) 個の要素を \(k\) 個の空でない部分集合に分割する方法は \(\genfrac{\{}{\}}{0pt}{}{n-1}{k}\) 通りあります。これらを乗算すると、この場合の方法の数は \(k \cdot \genfrac{\{}{\}}{0pt}{}{n-1}{k}\) となります。
以上より、漸化式が証明されました。
2.2. 数論変換(NTT)による計算
ある行の第二種スタirling数の最初のいくつかの項を求める方法（P5395）を考えてみましょう。
2つの多項式を構築します。
NTTを用いて \(C = A \ast B\) を求めると、\(C_0 \sim C_n\) が \(\genfrac{\{}{\}}{0pt}{}{n}{0} \sim \genfrac{\{}{\}}{0pt}{}{n}{n}\) となります。
なぜこのように構築できるのでしょうか？
まず、空の部分集合を許容する場合、\(n\) 個のラベル付き要素を \(k\) 個のラベル付き部分集合に分割する方法の数を考えてみましょう。この場合、各要素は \(k\) 個の部分集合のいずれかに自由に割り当てることができるため、答えは \(k^n\) となります。
次に、空の部分集合を除外した方法の数を考慮します。\(i\) 個の集合が空であると固定し、残りの部分集合は自由に配置します。空の集合を選ぶ方法は \(\binom{k}{i}\) 通りあり、各要素は残りの \(k-i\) 個の空でない部分集合のいずれかに自由に割り当てることができるため、空でない部分集合を選ぶ方法は \((k-i)^n\) 通りあります。これらを乗算すると、総方法数は \(\binom{k}{i} \cdot (k-i)^n\) となります。
ここで \(i\) を列挙し、包含排除の原理を適用すると、次の和の式が得られます。
ここで \(f(k)\) は、\(n\) 個のラベル付き要素を \(k\) 個のラベル付き部分集合に分割する方法の数を表します。第二種スタirling数の定義はラベルなしであるため、全順列数で割る必要があります。つまり：
ここで (1) 式が第二種スタirling数の一般項の式です。
\(\sum\) の後の2つの分数はそれぞれ \(A_{k-i}\) と \(B_i\) であることに注意してください。つまり：
したがって、我々の方法は正しいことがわかります。
この証明は主に代数的な導出によって行われたため、\(A\) と \(B\) には特別な直感的な組み合わせ論的意味はありません。
3. 性質
3.1. 性質 1
以下の性質を利用して、通常のべき乗から下降階乗への変換を実現できます。
組み合わせ論的証明で示します。左辺は \(n\) 個のラベル付き要素を \(x\) 個のラベル付き集合に配置する方法の数を表します。右辺の意味を考えてみましょう。
\(\genfrac{\{}{\}}{0pt}{}{n}{k}\) は \(n\) 個のラベル付き要素を \(k\) 個のラベルなしの空でない部分集合に分割する方法の数を表します。この上で、\(x^{\underline{k}}\) を順列数 \(A_x^k\) と見なすと、右辺は次のように解釈できます：まず \(n\) 個のラベル付き要素を \(k\) 個のラベルなしの集合に配置し、次に \(k\) 個のラベルなしの集合を \(x\) 個のラベル付きグループに配置し、異なる集合が異なるグループにある方法の数。
少し回りくどいですが、左右の両方の組み合わせ論的意味が等価であることがわかります。したがって、等式は成り立ちます。
4. 例題
4.1. 例 1：P5395 第二種スタirling数・行
上記で説明した方法を適用します。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = 998244353;
const int MAXN = 1e5 + 5;
ll fac[MAXN], inv_fac[MAXN];
ll qpow(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return res;
}
void init(int n) {
fac[0] = 1;
for (int i = 1; i <= n; ++i) fac[i] = fac[i-1] * i % MOD;
inv_fac[n] = qpow(fac[n], MOD-2);
for (int i = n; i >= 1; --i) inv_fac[i-1] = inv_fac[i] * i % MOD;
}
int main() {
ios_base::sync_with_stdio(false); cin.tie(nullptr);
int N; cin >> N;
init(N);
vector polyA(N+1), polyB(N+1);
for (int i = 0; i <= N; ++i) {
polyA[i] = qpow(i, N) * inv_fac[i] % MOD;
polyB[i] = (i % 2 == 1 ? MOD - inv_fac[i] : inv_fac[i]) % MOD;
}
// NTT multiplication would be performed here
// For simplicity, we assume a function `ntt_multiply` exists
vector result(N+1);
// ntt_multiply(polyA, polyB, result, N);
// In a real implementation, you would use an NTT library
// Here we just print the expected result for demonstration
for (int i = 0; i <= N; ++i) {
cout << result[i] << " ";
}
cout << "
";
return 0;
}
4.2. 例 2：P6620 [省選联考 2020 A 卷] 組合数問題
動的計画法で前計算したスタirling数を利用します。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = 1e9 + 7;
const int MAXN = 2e3 + 5;
ll dp[MAXN][MAXN];
ll qpow(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return res;
}
int main() {
ios_base::sync_with_stdio(false); cin.tie(nullptr);
int N, X, MOD_val, M;
cin >> N >> X >> MOD_val >> M;
vector coeffs(M+1);
for (int i = 0; i <= M; ++i) cin >> coeffs[i];
// Precompute Stirling numbers of the second kind
dp[0][0] = 1;
for (int i = 1; i < MAXN; ++i) {
dp[i][1] = 1;
for (int j = 2; j <= i; ++j) {
dp[i][j] = (dp[i-1][j-1] + j * dp[i-1][j]) % MOD_val;
}
}
vector transformed(M+1);
for (int i = 0; i <= M; ++i) {
ll sum = 0;
for (int j = i; j <= M; ++j) {
sum = (sum + dp[j][i] * coeffs[j]) % MOD_val;
}
transformed[i] = sum;
}
ll answer = 0;
for (int i = 0; i <= M; ++i) {
answer = (answer + transformed[i] * qpow(X, i) % MOD_val * qpow(X + 1, N - i) % MOD_val) % MOD_val;
}
cout << answer << "
";
return 0;
}
4.3. 例 3：P2791 幼児園バスケットボール問題
事前に計算したスタirling数をクエリごとに利用します。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MOD = 1e9 + 7;
ll qpow(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = res * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return res;
}
ll comb(ll n, ll k) {
if (k > n || k < 0) return 0;
ll res = 1;
for (ll i = 1; i <= k; ++i) {
res = res * (n - i + 1) % MOD;
res = res * qpow(i, MOD-2) % MOD;
}
return res;
}
int main() {
ios_base::sync_with_stdio(false); cin.tie(nullptr);
int total_elements, total_balls, test_cases, max_stirling;
cin >> total_elements >> total_balls >> test_cases >> max_stirling;
// Precompute Stirling numbers of the second kind up to max_stirling
vector> stirling(max_stirling + 1, vector(max_stirling + 1, 0));
stirling[0][0] = 1;
for (int i = 1; i <= max_stirling; ++i) {
stirling[i][1] = 1;
for (int j = 2; j <= i; ++j) {
stirling[i][j] = (stirling[i-1][j-1] + j * stirling[i-1][j]) % MOD;
}
}
while (test_cases--) {
int current_n, current_m, current_k;
cin >> current_n >> current_m >> current_k;
ll inv_comb = qpow(comb(current_n, current_k), MOD-2);
ll current_answer = 0;
for (int i = 0; i <= min({current_k, current_m, max_stirling}); ++i) {
ll term = stirling[max_stirling][i];
term = term * comb(current_m, i) % MOD;
term = term * comb(current_n - i, current_k - i) % MOD;
term = term * comb(max_stirling, i) % MOD; // jc[i] is C(max_stirling, i)
current_answer = (current_answer + term) % MOD;
}
current_answer = current_answer * inv_comb % MOD;
cout << current_answer << "
";
}
return 0;
}