各类数（四）—— Bell 数

最新推荐文章于 2022-04-24 20:28:51 发布

lxt_Lucia

最新推荐文章于 2022-04-24 20:28:51 发布

阅读量1.6k

点赞数

分类专栏： algorithm knowledge 文章标签：各类数 Bell数

knowledge 同时被 2 个专栏收录

42 篇文章 24 订阅

订阅专栏

algorithm

26 篇文章 28 订阅

订阅专栏

欢迎访问https://blog.csdn.net/lxt_Lucia～～

宇宙第一小仙女\(^o^)/～～萌量爆表求带飞=≡Σ((( つ^o^)つ~ dalao们点个关注呗～～

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 <cstdio>
#include <stdlib.h>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <string.h>
#include <vector>
#include <queue>
#include <stack>
#include <set>
#include <map>
#include <ctime>
#define maxn 1007
//#define N 100005
#define INF 0x3f3f3f3f
#define PI acos(-1)
#define lowbit(x) (x&(-x))
#define eps 0.000000001
#define read(x) scanf("%d",&x)
#define put(x) printf("%d\n",x)
#define Debug(x) cout<<x<<" "<<endl
using namespace std;
typedef long long ll;
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;
}

转自：https://blog.csdn.net/acdreamers/article/details/12309269