hdu477 Bell——求第n项贝尔数

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
const int maxn = 50;
const int mod = 95041567;

int m[5] = {31, 37, 41, 43, 47};
int Sti[maxn][maxn], bell[5][maxn];  //第二类司特林数、贝尔数

void Stirling2()
{
    Sti[0][0] = 1;
    for(int i = 0;i <= maxn;i++)
        for(int j = 1;j <= i;j++)
            Sti[i][j] = (Sti[i-1][j-1] + 1LL * j * Sti[i-1][j]) % mod;
}

void init()
{
    Stirling2();
    for(int i = 0;i < 5;i++)
    {
        bell[i][0] = 1;
        for(int j = 1;j <= m[i];j++)
        {
            bell[i][j] = 0;
            for(int k = 1;k <= j;k++)
                bell[i][j] = (bell[i][j] + Sti[j][k]) % m[i];
            //printf("%d\t%d\n",j,bell[i][j]);
        }
    }
}

struct matrix
{
    int r, c;
    int mat[maxn][maxn];
    matrix(){
        memset(mat, 0, sizeof(mat));
    }
};

matrix mul(matrix A, matrix B, int p)   //矩阵相乘
{
    matrix ret;
    ret.r = A.r; ret.c = B.c;
    for(int i = 0;i < A.r;i++)
        for(int k = 0;k < A.c;k++)
            for(int j = 0;j < B.c;j++)
            {
                ret.mat[i][j] = (ret.mat[i][j] + 1LL * A.mat[i][k] * B.mat[k][j]) % p;
            }
    return ret;
}

matrix mpow(matrix A, int n, int p)
{
    matrix ret;
    ret.r = A.r; ret.c = A.c;
    for(int i = 0;i < ret.r;i++)  ret.mat[i][i] = 1;
    while(n)
    {
        if(n & 1)  ret = mul(ret, A, p);
        A = mul(A, A, p);
        n >>= 1;
    }
    return ret;
}

int solve(int n, int p, int k)  //计算Bn % p
{
    matrix A;
    A.r = A. c = p;
    A.mat[0][p-1] = 1;
    for(int i = 1;i < p;i++)
        A.mat[i][i-1] = A.mat[i][i] = 1;

    matrix tmp = mpow(A, n/(p-1), p);

  //  int ret = 0;
//    for(int i = 0;i < p;i++)
//        ret = (ret + A.mat[0][i] * bell[k][i]) % p;

    int ans[p];
    for(int i = 0;i < p;i++)
    {
        ans[i] = 0;
        for(int j = 0;j < p;j++)
            ans[i] = (ans[i] + 1LL * bell[k][j] * tmp.mat[i][j]) % p;
    }


    return ans[n % (p-1)];
}

//ax + by = d，且|x|+|y|最小，其中d=gcd(a,b)
//即使a, b在int范围内，x和y也有可能超过int范围
void exgcd(ll a, ll b, ll &d, ll &x, ll &y)
{
    if (!b){ d = a; x = 1; y = 0;}
    else{ exgcd(b, a % b, d, y, x); y -= x * (a / b);}
}

//n个方程：x=a[i](mod m[i])
ll china(int n, int* a, int* m)
{
    ll M = 1, d, y, x = 0;
    for (int i = 0; i < n; i++)  M *= m[i];
    for (int i = 0; i < n; i++)
    {
        ll w = M / m[i];
        exgcd(m[i], w, d, d, y);        //d共用了
        x = (x + y * w * a[i]) % M;   //x相当于sum
    }
    return (x + M) % M;
}

int n;
int res[5];

int main()
{
    init();

    int T;
    scanf("%d", &T);
    while(T--)
    {
        scanf("%d", &n);
        for(int i = 0;i < 5;i++)  res[i] = solve(n, m[i], i);

//        for(int i = 0;i < 5;i++)  printf("%d ", res[i]);
//        printf("\n");

        int ans = china(5, res, m);
        printf("%d\n", ans);
    }
    return 0;
}