Bell数

Bell数的定义：第n个Bell数表示集合{1,2,3,...,n}的划分方案数，即：B[0] = 1;

 

 

每一个Bell数都是第二类Stirling数的和，即：

 

 

第二类Stirling数的意义是：S(n,k)表示将n个物体划分成k个非空的不可辨别的（可以理解为盒子没有编号）集合的方法

数。很明显，每一个Bell是对应的第二类Stirling数之和。

 

Bell数的指数生成函数是：

 

 

关于它的推导过程详见：http://ftiasch.github.io/useless/posts/2013-09-27-generating-function-of-bell-number.html

 

 

 

Bell三角形(构建方法类似于杨辉三角形)

 

Bell三角形的构造方法：

第一行第一个元素是1，即a[1][1] = 1

对于n>1，第n行第一项等于第n-1行最后一项，即a[n][1] = a[n-1][n-1];

对于m,n>1，第n行第m项等于它左边和左上方的两个数之和，即a[n][m] = a[n][m-1] + a[n-1][m-1];

 

如图：

 

可以看出，每行首项是贝尔数，每行之和是第二类Stirling数。

 

 

Bell还有两个重要的同余性质：

 

 

 

其中这里的p是不大于100的素数，这样，我们可以通过上面的性质来计算Bell数模小于100的素数值。

 

Bell数模素数p的周期为：

 

 

 

典型题目：HDU2512，HDU4767

 

Bell数的预处理：

void Bell(int T[],int MOD){    B[0] = 1;    B[1] = 1;    T[0] = 1;    for(int i=2;i<N;i++)    {        T[i-1] = B[i-1];        for(int j=i-2;j>=0;j--)            T[j] = (T[j]+T[j+1])%MOD;        B[i] = T[0];    }}

 

题目：http://acm.hdu.edu.cn/showproblem.php?pid=4767

 

题意：给定一个数n，范围是[1,2^31]，求Bell(n)(mod 95041567)

 

分析：注意95041567 = 31x37x41x43x47，那么我们先对每一个素数求出Bell(n)(mod p)，然后CRT合并即可。

在这里，我们有两种方法，第一种方法就是用

 

 

所以可以根据它构造50*50的矩阵。如图：

 

 

 

#include <iostream>#include <string.h>#include <stdio.h>using namespace std;const int N = 50;int a[5];int m[5] = {31,37,41,43,47};int B[5][N],T[N];void Bell(int T[],int index,int MOD){    B[index][0] = 1;    B[index][1] = 1;    T[0] = 1;    for(int i=2;i<N;i++)    {        T[i-1] = B[index][i-1];        for(int j=i-2;j>=0;j--)            T[j] = (T[j]+T[j+1])%MOD;        B[index][i] = T[0];    }}struct Matrix{    int m[N][N];};Matrix I;Matrix multi(Matrix a,Matrix b,int MOD){    int i,j,k;    Matrix c;    for(i=0;i<N;i++)    {        for(j=0;j<N;j++)        {            c.m[i][j] = 0;            for(k=0;k<N;k++)                c.m[i][j] = (c.m[i][j] + a.m[i][k]%MOD * b.m[k][j]%MOD)%MOD;            c.m[i][j] %= MOD;        }    }    return c;}Matrix power(Matrix A,int k,int MOD){    Matrix ans = I,p = A;    while(k)    {        if(k&1)        {            ans = multi(ans,p,MOD);            k--;        }        k>>=1;        p = multi(p,p,MOD);    }    return ans;}void extend_Euclid(int a,int b,int &x,int &y){    if(b==0)    {        x = 1;        y = 0;        return;    }    extend_Euclid(b,a%b,x,y);    int tmp = x;    x = y;    y = tmp - (a/b)*y;}int CRT(int a[],int m[],int n){    int x,y;    int Mi,M=1,ans = 0;    for(int i=0;i<n;i++)        M *= m[i];    for(int i=0;i<n;i++)    {        Mi = M/m[i];        extend_Euclid(Mi,m[i],x,y);        ans = (ans + Mi*x*a[i])%M;    }    if(ans < 0) ans += M;    return ans;}void Init(){    for(int i=0;i<5;i++)       Bell(T,i,m[i]);    for(int i=0;i<N;i++)       for(int j=0;j<N;j++)           I.m[i][j] = (i == j);}void Work(int n){    if(n < N)    {        for(int i=0;i<5;i++)            a[i] = B[i][n];        int ans = CRT(a,m,5);        printf("%d\n",ans);    }    else    {        Matrix A;        for(int i=0;i<5;i++)        {            for(int k=0;k<N;k++)               for(int j=0;j<N;j++)                  A.m[k][j] = 0;            for(int j=0;j<m[i]-1;j++)               A.m[j+1][j] = 1;            A.m[0][m[i]-1] = A.m[1][m[i]-1] = 1;            Matrix ans = power(A,n - m[i] + 1,m[i]);            a[i] = 0;            for(int j=0;j<m[i];j++)                a[i] = (a[i] + B[i][j]%m[i]*ans.m[j][m[i]-1]%m[i])%m[i];            a[i] %= m[i];        }        int ret = CRT(a,m,5);        printf("%d\n",ret);    }}int main(){    Init();    int T;    scanf("%d",&T);    while(T--)    {        int n;        scanf("%d",&n);        Work(n);    }    return 0;}

 

 

第二种方法就是利用公式：

 

 

思路：我们求Bell(n)(mod p)，那么先把n写成p进制，即：

 

 

先预处理b[i] = Bell(i)(mod p)，然后对于大数n求Bell(n)(mod p)就很容易写出代码：

 

for i = 0 to p    c[i] = b[i];for i = 1 to cnt-1    begin    for j = 1 to a[i] do        begin        for k = 0 to p-1    d[k] = (c[k+1] + i*c[k]) mod p;        d[p] = (d[1] + d[0]) mod p;        for k = 0 to p    c[k] = d[k];end    end

那么，d[a[0]]的值就是Bell(n)(mod p)

 

#include <iostream>#include <string.h>#include <stdio.h>using namespace std;const int N = 50;int a[5];int m[5] = {31,37,41,43,47};int B[5][N],T[N];void Bell(int T[],int index,int MOD){    B[index][0] = 1;    B[index][1] = 1;    T[0] = 1;    for(int i=2;i<N;i++)    {        T[i-1] = B[index][i-1];        for(int j=i-2;j>=0;j--)            T[j] = (T[j]+T[j+1])%MOD;        B[index][i] = T[0];    }}int Solve(int n,int index){    int p = m[index];    int a[N],c[N],d[N];    for(int i=0;i<=p;i++)        c[i] = B[index][i];    int cnt = 0;    while(n)    {        a[cnt++] = n%p;        n /= p;    }    for(int i=1;i<cnt;i++)    {        for(int j=1;j<=a[i];j++)        {            for(int k=0;k<p;k++)               d[k] = (c[k+1] + i*c[k])%p;            d[p] = (d[0] + d[1])%p;            for(int k=0;k<=p;k++)               c[k] = d[k];        }    }    return d[a[0]];}void extend_Euclid(int a,int b,int &x,int &y){    if(b==0)    {        x = 1;        y = 0;        return;    }    extend_Euclid(b,a%b,x,y);    int tmp = x;    x = y;    y = tmp - (a/b)*y;}int CRT(int a[],int m[],int n){    int x,y;    int Mi,M=1,ans = 0;    for(int i=0;i<n;i++)        M *= m[i];    for(int i=0;i<n;i++)    {        Mi = M/m[i];        extend_Euclid(Mi,m[i],x,y);        ans = (ans + Mi*x*a[i])%M;    }    if(ans < 0) ans += M;    return ans;}void Init(){    for(int i=0;i<5;i++)       Bell(T,i,m[i]);}void Work(int n){    if(n<50)    {        for(int i=0;i<5;i++)           a[i] = B[i][n];        int ans = CRT(a,m,5);        printf("%d\n",ans);    }    else    {        for(int i=0;i<5;i++)           a[i] = Solve(n,i);        int ans = CRT(a,m,5);        printf("%d\n",ans);    }}int main(){    Init();    int T;    scanf("%d",&T);    while(T--)    {        int n;        scanf("%d",&n);        Work(n);    }    return 0;}