题目链接: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了。