题意
设第 $n$ 个Bell数为 $B_n$，求 $B_n mod 95041567$.（$1 leq n leq 2^{31}$）
分析
贝尔数的概念和性质，维基百科上有，这里用到两点。
- 若 $p$ 是任意素数，有 $B_{p+n} = B_n + B_{n+1}(mod p)$
- 每个贝尔数都是相应第二类斯特林数之和，即 $B_n = sum_{k=1}^nS(n, k)$
贝尔数的这个递推式只适合质数，我们可以将模数拆成质数，然后用CRT合并。
$95041567 = 31 imes 37 imes 41 imes 43 imes 47$，所以预处理前50个，
对于 $n > 50$，使用递推式，递推式可转成矩阵乘法，如下：
$$egin{bmatrix} 0 & 0 & cdots & 1\
1 & 1 & cdots & 0\ vdots & vdots &ddots & vdots\ 0 & cdots & 1 & 1 end{bmatrix} imes
egin{bmatrix} B_n\ B_{n+1}\ vdots\ B_{n+p-1} end{bmatrix} =
egin{bmatrix} B_{n+p-1}\ B_{n+p}\ vdots \ B_{n+2p-2} end{bmatrix}$$
即 $B_{n+p-1} = A B_n$
设 $t = n / (p-1), k = n \% (p-1)$，
如果利用 $B_n = A^tB_k$，需要多预处理一倍，但计算时只需求第一个元素；
若利用 $B_{(p-1)t} = A^t B_0$，只需预处理前 $p-1$ 个，但是计算时需要算出第 $k$ 个。
反正两者时间也几乎一样。
时间复杂度为 $O(5cdot p^3 log{frac{n}{p}})$
#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[2*maxn][2*maxn], bell[5][2*maxn]; //第二类司特林数、贝尔数 void Stirling2() { Sti[0][0] = 1; for(int i = 0;i <= 2*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 <= 2*m[i];j++) { bell[i][j] = 0; //不知道为什么默认不是0 for(int k = 1;k <= j;k++) bell[i][j] = (bell[i][j] + Sti[j][k]) % m[i]; //printf("%d %d ",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 + tmp.mat[0][i] * bell[k][(n%(p-1))+i]) % p; return ret;} //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(" "); int ans = china(5, res, m); printf("%d ", ans); } return 0; }
第二种种写法：
#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 %d ",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(" "); int ans = china(5, res, m); printf("%d ", ans); } return 0; }
另外一种方法是利用公式：
$$B_{p^m+n} = mB_n + B_{n+1}$$
于是，我们求 $Bell(n)(mod p)$ 时，先把 $n$ 写成 $p$ 进制，即
$$n = a_kp^k +...+a_1p+a_0$$
先预处理 $i < p, bell[i] = Bell(i)(mod p)$，对于大于 $p$ 的用公式。
#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 < maxn;j++) { bell[i][j] = 0; //不知道为什么默认不是0 for(int k = 1;k <= j;k++) bell[i][j] = (bell[i][j] + Sti[j][k]) % m[i]; //printf("%d %d ",j,bell[i][j]); } } } int solve(int n, int p, int k) //计算Bn % p { int a[maxn], B[maxn], d[maxn]; for(int i = 0;i <= p;i++) B[i] = bell[k][i]; int cnt = 0; while(n) { a[cnt++] = n % p; n /= p; } for(int i = 1;i < cnt;i++) for(int j = 1;j <= a[i];j++) { for(int k = 0;k < p;k++) d[k] = (B[k+1] + i*B[k]) % p; d[p] = (d[0] + d[1]) % p; for(int k = 0;k <= p;k++) B[k] = d[k]; } //printf("after "); return d[a[0]]; } //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); if(n < maxn) { for(int i = 0;i < 5;i++) res[i] = bell[i][n]; } else { 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(" "); int ans = china(5, res, m); printf("%d ", ans); } return 0; }
这种方法的复杂度大约为 $O(p^2)$，比前面的方法快了不少。
%%%Acdreamers，链接
参考链接：
1. https://zh.wikipedia.org/w/index.php?title=%E8%B4%9D%E5%B0%94%E6%95%B0