HDU4767 Bell 中国剩余定理 贝尔数 第二类斯特灵数

这道题要是出在现场赛 估计做出的人要少很多很多，要查阅很多资料才做的出来，至少对于我这个菜狗来说 必须要查资料

http://zh.wikipedia.org/wiki/%E6%96%AF%E7%89%B9%E7%81%B5%E6%95%B0

首先要了解第二类斯特灵数，还有贝尔数，

接下来是解题的关键：

若p为任意质数，则有Bell(n+p) = Bell(n) + Bell(n+1) (mod  p)；

我们来分解题目给的模的数 95041567=31*37*41*43*47，这五个数都是素数，所以可以用上述的式子来解题目，同时又因为 这五个数两两互素，所以可以利用中国剩余定理来解 x=a1(mod m1)这个同余方程组，求出 mod（m1*m2...*mi）d的值，这样就可以得到答案了

但是在使用中国剩余定理之前 我们要求出 Bell（n）%各个素数的 值，这里需要矩阵快速幂

举个例子 若 素数为5此时可以构造这样一个矩阵

0 1 0 0 0
0 1 1 0 0
0 0 1 1 0
0 0 0 1 1
1 0 0 0 1

再列举一个3的矩阵，你就会发现一定的规律，这样就可以求出 Bell(n)%各个素数的值了

#include<iostream>#include<cstdio>#include<list>#include<algorithm>#include<cstring>#include<string>#include<queue>#include<stack>#include<map>#include<vector>#include<cmath>#include<memory.h>#include<set>#define ll long long#define LL __int64#define eps 1e-8//const ll INF=9999999999999;#define inf 0xfffffffusing namespace std;//vector<pair<int,int> > G;//typedef pair<int,int> P;//vector<pair<int,int>> ::iterator iter;////map<ll,int>mp;//map<ll,int>::iterator p;////vector<int>G[30012];int MOD[5]={31,37,41,43,47};int a[12][62][62],b[12][62];int c[62];int n;void getstrling()//获取第二类斯特灵数，可以看维基百科{for(int i=0;i<5;i++)a[i][0][0]=1;for(int i=1;i<=50;i++){for(int j=0;j<5;j++)a[j][i][1]=1;for(int j=1;j<=i;j++){for(int k=0;k<5;k++)a[k][i][j]=(a[k][i-1][j-1]+j*a[k][i-1][j])%MOD[k];}for(int j=0;j<5;j++){b[j][i]=0;for(int k=1;k<=i;k++)b[j][i]=(b[j][i]+a[j][i][k])%MOD[j];}}}int Matrixquick(int id,int mod){int x[62],xx[62],y[62][62],yy[62][62];int len=MOD[id]+1;for(int i=1;i<=len+1;i++)for(int j=1;j<=len+1;j++)y[i][j]=0;for(int i=1;i<len;i++)y[i+1][i]=1;y[2][len]=y[3][len]=1;for(int i=1;i<=len;i++)x[i]=b[id][i];//构造矩阵的过程，这个矩阵有点规律的，if(MOD[id]>=n)return x[n];int tempn=n-1;while(tempn)//再写一个函数反而觉得麻烦 直接讲快速幂套在里面了{if(tempn&1){for(int i=1;i<=len;i++){xx[i]=0;for(int j=1;j<=len;j++)xx[i]+=x[j]*y[j][i];}for(int i=1;i<=len;i++)x[i]=xx[i]%mod;}for(int i=1;i<=len;i++){for(int j=1;j<=len;j++){yy[i][j]=0;for(int k=1;k<=len;k++)yy[i][j]+=y[i][k]*y[k][j];yy[i][j]%=mod;}}for(int i=1;i<=len;i++)for(int j=1;j<=len;j++)y[i][j]=yy[i][j];tempn>>=1;}return x[1];}int exgcd(int a,int &x,int b,int &y){if(b==0){x=1;y=0;return a;}int r=exgcd(b,x,a%b,y);int tmp=x;x=y;y=tmp-a/b*y;return r;}int china()//中国剩余定理{int M=1;int ans=0;for(int i=0;i<5;i++)M*=MOD[i];for(int i=0;i<5;i++){int Mi=M/MOD[i];int x0,y0;int gcd=exgcd(Mi,x0,MOD[i],y0);ans=(ans+x0*Mi*c[i])%M;}return (ans+M)%M;}int main(void){getstrling();int t;scanf("%d",&t);while(t--){scanf("%d",&n);for(int i=0;i<5;i++)c[i]=Matrixquick(i,MOD[i]);int ans=china();printf("%d\n",ans);}}