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 pc[i] = b[i]; for i = 1 to cnt-1beginfor j = 1 to a[i] dobeginfor k = 0 to p-1d[k] = (c[k+1] + i*c[k]) mod p;d[p] = (d[1] + d[0]) mod p;for k = 0 to pc[k] = d[k];endend

那么，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; }

?

總結

以上是生活随笔為你收集整理的Bell数的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。