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;
}