矩阵乘法
1. 矩阵乘法原理
原理
- 矩阵乘法是大学线性代数中的内容，其定义是：
A m × p × B p × n = C m × n A_{m \times p} \times B_{p \times n} = C _ {m \times n} Am×p×Bp×n=Cm×n
其中C中的每一项如下（A的第i行分别和B的第j列相乘后再相加）：
C
[
i
]
[
j
]
=
∑
k
=
1
p
A
[
i
]
[
k
]
×
B
[
k
]
[
j
]
C[i][j] = \sum _ {k = 1} ^ p A[i][k] \times B[k][j]
C[i][j]=k=1∑pA[i][k]×B[k][j]
-
矩阵乘法有结合律，这是可以使用快速幂算法的本质。（如果一种运算满足结合律，则这种运算就可以使用快速幂算法）
-
矩阵乘法不满足交换律。
-
这类题目一般需要我们先挖掘出矩阵乘法（这是最难的地方），然后使用快速幂求解。
2. AcWing上的矩阵乘法题目
AcWing 1303. 斐波那契前 n 项和
问题描述
分析
以下两个问题都使用矩阵乘法的方式解决，因为矩阵大小分别为 2 × 2 2 \times 2 2×2或者 3 × 3 3 \times 3 3×3的，两矩阵乘法计算次数可以忽略不计，则时间复杂度为 O ( l o g ( n ) ) O(log(n)) O(log(n))的。
求斐波那契数列
- 设 F n = [ f n , f n + 1 ] F_n = [f_n, f_{n+1}] Fn=[fn,fn+1]，则有如下递推公式：
[ f n , f n + 1 ] [ 0 1 1 1 ] = [ f n + 1 , f n + 2 ] [f_n, f_{n+1}] \left[ \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right] = [f_{n+1}, f_{n+2}] [fn,fn+1][0111]=[fn+1,fn+2]
- 因此，有如下递推公式：
F n = F 1 [ 0 1 1 1 ] n − 1 F 1 = [ 1 , 1 ] F_n = F_1 \left[ \begin{matrix} 0 & 1 \\ 1 & 1 \end{matrix} \right] ^ {n - 1} \quad \quad F_1 = [1, 1] Fn=F1[0111]n−1F1=[1,1]
- 假如设
A表示那个变换矩阵，则我们对 A n − 1 A^{n-1} An−1进行快速幂求解即可。
求斐波那契数列前n项和
- 设
F
n
=
[
f
n
,
f
n
+
1
,
S
n
]
F_n = [f_n, f_{n+1}, S_n]
Fn=[fn,fn+1,Sn]，其中
S
n
S_n
Sn表示数列前
n项和，则有如下递推公式：
[ f n , f n + 1 , S n ] [ 0 1 0 1 1 1 0 0 1 ] = [ f n + 1 , f n + 2 , S n + 1 ] [f_n, f_{n+1}, S_n] \left[ \begin{matrix} 0 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 0 & 1 \end{matrix} \right] = [f_{n+1}, f_{n+2}, S_{n+1}] [fn,fn+1,Sn]⎣⎡010110011⎦⎤=[fn+1,fn+2,Sn+1]
- 因此，有如下递推公式：
F n = F 1 [ 0 1 0 1 1 1 0 0 1 ] n − 1 F 1 = [ 1 , 1 , 1 ] F_n = F_1 \left[ \begin{matrix} 0 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 0 & 1 \end{matrix} \right] ^ {n - 1} \quad \quad F_1 = [1, 1, 1] Fn=F1⎣⎡010110011⎦⎤n−1F1=[1,1,1]
- 假如设
A表示那个变换矩阵，则我们对 A n − 1 A^{n-1} An−1进行快速幂求解即可。
代码
- C++
#include <iostream>
#include <cstring>
using namespace std;
typedef long long LL;
const int N = 3;
int n, m; // 求前n项和对m取模
// c[] = a[] * b[][]
void mul(int c[], int a[], int b[][N]) {
int temp[N] = {0};
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
temp[i] = (temp[i] + (LL)a[j] * b[j][i]) % m;
memcpy(c, temp, sizeof temp);
}
// c[][] = a[][] * b[][]
void mul(int c[][N], int a[][N], int b[][N]) {
int temp[N][N] = {0};
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
temp[i][j] = (temp[i][j] + (LL)a[i][k] * b[k][j]) % m;
memcpy(c, temp, sizeof temp);
}
int main() {
cin >> n >> m;
int f1[3] = {1, 1, 1};
int a[N][N] = {
{0, 1, 0},
{1, 1, 1},
{0, 0, 1}
};
n--;
while (n) {
if (n & 1) mul(f1, f1, a); // f1 = f1 * a
mul(a, a, a); // a = a * a
n >>= 1;
}
cout << f1[2] << endl;
return 0;
}
AcWing 1304. 佳佳的斐波那契
问题描述
-
问题链接：AcWing 1304. 佳佳的斐波那契
分析
- 题目已经给出
S(n)和T(n)的表达式，我们可以构造P(n)：
P ( n ) = n × S ( n ) − T ( n ) = ( n − 1 ) × f 1 + ( n − 2 ) × f 2 + . . . + f n − 1 P(n) = n \times S(n) - T(n) = (n - 1) \times f_1 + (n - 2) \times f_2 + ... + f_{n-1} P(n)=n×S(n)−T(n)=(n−1)×f1+(n−2)×f2+...+fn−1
则有：
P
(
n
+
1
)
−
P
(
n
)
=
S
(
n
)
P(n+1) - P(n) = S(n)
P(n+1)−P(n)=S(n)
因此如果我们可以求出P(n)和S(n)，则就可以求出T(n)，其表达是为：
T
(
n
)
=
n
×
S
(
n
)
−
P
(
n
)
T(n) = n \times S(n) - P(n)
T(n)=n×S(n)−P(n)。
- 总结一下，有如下递推式：
f n = f n − 1 + f n − 2 S n = S n − 1 + f n P n = P n − 1 + S n − 1 f_n = f_{n-1} + f_{n-2} \\ S_n = S_{n-1} + f_n \\ P_n = P_{n-1} + S_{n-1} fn=fn−1+fn−2Sn=Sn−1+fnPn=Pn−1+Sn−1
- 设 F n = [ f n , f n + 1 , S n , P n ] F_n = [f_n, f_{n+1}, S_n, P_n] Fn=[fn,fn+1,Sn,Pn]（这里的 F n F_n Fn和题目中的不同， f n f_n fn和题目中的 F n F_n Fn相同），则有：
[ f n , f n + 1 , S n , P n ] [ 0 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 ] = [ f n + 1 , f n + 2 , S n + 1 , P n + 1 ] [f_n, f_{n+1}, S_n, P_n] \left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{matrix} \right] = [f_{n+1}, f_{n+2}, S_{n+1}, P_{n+1}] [fn,fn+1,Sn,Pn]⎣⎢⎢⎡0100110001100011⎦⎥⎥⎤=[fn+1,fn+2,Sn+1,Pn+1]
- 因此有如下递推式：
F n = F 1 [ 0 1 0 0 1 1 1 0 0 0 1 1 0 0 0 1 ] n − 1 F 1 = [ 1 , 1 , 1 , 0 ] F_n = F_1 \left[ \begin{matrix} 0 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1 \end{matrix} \right] ^ {n - 1} \quad \quad F_1 = [1, 1, 1, 0] Fn=F1⎣⎢⎢⎡0100110001100011⎦⎥⎥⎤n−1F1=[1,1,1,0]
- 假如设
A表示那个变换矩阵，则我们对 A n − 1 A^{n-1} An−1进行快速幂求解即可。
代码
- C++
#include <iostream>
#include <cstring>
using namespace std;
typedef long long LL;
const int N = 4;
int n, m;
// c = a * b
void mul(int c[][N], int a[][N], int b[][N]) {
static int t[N][N];
memset(t, 0, sizeof t);
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = (t[i][j] + (LL)a[i][k] * b[k][j]) % m;
memcpy(c, t, sizeof t);
}
int main() {
cin >> n >> m;
int f1[N][N] = {1, 1, 1, 0};
int a[N][N] = {
{0, 1, 0, 0},
{1, 1, 1, 0},
{0, 0, 1, 1},
{0, 0, 0, 1},
};
int k = n - 1;
while (k) {
if (k & 1) mul(f1, f1, a); // f1 = f1 * a
mul(a, a, a); // a = a * a
k >>= 1;
}
cout << (((LL)n * f1[0][2] - f1[0][3]) % m + m) % m << endl;
return 0;
}
AcWing 1305. GT考试
问题描述
-
问题链接：AcWing 1305. GT考试
分析
-
这一题是AcWing 1052. 设计密码的一道扩展题目，分析方式仍然是动态规划。扩展方式是数据量，AcWing 1052. 设计密码中的
n值最大为50，这里的n最大可以取到 1 0 9 10 ^ 9 109。这是一种扩展方式，还有另外一种扩展方式，不扩展n，而是让不能包含多个字符串，对应题目是：AcWing 1053. 修复DNA，可以使用AC自动机解决。 -
本题的分析如下：
- 通过上面的分析，我们根据状态计算可以得到第
i层和第i+1层之间的关系，即
f ( i + 1 , 0 ) = a 0 , 0 × f ( i , 0 ) + a 1 , 0 × f ( i , 1 ) + . . . + a m − 1 , 0 × f ( i , m − 1 ) f ( i + 1 , 1 ) = a 0 , 1 × f ( i , 0 ) + a 1 , 1 × f ( i , 1 ) + . . . + a m − 1 , 1 × f ( i , m − 1 ) . . . f ( i + 1 , m − 1 ) = a 0 , m − 1 × f ( i , 0 ) + a 1 , m − 1 × f ( i , 1 ) + . . . + a m − 1 , m − 1 × f ( i , m − 1 ) f(i+1, 0) = a_{0,0} \times f(i, 0) + a_{1,0} \times f(i, 1) + ... + a_{m-1,0} \times f(i, m - 1) \\ f(i+1, 1) = a_{0,1} \times f(i, 0) + a_{1,1} \times f(i, 1) + ... + a_{m-1,1} \times f(i, m - 1) \\ ... \\ f(i+1, m-1) = a_{0,m-1} \times f(i, 0) + a_{1,m-1} \times f(i, 1) + ... + a_{m-1,m-1} \times f(i, m - 1) f(i+1,0)=a0,0×f(i,0)+a1,0×f(i,1)+...+am−1,0×f(i,m−1)f(i+1,1)=a0,1×f(i,0)+a1,1×f(i,1)+...+am−1,1×f(i,m−1)...f(i+1,m−1)=a0,m−1×f(i,0)+a1,m−1×f(i,1)+...+am−1,m−1×f(i,m−1)
如果我们令：
F
(
i
+
1
)
=
[
f
(
i
+
1
,
0
)
,
f
(
i
+
1
,
1
)
,
.
.
.
,
f
(
i
+
1
,
m
−
1
)
]
A
=
[
a
0
,
0
a
0
,
1
.
.
.
a
0
,
m
−
1
a
1
,
0
a
1
,
1
.
.
.
a
1
,
m
−
1
.
.
.
.
.
.
.
.
.
.
.
.
a
m
−
1
,
0
a
m
−
1
,
1
.
.
.
a
m
−
1
,
m
−
1
]
F(i+1) = [f(i+1, 0), f(i+1, 1), ..., f(i+1, m-1)] \\ A = \left[ \begin{matrix} a_{0,0} & a_{0,1} & ... & a_{0,m-1} \\ a_{1,0} & a_{1,1} & ... & a_{1,m-1} \\ ... & ... & ... & ... \\ a_{m-1,0} & a_{m-1,1} & ... & a_{m-1,m-1} \end{matrix} \right]
F(i+1)=[f(i+1,0),f(i+1,1),...,f(i+1,m−1)]A=⎣⎢⎢⎡a0,0a1,0...am−1,0a0,1a1,1...am−1,1............a0,m−1a1,m−1...am−1,m−1⎦⎥⎥⎤
则有：
F
(
i
+
1
)
=
F
(
i
)
×
A
F(i+1) = F(i) \times A
F(i+1)=F(i)×A
展开为：
[
f
(
i
+
1
,
0
)
,
f
(
i
+
1
,
1
)
,
.
.
.
,
f
(
i
+
1
,
m
−
1
)
]
=
[
f
(
i
,
0
)
,
f
(
i
,
1
)
,
.
.
.
,
f
(
i
,
m
−
1
)
]
×
[
a
0
,
0
a
0
,
1
.
.
.
a
0
,
m
−
1
a
1
,
0
a
1
,
1
.
.
.
a
1
,
m
−
1
.
.
.
.
.
.
.
.
.
.
.
.
a
m
−
1
,
0
a
m
−
1
,
1
.
.
.
a
m
−
1
,
m
−
1
]
[f(i+1, 0), f(i+1, 1), ..., f(i+1, m-1)] = \\ [f(i, 0), f(i, 1), ..., f(i, m-1)] \times \left[ \begin{matrix} a_{0,0} & a_{0,1} & ... & a_{0,m-1} \\ a_{1,0} & a_{1,1} & ... & a_{1,m-1} \\ ... & ... & ... & ... \\ a_{m-1,0} & a_{m-1,1} & ... & a_{m-1,m-1} \end{matrix} \right]
[f(i+1,0),f(i+1,1),...,f(i+1,m−1)]=[f(i,0),f(i,1),...,f(i,m−1)]×⎣⎢⎢⎡a0,0a1,0...am−1,0a0,1a1,1...am−1,1............a0,m−1a1,m−1...am−1,m−1⎦⎥⎥⎤
- 根据上面的分析可知，矩阵
A只与不合法串S有关，因此A矩阵是不变的。根据上面递推式可知：
F ( n ) = F ( 0 ) × A n F ( 0 ) = [ 1 , 0 , 0 , . . . ] F(n) = F(0) \times A ^{n} \quad \quad F(0) = [1, 0, 0, ...] F(n)=F(0)×AnF(0)=[1,0,0,...]
- 如何求解数组
A呢？如果从f(i, j)可以转移到f(i+1, k)，则让a[j, k]++。即让f(i+1, k) += f(i, j)：
f ( i + 1 , k ) = a 0 , k × f ( i , 0 ) + a 1 , k × f ( i , 1 ) + . . . + a j , k × f ( i , j ) + . . . + a m − 1 , k × f ( i , m − 1 ) f(i+1, k) = a_{0,k} \times f(i, 0) + a_{1,k} \times f(i, 1) +... + a_{j,k} \times f(i, j) + ... + a_{m-1,k} \times f(i, m - 1) f(i+1,k)=a0,k×f(i,0)+a1,k×f(i,1)+...+aj,k×f(i,j)+...+am−1,k×f(i,m−1)
-
求出向量
F(n)后，最后的答案就是向量F(n)中所有的元素之和。 -
这是一类问题，凡是动态规划中两层之间的转移形式是乘以一个固定矩阵的，都可以使用快速幂优化。
代码
- C++
#include <iostream>
#include <cstring>
using namespace std;
const int N = 25;
int n, m, mod; // 准考证号为 n 位数, 不吉利数字为m位
char str[N]; // 不吉利数字串
int ne[N]; // KMP求str自身的ne
int a[N][N]; // 转移矩阵
void mul(int c[][N], int a[][N], int b[][N]) {
static int t[N][N];
memset(t, 0, sizeof t);
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
for (int k = 0; k < N; k++)
t[i][j] = (t[i][j] + a[i][k] * b[k][j]) % mod;
memcpy(c, t, sizeof t);
}
int qmi(int k) {
int f0[N][N] = {1};
while (k) {
if (k & 1) mul(f0, f0, a); // f0 = f0 * a
mul(a, a, a); // a = a * a;
k >>= 1;
}
int res = 0;
for (int i = 0; i < m; i++) res = (res + f0[0][i]) % mod;
return res;
}
int main() {
cin >> n >> m >> mod;
cin >> str + 1;
// KMP
for (int i = 2, j = 0; i <= m; i++) {
while (j && str[i] != str[j + 1]) j = ne[j];
if (str[i] == str[j + 1]) j++;
ne[i] = j;
}
// 初始化A[i][j]
for (int j = 0; j < m; j++)
for (int c = '0'; c <= '9'; c++) {
int k = j; // 原字符串后缀和str前缀匹配的长度
while (k && str[k + 1] != c) k = ne[k];
if (str[k + 1] == c) k++;
if (k < m) a[j][k]++;
}
// F[n] = F[0] * A^n
cout << qmi(n) << endl;
return 0;
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
