Bell
Time Limit:3000MS Memory Limit:32768KB 64bit IO Format:%I64d & %I64u
Description
What? MMM is learning Combinatorics!?
Looks like she is playing with the bell sequence now:
bell[n] = number of ways to partition of the set {1, 2, 3, ..., n}
Looks like she is playing with the bell sequence now:
bell[n] = number of ways to partition of the set {1, 2, 3, ..., n}
e.g. n = 3:
{1, 2, 3}
{1} {2 3}
{1, 2} {3}
{1, 3} {2}
{1} {2} {3}
so, bell[3] = 5
MMM wonders how to compute the bell sequence efficiently.
Input
T -- number of cases
for each case:
n (1 <= n < 2^31)
for each case:
n (1 <= n < 2^31)
Output
for each case:
bell[n] % 95041567
bell[n] % 95041567
Sample Input
6
1
2
3
4
5
6
Sample Output
1
2
5
15
52
203
首先，我很高兴终于把这题给解出来了！！！自己有几个难点一直没想通，今天无聊的时候一想。然后思路莫名奇妙的顺了，一切都变的理所当然。。。然后，我趁热打铁，把思路一理，就被我A了！
题意：就是求n个数的非空集合的划分方法数；
例如n = 3
{1, 2, 3}
{1} {2, 3}
{1, 2} {3}
{1, 3} {2}
{1} {2} {3}
所以Bell(3) = 5
给你一个n，求Bell(n) % 95041567的值
这道题涉及的知识比较多，，，本人觉得有一定的难度系数，没理解的，建议可以先从入门题刷起；
思路：
首先必须了解一个概念贝尔数（来自维基百科）& Stirling数（不懂的点链接）；
其次，如果你不懂Stirling数，那么自己动手推一下，其实这个还是简单的，杭电有一题简单的题，可以练习下：
http://acm.hdu.edu.cn/showproblem.php?pid=2512 相信很多同学都做过了，只是当时不知道是Stirling数，这里是第二类Stirling数；然后贝尔数呢，就是第二类Stirling数的和；
通过上面这个公式，我们可以利用矩阵计算大于n的贝尔数，方法是利用中国剩余定理，结合欧几里德和同余方程；首先我们输入n(假设，n>p，p是95041567的质因子，eg:n=50)，对n<50的Bell数进行预处理，构造一个p*p的矩阵；
利用公式求出，他们的余数存入数组at[i];质因子m[5]={31,37,41,43,47};利用中国剩余定理；求得余数！
转载请注明出处：![Bell(hdu4767+矩阵+中国剩余定理+bell数+Stirling数+欧几里德)-LMLPHP]( "Bell(hdu4767+矩阵+中国剩余定理+bell数+Stirling数+欧几里德)-LMLPHP")
寻找&星空の孩子
题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4767
#include<stdio.h>
#include<string.h>
#define LL long long
#define mod 95041567
#define MM 50 int m[]={,,,,};
int Sti[MM][MM],bell[MM][MM];
int at[];//各项余数 void stirling2()
{
memset(Sti,,sizeof(Sti));
Sti[][]=;
for(int i=;i<=MM;i++)
{
for(int j=;j<=i;j++)
{
Sti[i][j]=(int)(Sti[i-][j-]+((LL)j*(LL)Sti[i-][j])%mod)%mod;
}
}
}
void init()
{
stirling2();
for(int i=;i<;i++)
{
bell[i][]=;
for(int j=;j<m[i];j++)
{
bell[i][j]=;
for(int k=;k<=j;k++)
{
bell[i][j]=(bell[i][j]+Sti[j][k])%m[i];
}
// printf("%d\t%d\n",j,bell[i][j]);
}
}
}
struct Matrix
{
int mat[MM][MM];
}; Matrix multiply(Matrix a,Matrix b,int M)
{
Matrix c;
memset(c.mat,,sizeof(c.mat));
for(int i=;i<M;i++)
{
for(int j=;j<M;j++)
{
if(a.mat[i][j]==)continue;
for(int k=;k<M;k++)
{
if(b.mat[j][k]==)continue;
c.mat[i][k]=(c.mat[i][k]+a.mat[i][j]*b.mat[j][k])%M;
}
}
}
return c;
}
Matrix quickmod(Matrix a,int n,int M)
{
Matrix res; for(int i=;i<M;i++)
{
for(int j=;j<M;j++)
res.mat[i][j]=(i==j);
} while(n)
{
if(n&) res=multiply(res,a,M);
n>>=;
a=multiply(a,a,M);
}
return res;
} int work(int n,int M,int k)
{
Matrix per;//基础矩阵；
memset(per.mat,,sizeof(per.mat));
per.mat[][M-]=; for(int i=;i<M;i++)
{
per.mat[i][i]=per.mat[i][i-]=;
} Matrix tmp=quickmod(per,n/(M-),M); int ans[MM];
for(int i=;i<M;i++)
{
ans[i]=;
for(int j=;j<M;j++)
{
ans[i]=(ans[i]+bell[k][j]*tmp.mat[i][j])%M;
}
}
return ans[n%(M-)];
}
void exgcd(int a,int b,int& d,int& x,int &y)
{
if(!b){d=a;x=;y=;}
else
{
exgcd(b,a%b,d,y,x);
y-=x*(a/b);
}
}
int China(int r)
{
int Mc=;
int Mi,d,x,y,as=;
for(int i=;i<r;i++)
{
Mc*=m[i];
}
for(int i=;i<r;i++)
{
Mi=Mc/m[i];
exgcd(Mi,m[i],d,x,y);
as=(as+Mi*x*at[i])%Mc;
}
if(as<) as+=Mc;
return as;
}
int main()
{
int T,n;
scanf("%d",&T); init();
while(T--)
{
scanf("%d",&n);
for(int i=;i<;i++)
{
at[i]=work(n,m[i],i);
}
int sol=China();
printf("%d\n",sol);
} return ;
}