思路：矩阵快速幂+中国剩余定理！！
查资料得到2个公式：
1) B[n+p] = B[n] + B[n+1] mod p ;
2) B[p^m+n] = m*B[n] + B[n+1] mod p .
用这两个都可以解决这个问题，第二个可以在0ms解决。
质因数分解95041567=31*37*41*43*47，用矩阵快速幂分别求出B[n]%p (p是95041567的质因子)的结果。
这样就得到5个同余等式，在用中国剩余定理求解既可。
代码如下：
1 #include <iostream> 2 #include <cstring> 3 #include <cstdio> 4 #define M 51 5 #define mod 95041567 6 #define ll __int64 7 using namespace std; 8 int p[5]={31,37,41,43,47},a[5],bell[M][M]; 9 struct Mat 10 { 11 int m[M][M]; 12 }; 13 Mat Mul(Mat a,Mat b,int n) 14 { 15 Mat ans; 16 for(int i=0;i<n;i++) 17 for(int j=0;j<n;j++){ 18 ans.m[i][j]=0; 19 for(int k=0;k<n;k++) 20 ans.m[i][j]+=a.m[i][k]*b.m[k][j]; 21 ans.m[i][j]%=n; 22 } 23 return ans; 24 } 25 int Pow(int n,int m) 26 { 27 Mat ans,a; 28 memset(ans.m,0,sizeof(ans.m)); 29 memset(a.m,0,sizeof(a.m)); 30 a.m[0][m-1]=a.m[1][m-1]=1; 31 for(int i=0;i<m-1;i++) a.m[i+1][i]=1; 32 for(int i=0;i<m;i++) ans.m[0][i]=bell[i][i]%m; 33 while(n){ 34 if(n&1) ans=Mul(ans,a,m); 35 n>>=1; 36 a=Mul(a,a,m); 37 } 38 return ans.m[0][0]; 39 } 40 void init() 41 { 42 int i,j; 43 bell[0][0]=1;bell[1][1]=1; 44 for(i=2;i<M;i++){ 45 bell[i][1]=bell[i-1][i-1]; 46 for(j=2;j<=i;j++) bell[i][j]=(bell[i][j-1]+bell[i-1][j-1])%mod; 47 } 48 } 49 void gcd_extend(int a,int b,ll &x,ll &y) 50 { 51 if(b==0){ 52 x=1;y=0; 53 return ; 54 } 55 gcd_extend(b,a%b,x,y); 56 ll t=x; 57 x=y; 58 y=t-a/b*y; 59 } 60 int china() 61 { 62 ll ans=0,x,y; 63 for(int i=0;i<5;i++){ 64 int t=mod/p[i]; 65 gcd_extend(t,p[i],x,y); 66 ans=(ans+a[i]*t*((x%p[i]+p[i])%p[i]))%mod; 67 } 68 return ans; 69 } 70 int main(){ 71 int t,n; 72 init(); 73 scanf("%d",&t); 74 while(t--){ 75 scanf("%d",&n); 76 for(int i=0;i<5;i++) 77 a[i]=Pow(n,p[i]); 78 printf("%d ",china()); 79 } 80 return 0; 81 }