HDU 4767( china + 矩阵快速幂)

本题目的模数可以拆分啊，95041567 = 31*37*41*43*47；

有bell数有 对任意素数 b[ p+n ] = (b[ n ] + b[ n+1 ]) %p; 这个先用快速幂求出每个小素数模下的值，然后用China搞定。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
typedef long long ll;
typedef long long LL;
using namespace std;
#define rep1(i,x,y) for(int i=x;i<=y;i++)
#define repd(i,x,y) for(int i=x;i>=y;i--)
#define rep(i,n) for(int i=0;i<(int)n;i++)

const int MOD = 95041567;
const int N = 100;
int mod[]={31,37,41,43,47};
int c[5][50][50],A[5];
ll ans[5][50];
void init(){
   rep(k,5){
      memset(c[k],0,sizeof(c[k]));
      c[k][0][0]=1;
      for(int i=0;i<50;i++){
          c[k][i][0]=c[k][i][i]=1;
          for(int j=1;j<i;j++)
             c[k][i][j]=(c[k][i-1][j]+c[k][i-1][j-1])%mod[k];
      }
   }
   rep(k,5){
       ans[k][0]=ans[k][1]=1;
       rep1(i,2,49){
          ans[k][i]=0;
          rep1(j,0,i-1)
             ans[k][i]=(ans[k][i]+c[k][i-1][j]*ans[k][i-1-j])%mod[k];
       }
   }
}
void gcd(ll a,ll b,ll& d,ll& x,ll& y){
   if(!b){d=a; x=1; y=0;}
   else { gcd(b,a%b,d,y,x); y-=x*(a/b); }
}
LL china(int n,int* a,int* m){
   LL M = 1,d,y,x=0;
   rep(i,n) M*=m[i];
   for(int i=0;i<n;i++){
      LL w=M/m[i];
      gcd(m[i],w,d,d,y);
      x = (x+y*w*a[i])%M;
   }
   return (x+M)%M;
}
int n, modd;
struct Matrix{
    int mat[N][N];
    void show(){
      rep(i,n) rep(j,n) {cout<<mat[i][j]<<" "; if(j==n-1) cout<<endl; }
    }
}a;
Matrix Matrix_mul(Matrix a, Matrix b){
    Matrix ret;
    memset(ret.mat, 0, sizeof(ret.mat));
    for(int i = 0; i < n; i++)
    for(int j = 0; j < n; j++){
        if(a.mat[i][j]){
            for(int kk = 0; kk < n; kk++)
            ret.mat[i][kk] += (a.mat[i][j]*b.mat[j][kk])%modd;
        }
    }
    return ret;
}
Matrix Matrix_pow(Matrix tt, int nn){
    Matrix ret;
    Matrix temp = tt;
    memset(ret.mat, 0, sizeof(ret.mat));
    for(int i = 0; i< n; i++) ret.mat[i][i] = 1;
    while(nn){
        if(nn&1) ret = Matrix_mul(ret, temp);
        temp = Matrix_mul(temp, temp);
        nn >>= 1;
    }
    return ret;
}
int nn,te[N];
int main()
{
    init();
    int T;
    scanf("%d",&T);
    while(T--){
         scanf("%d",&nn);
         if(nn < 50){
            rep(i,5) A[i]=ans[i][nn];
            printf("%d\n",(int)china(5,A,mod));
         }
         else {
             rep(k,5){
                 n=modd=mod[k];
                 rep(i,mod[k]-1){
                     rep(j,n) a.mat[i][j] = 0;
                     a.mat[i][i+1]=1;
                 }
                 rep(i,mod[k]) a.mat[mod[k]-1][i] = 0;
                 a.mat[mod[k]-1][0]=a.mat[mod[k]-1][1]=1;
                 for(int i=0,j=1;i<n;i++,j++) te[i]=ans[k][j];
                 Matrix tt=Matrix_pow(a,nn-mod[k]);
                 A[k]=0;
                 rep(i,mod[k]) {A[k]=(A[k]+tt.mat[mod[k]-1][i]*te[i])%mod[k];}
             }
             printf("%d\n",(int)china(5,A,mod));
         }
    }
    return 0;
}