题意
设第 $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 \times 37 \times 41 \times 43 \times 47$，所以预处理前50个，
对于 $n > 50$，使用递推式，递推式可转成矩阵乘法，如下：
$$\begin{bmatrix} 0 & 0 & \cdots & 1\\
1 & 1 & \cdots & 0\\ \vdots & \vdots &\ddots & \vdots\\ 0 & \cdots & 1 & 1 \end{bmatrix} \times
\begin{bmatrix} B_n\\ B_{n+1}\\ \vdots\\ B_{n+p-1} \end{bmatrix} =
\begin{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(5\cdot 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\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 + 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("\n");
int ans = china(5, res, m);
printf("%d\n", 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\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;
}
参考链接：
1. https://zh.wikipedia.org/w/index.php?title=%E8%B4%9D%E5%B0%94%E6%95%B0
2. https://www.cnblogs.com/yuyixingkong/p/4489189.html
3.https://www.cnblogs.com/Chierush/p/3344661.html
原文地址：https://www.cnblogs.com/lfri/p/11546313.html
时间： 2024-11-13 08:14:43