hdu 4767 Bell

思路：矩阵快速幂+中国剩余定理！！

查资料得到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个同余等式，在用中国剩余定理求解既可。

代码如下：

#include <iostream>
#include <cstring>
#include <cstdio>
#define M 51
#define mod 95041567
#define ll __int64
using namespace std;
int p[5]={31,37,41,43,47},a[5],bell[M][M];
struct Mat
{
    int m[M][M];
};
Mat Mul(Mat a,Mat b,int n)
{
    Mat ans;
    for(int i=0;i<n;i++)
    for(int j=0;j<n;j++){
        ans.m[i][j]=0;
        for(int k=0;k<n;k++)
            ans.m[i][j]+=a.m[i][k]*b.m[k][j];
        ans.m[i][j]%=n;
    }
    return ans;
}
int Pow(int n,int m)
{
    Mat ans,a;
    memset(ans.m,0,sizeof(ans.m));
    memset(a.m,0,sizeof(a.m));
    a.m[0][m-1]=a.m[1][m-1]=1;
    for(int i=0;i<m-1;i++) a.m[i+1][i]=1;
    for(int i=0;i<m;i++) ans.m[0][i]=bell[i][i]%m;
    while(n){
        if(n&1) ans=Mul(ans,a,m);
        n>>=1;
        a=Mul(a,a,m);
    }
    return ans.m[0][0];
}
void init()
{
    int i,j;
    bell[0][0]=1;bell[1][1]=1;
    for(i=2;i<M;i++){
        bell[i][1]=bell[i-1][i-1];
        for(j=2;j<=i;j++) bell[i][j]=(bell[i][j-1]+bell[i-1][j-1])%mod;
    }
}
void gcd_extend(int a,int b,ll &x,ll &y)
{
    if(b==0){
        x=1;y=0;
        return ;
    }
    gcd_extend(b,a%b,x,y);
    ll t=x;
    x=y;
    y=t-a/b*y;
}
int china()
{
    ll ans=0,x,y;
    for(int i=0;i<5;i++){
        int t=mod/p[i];
        gcd_extend(t,p[i],x,y);
        ans=(ans+a[i]*t*((x%p[i]+p[i])%p[i]))%mod;
    }
    return ans;
}
int main(){
    int t,n;
    init();
    scanf("%d",&t);
    while(t--){
        scanf("%d",&n);
        for(int i=0;i<5;i++)
            a[i]=Pow(n,p[i]);
        printf("%d
",china());
    }
    return 0;
}