矩阵乘法入门
矩阵
类似一个二维数组吧。
矩阵的运算
矩阵的加法
\[C_{i,j} = A_{i,j} + B_{i,j} \]我不知道有什么用。
矩阵的减法
\[C_{i,j} = A_{i,j} - B_{i,j} \]我也不知道有什么用。
矩阵的乘法
\[C_{i,j} = \sum_k^{m} A_{i,k}B_{k,j} \]即答案矩阵的第 \(i\) 行第 \(j\) 列数是 \(A\) 矩阵第 \(i\) 行第 \(k\) 个数和 \(B\) 矩阵第 \(k\) 行第 \(j\) 列的数相乘再相加的。显然,\(A\) 矩阵的列数和 \(B\) 矩阵的行数是一样的。注意:矩阵的乘法不满足交换律,但是满足结合律。矩阵的乘法可以用来加速递推。
矩阵乘法代码
struct matrix {
long long num[3][3];
void init() {
memset(num, 0, sizeof(num));
}
matrix operator * (const matrix& x) const {
matrix res;
res.init();
for (int i = 1; i <= 2; i ++)
for (int j = 1; j <= 2; j ++)
for (int k = 1; k <= 2; k ++)
res.num[i][j] += num[i][k] * x.num[k][j]
return res;
}
};
题目
1. 斐波那契数列
大家都知道,斐波那契数列是满足如下性质的一个数列:
\[F_n = \left\{\begin{aligned} 1 \space (n \le 2) \\ F_{n-1}+F_{n-2} \space (n\ge 3) \end{aligned}\right. \]请你求出 \(F_n \bmod 10^9 + 7\) 的值。
【数据范围】
对于 \(60\%\) 的数据,\(1\le n \le 92\);
对于 \(100\%\) 的数据,\(1\le n < 2^{63}\)。
思路
斐波那契数列的递推式很简单,为 \(f_n = f_{n - 1} + f_{n - 2}\)。
对于 \(60\%\) 的数据,一个一个算就可以了。
对于 \(100 \%\) 的数据,一个一个算直接爆炸了,我们就需要矩阵加速。
求出 \(f_n\) 需要 \(f_{n - 1}\) 和 \(f_{n - 2}\) ,所以我们构造一个矩阵 :
\[\begin{bmatrix} f_{n - 1} & f_{n - 2} \end{bmatrix} \]我们还需要构造一个矩阵 \(B\) 使 $\begin{bmatrix} f_{n - 1} & f_{n - 2} \end{bmatrix} \times B = \begin{bmatrix} f_{n} & f_{n - 1} \end{bmatrix} $
可以看出来,\(B\) 一定是一个 \(2 \times 2\) 的矩阵。为了更直观的表现,我们在矩阵旁边标上对应的数。
\[\ \ \ \ \ \ \ \ \ \ \ \begin{bmatrix} f_n & f_{n - 1} \end{bmatrix} \]\[\begin{bmatrix} f_{n - 1} \\ f_{n - 2}\end{bmatrix}\begin{bmatrix} 1\ \ \ \ \ \ \ & 1 \\ 1\ \ \ \ \ \ \ & 0 \end{bmatrix} \]这样就很好理解了,\(f_n\),是由 \(f_{n - 1}\) 和 \(f_{n - 2}\) 相加而得的,所以 \(f_n\) 对应的列, \(f_{n - 1}\) 和 \(f_{n - 2}\) 对应的行上的数都是 $ 1 $ ,\(f_{n - 1}\) 就是 \(f_{n - 1}\),所以 \(f_{n-1}\) 对应的列和 \(f_{n - 1}\) 对应的行上的数是 \(1\),由于 \(f_{n-1}\) 中不含 \(f_{n-2}\) ,所以 \(f_{n-1}\) 对应的列和 \(f_{n - 2}\) 对应的列上的数是 \(0\)。
求出了 \(B\) 矩阵后我们只要把矩阵 \(\begin{bmatrix} f_2 & f_1\end{bmatrix}\) 乘上 \(n-2\) 个 \(B\) 矩阵,就能得到 \(\begin{bmatrix} f_n & f_{n-1}\end{bmatrix}\)。
为什么是 \(n - 2\) 次不是 \(n\) 次呢?手算一下发现:
\[\begin{bmatrix} f_2 & f_1\end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} f_3 & f_2\end{bmatrix} \]\[\begin{bmatrix} f_3 & f_2\end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} f_4 & f_3\end{bmatrix} \]乘 \(1\) 次得到 \(f_3\) ,乘 \(2\) 次得到 \(f_4\),乘 \(n-2\) 次就得到了 \(f_n\)。
所以:
\[\begin{bmatrix} f_2 & f_1\end{bmatrix}(\begin{bmatrix} 1 & 1\end{bmatrix}) \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times ... \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} f_n & f_{n-1}\end{bmatrix} \]\[\begin{bmatrix} f_2 & f_1\end{bmatrix} \times \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^{n-2} = \begin{bmatrix} f_n & f_{n-1}\end{bmatrix} \]因为 \(1\le n < 2^{63}\),直接乘肯定要爆炸,需要快速幂。
矩阵快速幂和数的快速幂非常相似,唯一不同的地方就是 \(res\) 需要初始化为单位矩阵,就是对角线上的数都是 \(1\) 的矩阵。\(A \times 单位矩阵 = A\), 矩阵中的单位矩阵就类似于数中的 \(1\)。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1000000007;
struct matrix {
LL a[4][4];
matrix() {memset(a, 0, sizeof(a));}
matrix operator * (const matrix& x) const {
matrix res;
for (int i = 1; i <= 3; i ++)
for (int j = 1; j <= 3; j ++)
for (int k = 1; k <= 3; k ++) {
res.a[i][j] = (res.a[i][j] + (a[i][k] % mod * x.a[k][j] % mod) % mod) % mod;
}
return res;
}
matrix operator ^ (LL p) const {
matrix res, x;
for (int i = 1; i <= 3; i ++) res.a[i][i] = 1;
for (int i = 1; i <= 3; i ++)
for (int j = 1; j <= 3; j ++)
x.a[i][j] = a[i][j];
while (p) {
if (p & 1) res = res * x;
x = x * x;
p >>= 1;
}
return res;
}
};
matrix m, p;
int main() {
LL n;
scanf("%lld", &n);
if (n <= 2) {
puts("1");
return 0;
}
p.a[1][1] = p.a[1][2] = p.a[2][1] = 1;
m.a[1][1] = m.a[1][2] = 1;
m = m * (p ^ (n - 2));
printf("%lld\n", m.a[1][1]);
return 0;
}
2. Fibonacci 前 n 项和
大家都知道 \(Fibonacci\) 数列吧。
现在问题很简单,输入 \(n\) 和 \(m\),求 \(f_n\) 的前 \(n\) 项和 \(S_n\) \(\text{mod}\) \(m\)。
思路
这道题就比上道题升级了,需要求前 \(n\) 项和。
还是构造一个矩阵:
\[\begin{bmatrix} S_{n-2} & f_{n-1} & f_{n-2}\end{bmatrix} \]与上一题同理,我们需要构造出一个矩阵 \(B\) ,使 \(\begin{bmatrix} S_{n-2} & f_{n-1} & f_{n-2}\end{bmatrix} \times B = \begin{bmatrix} S_{n-1} & f_{n} & f_{n-1}\end{bmatrix}\)
为了更直观的表现,我们在矩阵旁边标上对应的数太麻烦了,不标了。可以很简单的构造出:
\(S_{n-1}\) 需要 \(S_{n-2} + f_{n - 1}\), 所以 \(S_{n-1}\) 列上的 \(S_{n-2}\) 行和 \(f_{n-1}\) 行上的数都是 \(1\)。
剩下的 \(f_{n-1}\) 和 \(f_{n-2}\) 都与上一题一样,就略去了。
遇上一题一样,还是快速幂,不过这一次是 \(n - 1\) 次,因为我们初始矩阵中放的是 \(S_{n - 2}\) 。
初始矩阵:
\[\begin{bmatrix} 1 & 1 & 1\end{bmatrix} \]代码
#include <bits/stdc++.h>
using namespace std;
long long n, m;
struct matrix {
long long num[4][4];
void init() {
memset(num, 0, sizeof(num));
}
matrix operator * (const matrix& x) const {
matrix res;
res.init();
for (int i = 1; i <= 3; i ++)
for (int j = 1; j <= 3; j ++)
for (int k = 1; k <= 3; k ++)
res.num[i][j] = (res.num[i][j] + (num[i][k] % m * x.num[k][j] % m) % m) % m;
return res;
}
matrix operator ^ (long long p) const {
matrix res, a;
res.init();
for (int i = 1; i <= 3; i ++)
res.num[i][i] = 1;
for (int i = 1; i <= 3; i ++)
for (int j = 1; j <= 3; j ++)
a.num[i][j] = num[i][j];
while (p) {
if (p & 1) res = res * a;
a = a * a;
p >>= 1;
}
return res;
}
};
matrix a, b;
int main() {
scanf("%lld%lld", &n, &m);
a.num[1][1] = a.num[1][2] = a.num[1][3] = 1;
b.num[1][1] = b.num[2][1] = b.num[2][2] = b.num[2][3] = b.num[3][2] = 1;
a = a * (b ^ (n - 1));
printf("%lld\n", a.num[1][1]);
return 0;
}
3. fib
\[f_1=f_2=0, f_n=7f_{n-1}+6f_{n-2}+5n+4×3^n,求f_n \ \text{mod} \ 10000 \]思路
这道题就更加升级了,增加了一些系数和常数,我们把这些全部放入矩阵。
\[\begin{bmatrix} f_{n-1} & f_{n-2} & 5(n-1) & 5 & 3^{n - 1} \end{bmatrix} \]这样我们的矩阵 \(B\) 就变成了一个巨型矩阵:
\[\begin{bmatrix} 7 & 1 & 0 & 0 & 0 \\ 6 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 12 & 0 & 0 & 0 & 3\end{bmatrix} \]$ f_n=7f_{n-1}+6f_{n-2}+5n+4×3^n$ ,所以 \(f_n\) 列 \(f_{n-1}\) 行就是 \(7\), \(f_{n - 2}\) 行就是 \(6\),\(5n = 5(n - 1) + 5\),所以 \(5(n-1)\) 行和 \(5\) 行都是 \(1\) ,\(4 \times 3^n = 12 \times 3^{n-1}\),所以 \(3^{n-1}\) 行就是 \(12\)。
\(f_{n-1} = f_{n-1}\) 与前几题同理,略。
\(5n = 5(n-1)+5\),所以 \(5n\) 列 \(5(n-1)\) 行和 \(5\) 行都是 \(1\) 。
\(5 = 5\),略。
\(3^n = 3^{n - 1} \times 3\),所以 \(3^n\) 列,\(3^{n-1}\) 行就是 \(3\)。
剩下的就是常规的快速幂,\(n-2\) 次方。
代码
#include <bits/stdc++.h>
using namespace std;
long long m = 10000;
struct matrix {
long long num[6][6];
void init() {
memset(num, 0, sizeof(num));
}
matrix operator * (const matrix& x) const {
matrix res;
res.init();
for (int i = 1; i <= 5; i ++)
for (int j = 1; j <= 5; j ++)
for (int k = 1; k <= 5; k ++)
res.num[i][j] = (res.num[i][j] + (num[i][k] % m * x.num[k][j] % m) % m) % m;
return res;
}
matrix operator ^ (long long p) const {
matrix res, a;
res.init();
for (int i = 1; i <= 5; i ++)
res.num[i][i] = 1;
for (int i = 1; i <= 5; i ++)
for (int j = 1; j <= 5; j ++)
a.num[i][j] = num[i][j];
while (p) {
if (p & 1) res = res * a;
a = a * a;
p >>= 1;
}
return res;
}
};
matrix a, b;
long long n;
int main() {
while (scanf("%lld", &n) != EOF) {
if (n == -1) break;
if (n <= 2) {
puts("0");
continue;
}
a.num[1][1] = 0,
a.num[1][2] = 0,
a.num[1][3] = 10,
a.num[1][4] = 5,
a.num[1][5] = 9,
b.num[1][1] = 7,
b.num[2][1] = 6,
b.num[3][1] = 1,
b.num[4][1] = 1,
b.num[5][1] = 12,
b.num[1][2] = 1,
b.num[3][3] = 1;
b.num[4][3] = 1,
b.num[4][4] = 1,
b.num[5][5] = 3;
a = a * (b ^ (n - 2));
printf("%lld\n", a.num[1][1]);
}
return 0;
}
总结
矩阵的 \(Latex\) 太难了,手都废了。
转眼5个月没写博客了,今天开始继续写。
标签:begin,end,matrix,矩阵,num,bmatrix,乘法 From: https://www.cnblogs.com/maniubi/p/17550177.html