贝尔数 hdu4767 (矩阵快速幂+中国剩余定理+bell数+Stirling数+欧几里德)

#include <cstdio>
#include <cstring>
typedef long long llt ;
typedef long long int_t;
using namespace std;
const int MOD  = 95041567;
/*
 * 先写好三个模板, 矩阵快速幂 , 扩展欧几里德 , 中国剩余定理
 */
//非递归的扩展欧几里德算法
//返回a、b的gcd，同时x、y满足ax+by=gcd
int_t exEuclid(int_t a,int_t b,int_t&x,int_t&y){
    int_t x0 = 1, y0 = 0;
    int_t x1 = 0, y1 = 1;
    x = 0; y = 1;
    int_t r = a % b;
    int_t q = ( a - r ) / b;
    while( r ){
        x = x0 - q * x1;
        y = y0 - q * y1;
        x0 = x1; y0 = y1;
        x1 = x; y1 = y;
        a = b; b = r; r = a % b;
        q = ( a - r ) / b;
    }
    return b;
}
//求a相对于p的逆元，a、p互质才存在逆元
int_t inv(int_t a,int_t p){
    int_t x,y;
    int_t r = exEuclid(a,p,x,y);
    if ( r != 1 ) return 0;
    x = x % p;
    if ( x < 0 ) x += p;
    return x;
}

const int Cube_SIZE = 50;//矩阵大小
struct Cube{
    llt mat[Cube_SIZE][Cube_SIZE];
};
//单位矩阵
Cube _UnitCube;
//矩阵乘机
Cube _Multiply(Cube A,Cube B,llt mod){
    Cube _tmpCube;
    for (int i = 0;i < Cube_SIZE;++i)
    for (int j = 0;j < Cube_SIZE;++j){
        _tmpCube.mat[i][j] = 0;
        for (int k = 0;k < Cube_SIZE;++k){
            _tmpCube.mat[i][j] += A.mat[i][k]*B.mat[k][j];
            _tmpCube.mat[i][j] %= mod;
        }
    }
    return _tmpCube;
}
//矩阵快速幂
Cube power_Cube(Cube A,llt n,llt mod){
    Cube _tmpCube = _UnitCube;
    while(n){
        if(n & 1)
            _tmpCube = _Multiply(_tmpCube, A ,mod);
        A = _Multiply(A, A, mod);
        n >>= 1;
    }
    return _tmpCube;
}

// 中国剩余定理
int China( int n ,int  *remain, int *mod){

    llt M , x = 0;
    M = MOD;
    for ( int i = 0;i < 5;++i ){
        M /= mod[i];
        x = (x + remain[i]*M*inv(M,mod[i])) % MOD;
        M *= mod[i];
    }
    return x;
}


const int MOD = 95041567;
// 构建的一次矩阵
Cube BellCube;
//先计算出前50项
llt Bell[50];
llt striling2[50][50];

void GetStriling(){
    striling2[0][0] = 1;
    for(int i = 1;i < 50;++i ){
        striling2[i][0] = 0;
        for(int j = 1;j <= i; ++j)
            striling2[i][j] = ( striling2[i-1][j-1] + striling2[i-1][j]*j) % MOD;
    }
}
void GetBell(){
    memset(Bell,0,sizeof(Bell));
    Bell[0] = 1;
    for (int i = 1;i < 50;++i )
        for (int j = 0;j <= i;++j )
            Bell[i]  = ( Bell[i] + striling2[i][j] ) % MOD;
}

void init(){
    //得到bell数
    GetStriling();
    GetBell();
    //初始化单位矩阵
    memset(&_UnitCube, 0 ,sizeof(_UnitCube));
    for (int i = 0;i < 50;++i)
        _UnitCube.mat[i][i] = 1;

    //初始化一次矩阵
    memset(&BellCube , 0,sizeof(BellCube));
    for (int i = 1;i < 50 ;++i)
        BellCube.mat[i][i] = BellCube.mat[i][i-1] = 1;
}

// --------- 核心代码 ----------------//
// MOD的5个质因子
int factor[5] = { 31,37,41,43,47 };
// 用于记录5个余数
int leave[5];

int solve (int n ) {
    if ( n < 50 ) return Bell[n];
    llt ans[50];

    for (int i = 0;i < 5;++i ){
        Cube arr = BellCube;
        arr.mat[0][factor[i]-1] = 1;

        arr = power_Cube(arr,n/(factor[i]-1) , factor[i]);
        memset(ans , 0, sizeof(ans));

        for (int k = 0;k < factor[i];++k )
            for (int j = 0;j < factor[i];++j )
                ans[k] = (ans[k] + Bell[j]*arr.mat[k][j] ) % factor[i];

        leave[i] = ans[n%(factor[i]-1)];
    }
    // 利用5个因子与5个余数 求得 bell(n) 
    return China( n , leave , factor );
}

int main(){
    init();
    int t , n ;
    scanf("%d",&t);
    while ( t-- ) {
        scanf("%d",&n );
        printf("%d\n",solve(n));
    }
    return 0;
}