hdu 4767 bell 中国剩余定理＋矩阵快速幂

题目链接：点击打开链接

题意：求第n项bell数％95041567。

思路：

贝尔数有一些性质：点击打开链接

95041567可以分解质因数｛31，37，41，43，47｝，首先利用递推公式暴力求出前47项，然后利用同余公式矩阵快速幂求出bell(n)％31，37，41，43，47，之后用中国剩余定理合并。

cpp:点击打开链接

#include <cstdio>
#include <cstring>
using namespace std;
#define maxn 50
const int tp[5]={31,37,41,43,47};
int bell[maxn][5]={{1,1,1,1,1},{1,1,1,1,1},{2,2,2,2,2},{5,5,5,5,5},{15,15,15,15,15}};
struct Matrix {
	int n;
	int data[maxn][maxn];
	void set(int t,int tt){
		n=t;
		memset(data,0,sizeof(data));
		if(tt==1){
			for(int i=0;i<n;i++){
				data[i][i]=1;
			}
		}
	}
	void be_moded(int mod){
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				data[i][j]=data[i][j]%mod;
			}
		}
	}
	Matrix operator * (const Matrix& a){
		Matrix ret;
		ret.set(n,0);
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				for(int k=0;k<n;k++){
					ret.data[i][j]+=data[i][k]*a.data[k][j];
				}
			}
		}
		return ret;
	}
	void print(){
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				printf("%d ",data[i][j]);
			}
			puts("");
		}
	}
};
//矩阵快速幂
Matrix pow_mod(Matrix a,int n,int mod){
	Matrix ret;
	ret.set(a.n,1);
	while (n){
		if(n%2==1) ret=ret*a,ret.be_moded(mod);
		a=a*a;
		a.be_moded(mod);
		n/=2;
	}
	return ret;
}
//预处理出前47项bell数
void prepare(){
	int i;
	int c[maxn][maxn][5];
	memset(c,0,sizeof(c));
	for(int t=0;t<5;t++){
		c[1][0][t]=1;
		c[1][1][t]=1;
		for(i=2;i<maxn;i++){
			c[i][0][t]=1;
			for(int j=1;j<=i;j++){
				c[i][j][t]=(c[i-1][j-1][t]+c[i-1][j][t])%tp[t];
			}
		}
		for(i=4;i<maxn;i++){
			for(int j=0;j<=i;j++){
				bell[i+1][t]+=(c[i][j][t]*bell[j][t])%tp[t];
			}
		}
	}
}
//中国剩余定理
void china(int* f){
	long long A[5]={55185426,64217275,39407479,11051345,20221610};
	int ans=0;
	for(int i=0;i<5;i++){
		ans=(ans+(A[i]*f[i])%95041567)%95041567;	
	}
	printf("%d\n",ans);
}
//slove
void slove(int n){
	int ans[5];
	//求bell[n]%31,37,41,43,47
	for(int t=0;t<5;t++){
		//构造递推矩阵
		Matrix mat;
		mat.set(tp[t],0);
		for(int i=0;i+1<tp[t];i++){
			mat.data[i][i+1]=1;
		}
		mat.data[tp[t]-1][0]=1;
		mat.data[tp[t]-1][1]=1;
		int tpe=n-tp[t];
		if(tpe<=0){
			ans[t]=bell[n][t];
		}
		else {
			ans[t]=0;
			Matrix tp_mat=pow_mod(mat,tpe,tp[t]);
			for(int i=0;i<tp[t];i++){
				ans[t]+=(tp_mat.data[tp[t]-1][i]*bell[i+1][t]);
			}
		}
	}
	china(ans);
}
int main(){
	int T,n;
	prepare();
	scanf("%d",&T);
	while (T--){
		scanf("%d",&n);
		slove(n);
	}
	return 0;
}