ベル数
#数え上げ #写像12相 #形式的冪級数
区別できる$ n個の玉を、区別できない$ k個以下の箱に入れる方法を$ B_{n,k}と表す。
また、箱の個数を指定しない場合は$ B_n = B_{n, n}と表す。
式
$ B_nは以下の漸化式で表される。
$ B_0 = 1
$ B_{n + 1} = \sum_{k = 0}^n \binom{n}{k} B_k
一般項は以下のように表される。なお、$ S(n, k)は第2種スターリング数。
$ B_{n, k} = \sum_{i=0}^k S(n, i) = \sum_{i=0}^k \frac{i^n}{i!} \sum_{j=0}^{k-i} \frac{(-1)^{j}}{j!}
指数型母関数により以下のように表される。
$ B(x) = \sum_{n = 0}^\infty B_n \frac{x^n}{n!} = e^{e^x-1}
実装
DPにより、すべての$ n \in [0, N], k \in [0, K] について$ B_{n,k} を$ O(NK)で求めることができる。
code:Rust
use ac_library::modint::ModIntBase;
fn bell_table<M: ModIntBase>(n_max: usize, k_max: usize) -> Vec<Vec<M>> {
let mut table = vec![vec!M::from(0); k_max + 1; n_max + 1];
table0.fill(M::from(1));
for n in 1 ..= n_max {
tablen1 = M::from(1);
if n <= k_max {
tablenn = M::from(1);
}
for k in 2 ..= (n - 1).min(k_max) {
let x = tablenk - 1 + M::from(k) * tablen - 1k - M::from(k - 1) * tablen - 1k - 1 - tablen - 1k - 2;
tablenk += x;
}
for k in n ..= k_max {
let x = tablenk - 1;
tablenk += x;
}
}
table
}
一般項の内側のシグマを累積和で処理することで、特定の$ n, kについて$ B_{n,k}を$ O(\min\{n, k\} \log n)で求めることができる。
code:Rust
use ac_library::modint::ModIntBase;
fn bell_number<M: ModIntBase>(n: usize, k: usize) -> M {
if k > n {
return bell_number(n, n);
}
let mut fact_inv = vec!M::from(1); k + 1;
fact_invk = (1 ..= k).fold(M::from(1), |f, i| f * M::from(i) ).inv();
for k in (1 .. k).rev() {
fact_invk = fact_invk + 1 * M::from(k + 1);
}
let mut inner = vec!M::from(1); k + 1;
for j in 1 ..= k {
innerj = innerj - 1 + if j % 2 == 0 { fact_invj } else { -fact_invj };
}
let mut ans = M::from(0);
for i in 0 ..= k {
ans += M::from(i).pow(n as u64) * fact_invi * innerk - i;
}
ans
}
指数型母関数を用いて、すべての$ n \in [0, N] について$ B_n を$ O(N \log N) で求めることができる。
code:Rust
fn bell_fps<M: ModIntBase>(n_max: usize) -> Vec<M> {
let mut f = FPS::<M>::from(vec!M::from(0); n_max + 1);
f1 += M::from(1);
f = f.exp();
f0 -= M::from(1);
f = f.exp();
let mut fact = vec!M::from(1); n_max + 1;
for n in 1 ..= n_max {
factn = factn - 1 * M::from(n);
}
(0 ..= n_max).map(|n| fn * factn ).collect()
}
https://ja.wikipedia.org/wiki/ベル数
https://qiita.com/drken/items/f2ea4b58b0d21621bd51
https://algo-logic.info/bell_number/
https://ngtkana.hatenablog.com/entry/2022/07/05/201141
https://lib.cp-algorithms.com/verify/poly/bell.test.cpp.html