贝尔数（Bell Number）
贝尔数表示的是基数为\(n\)的集合的划分成两两不想交的子集的划分方案数，记作\(B_n\)
贝尔数的前几项为：\(1,1,2,5,15,52,203..........\)
1、组合意义推导
我们考虑第\(n+1\)个元素\(b_{n+1}\)，那么可以对其分类讨论：
（1）、如果\(b_{n+1}\)单独分为一个子集，则不用从\(n\)个元素中选择元素出来和其组成一个子集，故此时的方案数为：\(C_n^0B_n\)
（2）、如果\(b_{n+1}\)和一个元素组成一个子集，则需要从\(n\)个元素中选出一个元素与之组成一个子集，故此时的方案数为：\(C_n^1B_{n-1}\)
（3）、如果\(b_{n+1}\)和两个元素组成一个子集，则需要从\(n\)个元素中选出两个元素与之组成一个子集，故此时的方案数为：\(C_n^2B_{n-2}\)
\(..........\)
以此类推，那么我们可以得到递推式： \[ \begin{aligned} B_{n+1}&=\sum_{k=0}^{n}C_n^{k}B_{n-k}\\ &=\sum_{k=0}^{n}C_n^{k}B_k \end{aligned} \] 另一种意义：
我们还可以考虑另外一种组合意义：
（1）、将\(n\)个元素分成\(1\)个子集，此时方案数：\(\begin{Bmatrix}n\\\\1\end{Bmatrix}\)
（2）、将\(n\)个元素分成\(2\)个子集，此时方案数：\(\begin{Bmatrix}n\\\\2\end{Bmatrix}\)
（3）、将\(n\)个元素分成\(3\)个子集，此时方案数：\(\begin{Bmatrix}n\\\\3\end{Bmatrix}\)
\(........\)
以此类推，我们可以得到新的递推式： \[ B_n=\sum_{k=0}^n\begin{Bmatrix}n\\\\k\end{Bmatrix} \] 这种推导方式也推导出了，贝尔数与第二类斯特林数的关系
2、生成函数推导
\[ 设b(x)为序列【B_1,B_2,B_3....B_n】的生成函数 \\ 则b(x)=\sum_{i=0}^{\infty}B_i\frac{x^i}{i!}\\ \therefore b(x)=\sum_{i=0}^{\infty}\sum_{k=0}^{i-1}C_{i-1}^{k}B_k\frac{x^i}{i!}\\改变i与k的求和顺序\\ \begin{aligned} \therefore b(x)&=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}C_{i-1}^{k}\frac{x^i}{i!}\\ &=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{(i-1)!}{(i-k-1)!k!}\frac{x^i}{i!}\\ &=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{1}{i(i-k-1)!k!}x^i\\ \end{aligned} \\ 对b(x)求关于x的导数\\ \begin{aligned} \therefore b^{`}(x)&=\sum_{k=0}^{\infty}B_k\sum_{i=k+1}^{\infty}\frac{1}{(i-k-1)!k!}x^{i-1}\\ &=\sum_{k=0}^{\infty}\frac{B_k}{k!}\sum_{i=k+1}^{\infty}\frac{1}{(i-k-1)!}x^{i-1}\\ &=\sum_{k=0}^{\infty}\frac{B_k}{k!}\sum_{i=0}^{\infty}\frac{1}{i!}x^{i-1}x^{k+1}\\ &=\sum_{k=0}^{\infty}\frac{B_kx^k}{k!}\sum_{i=0}^{\infty}\frac{1}{i!}x^i\\ &=b(x)e^x\\ \end{aligned}\\ \therefore \frac{b^{`}(x)}{b(x)}=e^x \\ 又\because ln(f(x))=\frac{f^{`}(x)}{f(x)}\\ \therefore ln(b(x))=e^x \\ \therefore b(x)=e^{e^x+c},c为常数 \\代入计算得C=-1\\ \therefore b(x)=e^{e^x-1} \]
3、性质
引理：
\[ B_{n+m}=\sum_{i=0}^{n}C_n^iB_i\sum_{j=0}^mj^{n-i}\begin{Bmatrix}m\\\\j\end{Bmatrix} \]
组合意义证明：
把\(n+m\)个元素分成\(n\)个元素和\(m\)个元素；首先枚举\(m\)个元素划分成\(j\)个集合，此时有\(\begin{Bmatrix}m\\\\j\end{Bmatrix}\);再枚举\(n\)个元素中选\(i\)个出来任意划分 成子集，此时有\(C_n^iB_i\)；还剩下\(n-i\)个元素，那么这\(n-i\)个元素必然在\(j\)个集合中，此时有\(j^{n-i}\)，然后根据乘法原理将其相乘，便是最后的结果
3.1、Dobinski公式：
\[ B_n=\frac{1}{e}\sum_{k=0}^{\infty}\frac{k^n}{k!} \]
期望为\(n\)的泊松分数的\(n\)次矩
3.2、Touchard同余
\[ 若对\forall p\in primes,有B_{p+n}\equiv B_n+B_{n+1}(mod\space p) \]
证明： \[ \begin{aligned} B_{n+p}&=\sum_{i=0}^{n}C_n^iB_i\sum_{j=0}^pj^{n-i}\begin{Bmatrix}p\\\\j\end{Bmatrix}\\ &=\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\sum_{i=0}^{n}C_n^iB_ij^{n-i}\\ 又\because 对于\forall &j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p) \\ \therefore B_{n+p}\space mod\space p&=\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p\\ &=(\sum_{j=0}^p\begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p)(\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p)\space mod\space p\\ &=(\begin{Bmatrix}p\\\\1\end{Bmatrix}+0+\begin{Bmatrix}p\\\\p\end{Bmatrix})(\sum_{i=0}^{n}C_n^iB_ij^{n-i}\space mod\space p)\space mod\space p\\ &=(\begin{Bmatrix}p\\\\1\end{Bmatrix}(\sum_{i=0}^{n}C_n^iB_i1^{n-i}\space mod\space p)\space mod\space p+\begin{Bmatrix}p\\\\p\end{Bmatrix}(\sum_{i=0}^{n}C_n^iB_ip^{n-i}\space mod\space p)\space mod\space p)\space mod\space p\\ &=(\sum_{i=0}^{n}C_n^iB_i1^{n-i}+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ &=(\sum_{i=0}^{n}C_n^iB_i+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ &=(B_{n+1}+\sum_{i=0}^{n}C_n^iB_ip^{n-i})\space mod\space p\\ 又\because 对\forall i>0,&p^i\equiv 0(mod \space p )\\ \therefore \sum_{i=0}^{n}C_n^i&B_ip^{n-i}\space mod\space p=C_n^nB_np^{n-n}\space mod\space p=B_n\space mod\space p\\ \therefore B_{n+p}&\equiv B_{n+1}+B_n(\space mod\space p)\\故：结论得证 \end{aligned} \] TIPS：
证明：\(对于\forall j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p)\) \[ 由第二类斯特林数的通项公式得：\begin{Bmatrix}p\\\\j\end{Bmatrix}=\sum_{k=0}^m\frac{(-1)^k}{k!}\frac{(j-k)^p}{(j-k)!}\\ \begin{aligned} \therefore \begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p&=(\sum_{k=0}^j\frac{(-1)^k}{k!}\frac{(j-k)^p}{(j-k)!})\space mod\space p\\ &=(\sum_{k=0}^j(\frac{(-1)^k}{k!}\space mod\space p)(\frac{(j-k)^p}{(j-k)!}\space mod\space p)\space mod\space p)\space mod\space p\\ \end{aligned} \\ 由费马小定理可得:对于\forall p\in primes且\forall a 满足gcd(a,p)=1，都有a^{p-1}\equiv 1(\space mod\space p)\\ \therefore \frac{(j-k)^p}{(j-k)!}\space mod\space p=\frac{(j-k)}{(j-k)!}\space mod\space p=\frac{1}{(j-k-1)!}\space mod\space p\\ \therefore \begin{Bmatrix}p\\\\j\end{Bmatrix}\space mod\space p=(\sum_{k=0}^j\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!})\space mod\space p\\ 这个时候不妨考虑生成函数A(x)=\sum_{i=0}^{\infty}(\sum_{k=0}^i\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!})x^i\\ 通过观察发现A(x)其实是卷积的形式\\ 不妨设B(x)=\sum_{i=0}^{\infty}\frac{(-1)^i}{i!}x^i,C(x)=\sum_{i=0}^{\infty}\frac{1}{(i-1)!})x^i \\ 则C(x)=\sum_{i=0}^{\infty}\frac{1}{(i-1)!})x^i=x\sum_{i=0}^{\infty}\frac{1}{(i-1)!}x^{i-1}=x\sum_{i=0}^{\infty}\frac{1}{i!}x^i\\ 由e^x的泰勒展开可得e^x=\sum_{i=0}^{\infty}\frac{1}{i!}x^i,e^{-x}=\frac{(-1)^i}{i!}x^i \\ 那么可得A(x)=B(x)C(x)=e^{-x}xe^x=x \\ \therefore A(x)=x\\ \therefore 对于i>1,\sum_{k=0}^i\frac{(-1)^k}{k!}\frac{1}{(j-k-1)!}=0 \\ \therefore对于\forall j\in(1,p)且p\in primes,有 \begin{Bmatrix}p\\\\j\end{Bmatrix}\equiv0(mod\space p) \\结论得证 \]
4、贝尔三角形
构造三角矩阵
满足一下递推式： \[ a_{n,m}= \begin{cases} 1,n=1,m=1\\ a_{n-1,n-1},n>1,m=1\\ a_{n,m-1}+a_{n-1,m-1},n>1,m>1 \end{cases} \]
\[ \begin{aligned} & 1 \\ & 1\quad\qquad 2 \\ & 2\quad\qquad 3\quad\qquad 5 \\ & 5\quad\qquad 7\quad\qquad 10\qquad 15 \\ & 15\qquad 20\qquad 27\qquad 37\qquad 52 \\ & 52\qquad 67\qquad 87\qquad 114\qquad 151\qquad 203\\ & 203\qquad 255\qquad 322\qquad 409\qquad 523\qquad 674\qquad 877 \\ \end{aligned} \] 每行的首项为贝尔数
5、快速求贝尔数在模意义下的第\(n\)项
题目链接：https://vjudge.csgrandeur.cn/problem/HDU-4767
题意：求\(B_n \space mod \space 95041567\)
思路：
首先，我们发现\(95041567\)并不是一个素数；但是我们发现\(95041567\)是多个不相同的素数的乘积，因此我们可以考虑CRT合并，所以我们可以列出以下方程： \[ 95041567=\Pi_ip_i\\ \begin{cases} x\equiv a_1(\space mod \space p_1)\\ x\equiv a_2(\space mod \space p_2)\\ ....\\ x\equiv a_n(\space mod \space p_k)\\ \end{cases} \] 因此我们将问题转化为了求出一组\(a_i=B_n \space mod \space p_i\\\)
而对于每一组素数\(p_i\)，我们可以利用\(Touchard同余\) \[ 将 B_{n+p}\equiv B_{n+1}+B_n(\space mod\space p)改写成模意义下的矩阵相乘的形式: \\ \begin{bmatrix} B_{n+p-1}\\ B_{n+p}\\ B_{n+p+1}\\ \vdots \\ B_{n+p+p-1} \end{bmatrix}= \begin{bmatrix} &0 &0 &0 &\dots&0&1\\ &1 &1 &0 &\dots&0&0\\ &0 &1 &1 &\dots&0&0\\ &\vdots &\vdots &\vdots &\ddots &\vdots&\vdots\\ &0 &0 &0 &\dots&1&1\ \end{bmatrix}\times \begin{bmatrix} B_{n}\\ B_{n+1}\\ B_{n+2}\\ \vdots \\ B_{n+p-1} \end{bmatrix} \] 那么可以利用矩阵的递推推出第\(n\)项与前几项的关系： \[ \begin{bmatrix} B_{n}\\ B_{n+1}\\ B_{n+2}\\ \vdots \\ B_{n+p-1} \end{bmatrix}= \begin{bmatrix} &0 &0 &0 &\dots&0&1\\ &1 &1 &0 &\dots&0&0\\ &0 &1 &1 &\dots&0&0\\ &\vdots &\vdots &\vdots &\ddots &\vdots&\vdots\\ &0 &0 &0 &\dots&1&1\ \end{bmatrix}^{\frac{n}{p-1}}\times \begin{bmatrix} B_{0}\\ B_{1}\\ B_{2}\\ \vdots \\ B_{p-1} \end{bmatrix} \] 我们发现中间那个矩阵可以利用矩阵快速幂求出，那么我们只需要预处理出前\(p_i\)个\(B_i\)，然后利用矩阵乘法就可以得到\(B_n\space mod \space p_i\)
现在我们已经知道一组\(a_i\),那么利用CRT合并得到的结果，便是最后的答案！
#include <iostream>
#include <algorithm>
#include <cmath>
#include <string.h>
#include <vector>
using namespace std;
#define int long long
#define endl "\n"
const int maxn = 1e2 + 10;
const int p = 95041567;
int bell[maxn][maxn], ans[maxn];
//矩阵快速幂
int cnt;
struct jz
{
int m[maxn][maxn];
jz()
{
memset(m, 0, sizeof(m));
}
} ot;
jz mul(jz a, jz b, int mod)
{
jz c;
for (int i = 1; i <= cnt; i++)
{
for (int j = 1; j <= cnt; j++)
{
c.m[i][j] = 0;
for (int k = 1; k <= cnt; k++)
c.m[i][j] = (c.m[i][j] + a.m[i][k] * b.m[k][j] % mod) % mod;
}
}
return c;
}
jz ksm(jz a, int b, int mod)
{
jz ans;
for (int i = 1; i <= cnt; i++)
for (int j = 1; j <= cnt; j++)
ans.m[i][j] = 1 * (i == j);
while (b)
{
if (b % 2)
ans = mul(ans, a, mod);
a = mul(a, a, mod);
b >>= 1;
}
return ans;
}
//模数分解
vector<int> mod;
vector<int> a;
void div(int n)
{
for (int i = 2; i * i <= n; i++)
{
if (n % i == 0)
{
mod.push_back(i);
while (n % i == 0)
{
n /= i;
}
}
}
if (n > 1)
{
mod.push_back(n);
}
}
//CRT合并
void exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return;
}
exgcd(b, a % b, x, y);
int temp = x;
x = y;
y = temp - a / b * y;
}
int CRT(vector<int> a, vector<int> mod)
{
int n = mod.size();
int M = 1, ans = 0;
for (int i = 0; i < n; i++)
{
M *= mod[i];
}
for (int i = 0; i < n; i++)
{
int m = M / mod[i];
int b, y;
exgcd(m, mod[i], b, y);
ans = (ans + a[i] * m * b % M) % M;
}
return (ans % M + M) % M;
}
//贝尔三角形递推前几项贝尔数
void init()
{
bell[1][1] = 1;
for (int i = 2; i < maxn; i++)
{
for (int j = 1; j < maxn; j++)
{
if (j == 1)
{
bell[i][j] = bell[i - 1][i - 1];
}
else
{
bell[i][j] = (bell[i][j - 1] + bell[i - 1][j - 1]) % p;
}
}
}
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
div(95041567);
init();
int t;
cin >> t;
while (t--)
{
int n;
cin >> n;
a.clear();
for (auto mo : mod)
{
cnt = mo;
jz temp;
temp.m[1][mo] = 1;
for (int i = 2; i <= mo; i++)
{
for (int j = 1; j <= mo; j++)
{
if (i == j || i == j + 1)
{
temp.m[i][j] = 1;
}
}
}
ot = ksm(temp, n / (mo - 1), mo);
memset(ans, 0, sizeof(ans));
for (int i = 1; i <= mo; i++)
{
for (int j = 1; j <= mo; j++)
{
ans[i] = (ans[i] + ot.m[i][j] * bell[j][1] % mo) % mo;
}
}
a.push_back(ans[n % (mo - 1) + 1]);
}
int ans = CRT(a, mod);
cout << ans << endl;
}
return 0;
}
