【ACM数论】贝尔数

贝尔数（Bell Number）

贝尔数表示的是基数为$n$的集合的划分成两两不想交的子集的划分方案数，记作$B_n$

贝尔数的前几项为：$1,1,2,5,15,52,\mathrm{203..........}$

1、组合意义推导

我们考虑第$n+1$个元素$b_{n+1}$，那么可以对其分类讨论：

（1）、如果$b_{n+1}$单独分为一个子集，则不用从$n$个元素中选择元素出来和其组成一个子集，故此时的方案数为：$C_n^0{B}_n$

（2）、如果$b_{n+1}$和一个元素组成一个子集，则需要从$n$个元素中选出一个元素与之组成一个子集，故此时的方案数为：$C_n^1{B}_{n-1}$

（3）、如果$b_{n+1}$和两个元素组成一个子集，则需要从$n$个元素中选出两个元素与之组成一个子集，故此时的方案数为：$C_n^2{B}_{n-2}$

$..........$

以此类推，那么我们可以得到递推式：
$$\begin{array}{rl}{B}_{n+1}& =\sum _{k=0}^n{C}_n^k{B}_{n-k}\\ & =\sum _{k=0}^n{C}_n^k{B}_k\end{array}$$
另一种意义：

我们还可以考虑另外一种组合意义：

（1）、将$n$个元素分成$1$个子集，此时方案数：$\left\{\begin{array}{c}n\\ \\ 1\end{array}\right\}$

（2）、将$n$个元素分成$2$个子集，此时方案数：$\left\{\begin{array}{c}n\\ \\ 2\end{array}\right\}$

（3）、将$n$个元素分成$3$个子集，此时方案数：$\left\{\begin{array}{c}n\\ \\ 3\end{array}\right\}$

$........$

以此类推，我们可以得到新的递推式：
$$B_n=\sum _{k=0}^n\left\{\begin{array}{c}n\\ \\ k\end{array}\right\}$$
这种推导方式也推导出了，贝尔数与第二类斯特林数的关系

2、生成函数推导

$$\begin{gathered} 设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{array}{rl}\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{array} \\ 对b(x)求关于x的导数 \\ \begin{array}{rl}\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_k{x}^k}{k!}\sum _{i=0}^{\infty }\frac{1}{i!}{x}^i\\ & =b(x)e^x\end{array} \\ \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} \end{gathered}$$

3、性质

引理：

$$B_{n+m}=\sum _{i=0}^n{C}_n^i{B}_i\sum _{j=0}^m{j}^{n-i}\left\{\begin{array}{c}m\\ \\ j\end{array}\right\}$$

组合意义证明：

​ 把$n+m$个元素分成$n$个元素和$m$个元素；首先枚举$m$个元素划分成$j$个集合，此时有$\left\{\begin{array}{c}m\\ \\ j\end{array}\right\}$;再枚举$n$个元素中选$i$个出来任意划分 成子集，此时有$C_n^i{B}_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}(modp)$$

​ 证明：
$$\begin{array}{rl}{B}_{n+p}& =\sum _{i=0}^n{C}_n^i{B}_i\sum _{j=0}^p{j}^{n-i}\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\\ & =\sum _{j=0}^p\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\sum _{i=0}^n{C}_n^i{B}_i{j}^{n-i}\\ 又\because 对于\forall & j\in (1,p)且p\in primes,有\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\equiv 0(modp)\\ \therefore B_{n+p}modp& =\sum _{j=0}^p\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\sum _{i=0}^n{C}_n^i{B}_i{j}^{n-i}modp\\ & =(\sum _{j=0}^p\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}modp)(\sum _{i=0}^n{C}_n^i{B}_i{j}^{n-i}modp)modp\\ & =(\left\{\begin{array}{c}p\\ \\ 1\end{array}\right\}+0+\left\{\begin{array}{c}p\\ \\ p\end{array}\right\})(\sum _{i=0}^n{C}_n^i{B}_i{j}^{n-i}modp)modp\\ & =(\left\{\begin{array}{c}p\\ \\ 1\end{array}\right\}(\sum _{i=0}^n{C}_n^i{B}_i{1}^{n-i}modp)modp+\left\{\begin{array}{c}p\\ \\ p\end{array}\right\}(\sum _{i=0}^n{C}_n^i{B}_i{p}^{n-i}modp)modp)modp\\ & =(\sum _{i=0}^n{C}_n^i{B}_i{1}^{n-i}+\sum _{i=0}^n{C}_n^i{B}_i{p}^{n-i})modp\\ & =(\sum _{i=0}^n{C}_n^i{B}_i+\sum _{i=0}^n{C}_n^i{B}_i{p}^{n-i})modp\\ & =(B_{n+1}+\sum _{i=0}^n{C}_n^i{B}_i{p}^{n-i})modp\\ 又\because 对\forall i>0,& p^i\equiv 0(modp)\\ \therefore \sum _{i=0}^n{C}_n^i& B_i{p}^{n-i}modp=C_n^n{B}_n{p}^{n-n}modp=B_{n}modp\\ \therefore B_{n+p}& \equiv B_{n+1}+B_n(modp)\\ 故：结论得证& \end{array}$$
TIPS：

证明：$对于\forall j\in (1,p)且p\in primes,有\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\equiv 0(modp)$
$$\begin{gathered} 由第二类斯特林数的通项公式得：\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}=\sum _{k=0}^m\frac{(-1{)}^k}{k!}\frac{(j-k{)}^p}{(j-k)!} \\ \begin{array}{rl}\therefore \left\{\begin{array}{c}p\\ \\ j\end{array}\right\}modp& =(\sum _{k=0}^j\frac{(-1{)}^k}{k!}\frac{(j-k{)}^p}{(j-k)!})modp\\ & =(\sum _{k=0}^j(\frac{(-1{)}^k}{k!}modp)(\frac{(j-k{)}^p}{(j-k)!}modp)modp)modp\end{array} \\ 由费马小定理可得:对于\forall p\in primes且\forall a满足gcd(a,p)=1，都有a^{p-1}\equiv 1(modp) \\ \therefore \frac{(j-k{)}^p}{(j-k)!}modp=\frac{(j-k)}{(j-k)!}modp=\frac{1}{(j-k-1)!}modp \\ \therefore \left\{\begin{array}{c}p\\ \\ j\end{array}\right\}modp=(\sum _{k=0}^j\frac{(-1{)}^k}{k!}\frac{1}{(j-k-1)!})modp \\ 这个时候不妨考虑生成函数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}x{e}^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,有\left\{\begin{array}{c}p\\ \\ j\end{array}\right\}\equiv 0(modp) \\ 结论得证 \end{gathered}$$

4、贝尔三角形

构造三角矩阵

满足一下递推式：
$$a_{n,m}=\{\begin{array}{l}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{array}$$

$$\begin{array}{rl}& 1\\ & 12\\ & 235\\ & 571015\\ & 1520273752\\ & 526787114151203\\ & 203255322409523674877\end{array}$$

每行的首项为贝尔数

5、快速求贝尔数在模意义下的第$n$项

题目链接：https://vjudge.csgrandeur.cn/problem/HDU-4767

题意：求$B_{n}mod95041567$

思路：

​ 首先，我们发现$95041567$并不是一个素数；但是我们发现$95041567$是多个不相同的素数的乘积，因此我们可以考虑CRT合并，所以我们可以列出以下方程：
$$\begin{gathered} 95041567={\Pi }_i{p}_i \\ \{\begin{array}{l}x\equiv a_1(mod{p}_1)\\ x\equiv a_2(mod{p}_2)\\ ....\\ x\equiv a_n(mod{p}_k)\end{array} \end{gathered}$$
​ 因此我们将问题转化为了求出一组$a_i=B_n \space mod \space p_i\$

​ 而对于每一组素数$p_i$，我们可以利用$Touchard同余$
$$\begin{gathered} 将B_{n+p}\equiv B_{n+1}+B_n(modp)改写成模意义下的矩阵相乘的形式: \\ \left[\begin{array}{c}{B}_{n+p-1}\\ B_{n+p}\\ B_{n+p+1}\\ ⋮\\ B_{n+p+p-1}\end{array}\right]=\left[\begin{array}{ccccccc}& 0& 0& 0& \dots & 0& 1\\ & 1& 1& 0& \dots & 0& 0\\ & 0& 1& 1& \dots & 0& 0\\ & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ & 0& 0& 0& \dots & 1& 1\end{array}\right]\times \left[\begin{array}{c}{B}_n\\ B_{n+1}\\ B_{n+2}\\ ⋮\\ B_{n+p-1}\end{array}\right] \end{gathered}$$
那么可以利用矩阵的递推推出第$n$项与前几项的关系：
$$\left[\begin{array}{c}{B}_n\\ B_{n+1}\\ B_{n+2}\\ ⋮\\ B_{n+p-1}\end{array}\right]={\left[\begin{array}{ccccccc}& 0& 0& 0& \dots & 0& 1\\ & 1& 1& 0& \dots & 0& 0\\ & 0& 1& 1& \dots & 0& 0\\ & ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ & 0& 0& 0& \dots & 1& 1\end{array}\right]}^{\frac{n}{p-1}}\times \left[\begin{array}{c}{B}_0\\ B_1\\ B_2\\ ⋮\\ B_{p-1}\end{array}\right]$$
我们发现中间那个矩阵可以利用矩阵快速幂求出，那么我们只需要预处理出前$p_i$个$B_i$，然后利用矩阵乘法就可以得到$B_{n}mod{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;
}