HDUOJ 4767 Bell 长春网络赛1009

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4767

Bell 数背景知识：http://zh.wikipedia.org/wiki/%E8%B4%9D%E5%B0%94%E6%95%B0

当时比赛临时google的，发现有一个递推式可以用：B[n+p] = B[n] +B[n+1] (mod p)，然后又因为题目中给的大数可以分解为：31*37*41*43*47

所以我们可以对这5个数分别求出第n项模这个质数的结果，（构造31 37 41 43 47 维矩阵）矩阵方法需要前50项（大概数目，每个模不一样，所以我用含组合数的递推公式求了bell数前50的模31 37 41 43 47的结果，组合数当然也要分别求5组关于5个模的结果）最后利用中国剩余定理解决，代码细节如下：

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<iostream>
#define LL int
#define ll int 
#pragma warning (disable : 4996)
using namespace std;
#define N 50
#define mem(a) memset(a, 0, sizeof(a))
#define Mod 95041567

LL a[N], m[N] = {31, 37, 41, 43, 47}, bell[N][5], c[N][N][5];

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 -= (a / b) * x;
}

LL china (ll n, int *a, int *m){
    LL re = 0, i, j, k, M = 1, d, x, y;
    for (i = 0; i < n; ++i) M *= m[i];
    for (i = 0; i < n; ++i){
        exgcd(M / m[i], m[i], d, x, y);
        re = (re + x * M / m[i] * a[i]) % M;
    }
    return (re + M) % M;
}

struct Matrix {
    int v[N][N];
};

Matrix ori, ans, I;

Matrix mul(Matrix m1,Matrix m2, LL M){
    int i,j,k;
    Matrix c;
    mem(c.v);
    for(i=0;i<N;i++)for(j=0;j<N;j++){
        for(k=0;k <N;k++)
            c.v[i][j]+=(m1.v[i][k]*m2.v[k][j]) % M;
    }
    return c;
}

Matrix Mpow(Matrix A, int n, LL M){
    Matrix x=A,c=I;
    while(n>=1){
        if(n&1) c=mul(c,x,M);
        n>>=1;
        x=mul(x,x,M);
    }
    return c;
}

int main(){
    /*freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);*/
    LL n, i, j, k, ca = 1, re, t;
    for (i = 0; i < 5; ++i) bell[0][i] = bell[1][i] = 1;
    for (i = 0; i < N; ++i) for (k = 0; k < 5; ++k) c[i][0][k] = 1, c[i][i][k] = 1;
    for (i = 1; i < N; ++i)
        for (j = 1; j < i; ++j) for (k = 0; k < 5; ++k) c[i][j][k] = (c[i - 1][j - 1][k] + c[i - 1][j][k]) % m[k];
    for (i = 2; i < N; ++i)
        for (k = 0; k < 5; ++k)
        for (j = 0; j < i; ++j) 
            bell[i][k] = (bell[i][k] + bell[j][k] * c[i - 1][j][k]) % m[k];
    for (i = 0; i < N; ++i) I.v[i][i] = 1;
    cin >> t;
    while (t--){
        cin >> n;
        if (n < 50) {
            for (i = 0; i < 5; ++i) a[i] = bell[n][i];
            re = china(5, a, m);
            printf("%d\n", re);
        }
        else{
			for (j = 0; j < 5; ++j) {
				mem(ori.v);
				for (i = 0; i < m[j] - 1; ++i) ori.v[i + 1][i] = 1;
				ori.v[0][m[j] - 1] = ori.v[1][m[j] - 1] = 1;
				ans = Mpow(ori, n - m[j] + 1, m[j]);
				 for (i = 0, a[j] = 0; i < m[j]; ++i) a[j] += ((bell[i][j] * ans.v[i][m[j] - 1]) % m[j]);
				 a[j] %= m[j];
			}
            re = china (5, a, m);
            printf("%d\n", re);
        }
    }
    
    return 0;
}

哎总之比赛时能过略幸运，因为这题debug占了好多时间，结果队友没法拍题了，略感亚历山大，还好最后AC了。