矩阵快速幂以及其优化【华东交大课程】

快速幂基础：C++快速幂_Kicamon的博客-CSDN博客

矩阵快速幂就是在快速幂的基础上结合矩阵运算的用法，其用途较为广泛，可以很大程度上优化代码。

一、矩阵快速幂的实现

在矩阵快速幂中，我们利用类来存储矩阵及其运算法则。

我们先来看看代码的具体实现，再来讲解其语法。

struct Matrix
{ll a[N][N];Matrix()//构造函数{memset(a,0,sizeof a);}void init()//单元矩阵{for(int i = 0;i <= n;i++)a[i][i] = 1;}Matrix operator*(const Matrix& B)const//运算符*的重构{Matrix C;for(int i = 0;i < n;i++)for(int k = 0;k < n;k++)for(int j = 0;j < n;j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;return C;}Matrix operator^(const int& t)const{Matrix A = (*this),res;res.init();int p = t;while(p){if(p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};

该结构体使用到的陌生语法有：构造函数、运算符重载和this指针，下边我们来看看着两个语法：
- 构造函数：
  
  构造函数的特点是在结构体或者类中，函数名与结构体名称相同且其声明前方没有返回值。它会在结构体创建的时候自动调用，完成对结构体的初始化。与之对应的是析构函数，作用是对结构体的清除，感兴趣的同学可以自行了解。
- 运算符重载：
  - 运算发重载的含义：指的是将已有的运算符赋予新的功能，并通过识别使用是的数据类型和数目来判断采用哪种操作。c/c++中有许多运算符其实已经经过了重构，例如*，在将其作用于地址的时候，它的作用是解引用；将其用于两个数字的时候，它的作用是得到两数的乘积。运算符重载的存在，使得代码可以变得更加简洁易懂。
  - 在这里，我们将 * 运算符重载出求两矩阵的积的作用，将 ^ 运算符重载出求矩阵的次方的作用。当然，在 ^ 中我们使用了重载之后的 * 运算符。
- this指针
  - this指针存在于成员函数中，是一个用来指向对象本身地址的指针。例如上面代码中，this指针将该指针所在的结构体本身赋值给了结构体A。

二、例题讲解

1、利用矩阵快速幂推导斐波那契数列

在这之前，我们推导斐波那契数列的方法是这样的：

f[1] = f[2] = 1;
for(int i = 3;i <= n;i++)f[i] = f[i - 1] + f[i - 2];

但是，学了矩阵快速幂之后，我们可以利用另一种方法来推导它：

所以我们可以利用矩阵快速幂先得到状态转换矩阵的n次方，再将其于初始矩阵相乘

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;struct Matrix
{ll a[2][2];Matrix(){memset(a, 0, sizeof a);}void init(){for (int i = 0; i < 2; i++)a[i][i] = 1;}Matrix operator*(const Matrix& B)const{Matrix C;for (int i = 0; i < 2; i++)for (int k = 0; k < 2; k++)for (int j = 0; j < 2; j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;return C;}Matrix operator^(const int& t)const{Matrix A = (*this), res;res.init();int p = t;while (p){if (p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};Matrix bace;int main()
{bace.a[0][0] = 0, bace.a[0][1] = 1;bace.a[1][0] = 1, bace.a[1][1] = 1;int n;cin >> n;if (n == 1 || n == 2)puts("1");else{Matrix ans = bace ^ n;cout << ans.a[0][1] << endl;//f[n] = ans.a[0][0] * f[0] + ans.a[0][1] * f[1],f[0] = 0,f[1] = 1}return 0;
}

2、求S_n

我们先来看看式子的推导过程：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 2, MOD = 1e9 + 7;struct Matrix
{ll a[2][2];Matrix(){memset(a, 0, sizeof a);}void init(){for (int i = 0; i < 2; i++)a[i][i] = 1;}Matrix operator*(const Matrix& B)const{Matrix C;for (int i = 0; i < 2; i++)for (int k = 0; k < 2; k++)for (int j = 0; j < 2; j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;return C;}Matrix operator^(const int& t)const{Matrix A = (*this), res;res.init();int p = t;while (p){if (p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};Matrix bace;int main()
{int a, b, n, m;cin >> a >> b >> n >> m;bace.a[0][0] = a, bace.a[0][1] = b;bace.a[1][0] = 1, bace.a[1][1] = a;Matrix t = bace ^ n;int Xn = t.a[0][0];cout << (Xn * 2) % MOD << endl;return 0;
}

3、梦想之弧

先看看推导过程：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int MOD = 1e9 + 7;int A0, AX, AY, B0, BX, BY;struct Matrix
{int a[6][6];Matrix(){memset(a, 0, sizeof a);}void init(){for (int i = 0; i < 5; i++)a[i][i] = 1;}Matrix operator*(const Matrix& B)const{Matrix C;for (int i = 0; i < 5; i++)for (int k = 0; k < 5; k++)for (int j = 0; j < 5; j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;return C;}Matrix operator^(const int& t)const{Matrix A = (*this), res;res.init();int p = t;while (p){if (p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};Matrix bace, once;int main()
{int n;cin >> n;cin >> A0 >> AX >> AY;cin >> B0 >> BX >> BY;bace.a[0][0] = AX * BX, bace.a[0][4] = 1;bace.a[1][0] = AX * BY, bace.a[1][1] = AX;bace.a[2][0] = AY * BX, bace.a[2][2] = BX;bace.a[3][0] = AY * BY, bace.a[3][1] = AY, bace.a[3][2] = BY, bace.a[3][3] = 1;bace.a[4][4] = 1;int a1 = A0 * AX + AY, b1 = B0 * BX + BY;once.a[0][0] = a1 * b1, once.a[0][1] = a1, once.a[0][2] = b1, once.a[0][3] = 1, once.a[0][4] = A0 * B0;Matrix ans = once * (bace ^ (n - 1));cout << ans.a[0][4] << endl;return 0;
}

三、矩阵乘法的优化

在上边我们用到的矩阵乘法是三重循环的，时间复杂度为O(n³)，在进行较大的矩阵运算时便会显得吃力。下边我们介绍两种方法来对矩阵的乘法进行优化。

1、省去为0的运算

Matrix operator*(const Matrix& B)const
{Matrix C;for(int i = 0;i < n;i++){for(int k = 0;k < n;k++){if(!a[i][k])continue;//为0时乘积为0，直接跳过for(int j = 0;j < n;j++){C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;}}}
}

2、利用二进制优化

这个方法需要用到模板类bitset，我们先来学习一下待会要用到的语法。

bitset是将数据用二进制储存下来并且进行编辑。在定义的时候，我们需要规定其长度，例如bitset<10>a，就是创建一个十位的二进制空间a。

bitset<10>a;//创建，a的初始状态为：0000000000a.set(3);//将第3位变成1，此时a的状态为：0000000100a.reset();//将所有位清零，此时a的状态为：0000000000a.set(4),a.set(5),a.set(2);a.count();//统计a中的1的数量，此时a的状态为：0000011010，a.count() = 3

在掌握以上语法之后，我们就可以来看看它的优化方法了：

首先我们创建两个数组：

bitset<1000>bt[805][3],ct[805][3];

然后，我们来看看它的运算过程：

如此如此便可得到结果，我们来看看最终的代码，这时候的时间复杂度为O(n²)：

#include<iostream>
#include<bitset>
using namespace std;
const int N = 805;
bitset<1000>bt[N][3], ct[N][3];int main()
{int n, c;cin >> n;for (int i = 1; i <= n; i++)for (int j = 0; j < 3; j++)bt[i][j].reset(), ct[i][j].reset();for(int i = 1;i <= n;i++)for (int j = 1; j <= n; j++){cin >> c;bt[i][c % 3].set(j);}for(int i = 1;i <= n;i++)for (int j = 1; j <= n; i++){cin >> c;ct[j][c % 3].set(i);}for (int i = 1; i <= n; i++){for (int j = 1; j <= n; j++){int c = ((bt[i][1] & ct[j][1]).count()+ 2 * (bt[i][1] & ct[j][2]).count()+ 2 * (bt[i][2] & ct[j][1]).count()+ 2 * (bt[i][2] & ct[j][2]).count()) % 3;cout << c << ' ';}cout << endl;}return 0;}

四、矩阵快速幂加速DP

1、Clarke and digits

这道题可是重量级。它运用的是矩阵快速幂和数位DP（数位DP将在之后吉吉会详细地讲解）。

首先，我们要考虑的是如何表示状态，根据题干，我们可以提取三个信息：长度、模数、数位。

我们创建一个三位数组，每一维度分别表示一条信息。如==a[i] [j] [k]==表示的长度为i、模数为j且尾数为k的数的数量。

下面来看看它的表示原理和状态更新的方式：

但是矩阵应该由二维数组表示，那我们开看看，能不能通过某种方式来减去一个维度。答案当然是可以的，先来看看压缩之后的结果：

a[i] [j] [k] = a[i] [j * 10 + k]

原理何在？我们知道，j是对7取模之后的结果，所以它的范围为[0,7]；k是尾数，尾数即最后一位数，所以它的范围为[0,9]。所以数组a的两维可以表示为一个7 * 10的矩阵，而上图的压缩方法是将其变成了一个1 * 70的矩阵，0 ~ 9表示原来的第一行的状态，10 ~ 19表示原来的第二行的状态，以此类推。那么，我们就得到一个二维数组，即我们需要的矩阵了。

学到这里，结果就明了了，构建完状态矩阵和初始矩阵之后，其他的套模板即可：

状态矩阵：

void init(int k)
{memset(bace.a,0,sizeof bace.a);for(int i = 0;i < 70;i++){int x = i / 10,y = i % 10;//余数和尾数for(int j = 0;j < 70;j++){int a = j / 10,b = j % 10;//上一状态的余数和尾数if((a * 10 + y) == x && y + b != k)bace.a[j][i] = 1;}}
}

初始矩阵：

for(int i = 1;i < 10;i++)once.a[0][(i % 7) * 10 + i] = 1;

完整代码如下：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 71, MOD = 1e9 + 7;struct Matrix
{int a[N][N];Matrix(){memset(a, 0, sizeof a);}void init(){for (int i = 0; i < N; i++)for (int j = 0; j < N; j++)a[i][j] = (i == j);}Matrix operator*(const Matrix& B)const{Matrix C;for (int i = 0; i < N; i++)for (int k = 0; k < N; k++)for (int j = 0; j < N; j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % MOD;return C;}Matrix operator^(const int& t)const{Matrix A = (*this), res;res.init();int p = t;while (p){if (p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};Matrix bace, once;void init(int k)
{memset(bace.a, 0, sizeof bace.a);for (int i = 0; i < 70; i++){int x = i / 10, y = i % 10;for (int j = 0; j < 70; j++){int a = j / 10, b = j % 10;if ((a * 10 + y) % 7 == x && y + b != k)bace.a[j][i] = 1;}}for (int i = 0; i < 10; i++)bace.a[i][70] = 1;bace.a[70][70] = 1;
}int main()
{memset(once.a, 0, sizeof once.a);for (int i = 1; i < 10; i++)once.a[0][(i % 7) * 10 + i] = 1;int t;cin >> t;while (t--){int l, r, k;cin >> l >> r >> k;init(k);Matrix ans1 = once * (bace ^ (l - 1));Matrix ans2 = once * (bace ^ r);ll ans = ans2.a[0][70] - ans1.a[0][70];ans = (ans % MOD + MOD) % MOD;cout << ans << endl;}return 0;
}

2、我也不知道这题叫啥，因为找不到

这里没有什么变编程上的语法知识，我们直接看看推导过程即可：

以下是费马小定理的衍生：

我们在求K_n时便可对p-1取模，防止数据过大

代码如下：

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N = 3;
int n, a, b, c, p;struct Matrix
{int a[N][N];Matrix(){memset(a, 0, sizeof a);}void init(){for (int i = 0; i < N; i++)a[i][i] = 1;}Matrix operator*(const Matrix& B)const{Matrix C;for (int i = 0; i < N; i++)for (int k = 0; k < N; k++)for (int j = 0; j < N; j++)C.a[i][j] = (C.a[i][j] + a[i][k] * B.a[k][j]) % p;return C;}Matrix operator^(const int& t)const{Matrix A = (*this), res;int p = t;res.init();while (p){if (p & 1)res = res * A;A = A * A;p >>= 1;}return res;}
};ll pow(int a, int b)
{int ans = 1;while (b){if (b & 1)ans = (ans * a) % p;a = (a * a) % p;b >>= 1;}return 1;
}Matrix bace, once;void init()
{memset(bace.a, 0, sizeof bace.a);memset(once.a, 0, sizeof once.a);bace.a[0][0] = c, bace.a[0][1] = 1, bace.a[0][2] = 1;bace.a[1][0] = 1, bace.a[2][2] = 1;once.a[0][0] = b, once.a[1][0] = 0, once.a[2][0] = b;
}int main()
{cin >> n >> a >> b >> c >> p;if (n == 1)puts("1");else if (n == 2)cout << pow(a, b) << endl;else{p--;Matrix ans = bace ^ (n - 2);p++;ans = once * ans;cout << ans.a[0][0] << endl;}return 0;
}