贝尔数

题意：

给你两个公式：

对于任意质数p:

让你求Bell_n模95041567，n<=2147483647；

其中95041567=31*37*41*43*47;

考虑这些质数都很小，那么我们可以考虑对每一个质数分别处理，然后用CRT合并；

然后对于每个质数，我们发现第二个公式是一个类似斐波那契的递推，我们考虑用一个50*50的矩阵乘法优化递推；

下面给出CRT合并的公式(蒯自ACdreamer)：

设正整数两两互素，则同余方程组

 

                            

 

有整数解。并且在模下的解是唯一的，解为

 

                              

 

其中，而为模的逆元。

 

//MADE BY QT666
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=100050;
const int Mod=95041567;
int m[]={31,37,41,43,47};
ll c[55][55],bell[55],ans[3][55][55],s[55][55],a[N];
void work(int a,int b,int mo){
    for(int i=1;i<=50;i++){
	for(int j=1;j<=50;j++){
	    for(int k=1;k<=50;k++){
		(s[i][j]+=(ans[a][i][k]*ans[b][k][j])%mo)%=mo;
	    }
	}
    }
    for(int i=1;i<=50;i++){
	for(int j=1;j<=50;j++){
	    ans[a][i][j]=s[i][j],s[i][j]=0;
	}
    }
}
ll calc(int k,int y){
    memset(ans,0,sizeof(ans));
    memset(bell,0,sizeof(bell));
    memset(c,0,sizeof(c));
    for(int i=0;i<=50;++i) c[i][0]=1;
    for(int i=1;i<=50;++i)
	for(int j=1;j<=i;++j){
	    c[i][j]=(c[i-1][j-1]+c[i-1][j])%m[k];
	}
    bell[0]=1;
    for(int i=1;i<=50;i++){
	for(int j=0;j<i;j++) (bell[i]+=(bell[j]*c[i-1][j])%m[k])%=m[k];
    }
    for(int i=1;i<=50;i++) ans[1][1][i]=bell[i];
    ans[2][50-m[k]+1][50]=1,ans[2][50-m[k]+2][50]=1;
    for(int i=1;i<=49;i++) ans[2][i+1][i]=1;
    if(y<=50) return ans[1][1][y];
    else{
	y-=50;
	while(y){
	    if(y&1) work(1,2,m[k]);
	    work(2,2,m[k]);y>>=1;
	}
	return ans[1][1][50];
    }
}
void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b) {x=1,y=0;return;}
    exgcd(b,a%b,y,x),y-=x*(a/b);
}
ll qpow(ll x,ll y,ll mo){
    ll ret=1;
    while(y){
	if(y&1) (ret*=x)%=mo;
	(x*=x)%=mo;y>>=1;
    }
    return ret;
}
void Crt(int n){
    ll Ans=0;
    for(int i=0;i<5;i++){
	ll M=Mod/m[i];
	(Ans+=qpow(M,m[i]-2,m[i])*M%Mod*calc(i,n)%Mod)%=Mod;
    }
    printf("%lld\n",Ans);
}
int main(){
    freopen("bell.in","r",stdin);
    freopen("bell.out","w",stdout);
    int T;scanf("%d",&T);
    while(T--){
	int n;scanf("%d",&n);Crt(n);
    }
    return 0;
}