矩阵快速幂
快速幂
为了更快的进行幂运算,我们可以使用这个方法。原理太简单了,不解释。
看一下快速幂的非递归实现:
typedef long long ll;
/*
base 基数
p 指数
m 模数(看要求,决定要不要)
*/
ll fastPower(ll base, ll p, ll m)
{
ll res = 1;
while(p)
{
if(p & 1) res = res * base % m;
p >>= 1;
base = base * base % m;
}
return res;
}
矩阵快速幂
矩阵快速幂是为了加速矩阵的乘法运算。我在大一的线性代数课上学到了这个矩阵的相关知识,当时没觉得有什么用,现在发现...有点东西。
矩阵相乘的算法呢,百度一下有很多,我在这贴一个 OI Wiki 上的C++实现。
struct mat {
LL a[sz][sz];
inline mat() { memset(a, 0, sizeof a); }
inline mat operator-(const mat& T) const {
mat res;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) {
res.a[i][j] = (a[i][j] - T.a[i][j]) % MOD;
}
return res;
}
inline mat operator+(const mat& T) const {
mat res;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) {
res.a[i][j] = (a[i][j] + T.a[i][j]) % MOD;
}
return res;
}
inline mat operator*(const mat& T) const {
mat res;
int r;
for (int i = 0; i < sz; ++i)
for (int k = 0; k < sz; ++k) {
r = a[i][k];
for (int j = 0; j < sz; ++j)
res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
}
return res;
}
inline mat operator^(LL x) const {
mat res, bas;
for (int i = 0; i < sz; ++i) res.a[i][i] = 1;
for (int i = 0; i < sz; ++i)
for (int j = 0; j < sz; ++j) bas.a[i][j] = a[i][j] % MOD;
while (x) {
if (x & 1) res = res * bas;
bas = bas * bas;
x >>= 1;
}
return res;
}
};
那它有什么用呢,目前我知道的是,它可以加快一些递推公式的递推速度。
比如说呢,斐波那契数列:
f[0] = f[1] = 1,
f[n] = f[n - 1] + f[n - 2].
也就是$$a_0$$ = 1, $$a_1$$ = 1, $$a_n$$ = $$a_{n-1}$$ + $$a_{n - 2}$$。
如果我们要求$$a_n$$的话,我们需要从$$a_0$$ = 1, $$a_1$$ = 1开始递推到$$a_n$$,时间复杂度是O(N).但是如果我们用到矩阵加速的话,时间复杂度可以小很多具体多少额也不清楚。看看下面这个式子:
$$\begin{bmatrix}{F}_{2}&{F}_{3}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}\end{bmatrix} * \begin{bmatrix}0&1\\1&1\end{bmatrix}$$
它就是$$F_3 = F_1 + F_2$$这个过程,那就好说了。
$$\begin{bmatrix}{F}_{n}&{F}_{n + 1}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}\end{bmatrix} * {\begin{bmatrix}0&1\\1&1\end{bmatrix}}^{n - 1}$$
现在我们要求$$F_n$$的话,只需要对右边那个矩阵进行矩阵加速,得到它的n - 1次方即可。
我写了个模板:
#define MATRIXSIZE 2
//a b矩阵相乘,结果返回给a
void matrixMultiply(int a[][MATRIXSIZE], int la, int ca, int b[][MATRIXSIZE], int lb, int cb)
{
int i, j, k;
int c[MATRIXSIZE][MATRIXSIZE] = {0};
for(i = 0; i < la; i++)
for(j = 0; j < cb; j++)
for(k = 0; k < ca; k++)
c[i][j] += a[i][k] * b[k][j];
for(i = 0; i < la; i++)
for(j = 0; j < ca; j++)
a[i][j] = c[i][j];
}
//res中传入初始矩阵,base中传入基数矩阵,p为指数
void matrixFastPower(int base[][MATRIXSIZE], int p, int res[][MATRIXSIZE])
{
while(p)
{
if(p & 1) matrixMultiply(res, 1, MATRIXSIZE, base, MATRIXSIZE, MATRIXSIZE); //矩阵大小手动改一下
p >>= 1;
matrixMultiply(base, MATRIXSIZE, MATRIXSIZE, base, MATRIXSIZE, MATRIXSIZE);
}
}
例题
贴个题运用一下: 洛谷P1939, 点我看题
这个题的递推公式是隔一个相加,我当时愣是没想出来怎么构建这个矩阵,实在是太菜了。看完题解之后,给了自己一个大嘴巴子。
$$\begin{bmatrix}{F}_{n}&{F}_{n + 1}&{F}_{n + 2}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}&{F}_{3}\end{bmatrix} * {\begin{bmatrix}0&0&1\\1&0&0\\0&1&1\end{bmatrix}}^{n - 1}$$
看到这个矩阵递推公式就不用我多说什么了吧...
上AC代码:
//洛谷P1939
#include <stdio.h>
#include <string.h>
int tempa[3][3], tempb[3][3], M = 1000000007;
void matrixMultiply(int a[][3], int la, int ca, int b[][3], int lb, int cb)
{
int i, j, k;
for(i = 0; i < la; i++)
for(j = 0; j < ca; j++)
tempa[i][j] = a[i][j];
for(i = 0; i < lb; i++)
for(j = 0; j < cb; j++)
tempb[i][j] = b[i][j];
//清零
for(i = 0; i < la; i++)
for(j = 0; j < ca; j++)
a[i][j] = 0;
for(i = 0; i < la; i++)
for(j = 0; j < cb; j++)
for(k = 0; k < lb; k++){
a[i][j] += (tempa[i][k] * (long long)(tempb[k][j])) % M;
a[i][j] %= M;
}
}
void fastPower(int base[][3], int p, int res[][3])
{
while(p)
{
if(p & 1) matrixMultiply(res, 1, 3, base, 3, 3);
p >>= 1;
matrixMultiply(base, 3, 3, base, 3, 3);
}
}
int main(void)
{
int t, n;
scanf("%d", &t);
while(t--)
{
int base[3][3] = {
{0, 0, 1},
{1, 0, 0},
{0, 1, 1}
};
int res[3][3] = {1, 1, 1};
scanf("%d", &n);
if(n >= 4){
fastPower(base, n - 1, res);
printf("%d\n", res[0][0]);
}else
printf("1\n");
}
return 0;
}
本来还想再贴一题的,时间不允许。
大家一起加油!
Comments | NOTHING