快速幂与矩阵快速幂
快速幂
如果要计算 \(a\) 的 \(k\) 次，即 \(a^k\) 的结果，当然结果很一个很大的数，因此一般都会对结果模上 \(p\)。
如果用最朴素的写法是这样的：
long long res = 1; //res存结果
for(int i = 0; i < k; i ++)
{
res = res * a % p;
}
很显然这样的时间复杂度是 O(k) 的，完全达不到比赛的时限要求。
这时候快速幂就登场了！
-
快速幂的最基本的思路是：
-
预处理出 \(a^{2^0},a^{2^1},a^{2^2},...,a^{2^{\log k}}\) 总共 \(\log k\) 个数
-
将 \(a^k\) 用以上 \(logk\)个数来组合，那么就可以组合成
\[a^k = a^{2^{x_{1}}} \times a^{2^{x_{2}}} \times a^{2^{x_{3}}} \times ... \times a^{2^{x_{t}}} = a^{2^{x_{1}} + 2^{x_{2}} + 2^{x_{3}} + ... + 2^{x_{t}}} \] -
使得 \(2^{x_{1}} + 2^{x_{2}} + 2^{x_{3}} + ... + 2^{x_{t}} = k\)
很显然只要将 \(k\) 用二进制表示即可。 -
并且 \(a^{2^0} * a^{2^0} = a^{2^1}, a^{2^1} * a^{2^1} = a^{2^2}\) 以此类推...
-
这个时候我们就能发现在选择需要哪些数取决于 \(k\) 的二进制表示的当前位为 \(1\) 还是 \(0\)，因此我们可以一边判断 \(k\) 的每一位一边处理 \(a^{2^{x_{t}}}\)。
写成代码就是这个样子：
long long power(long long a, long long k, long long p)
{
long long res = 1; // 用于存结果
while(k)
{
if(k & 1) res = res * a % p; // k & 1用于判断k的最末位为1还是0，等价于k % 2 == 1
a = a * a % p;
K >>= 1; //等价于k /= 2
}
}
快速幂的时间复杂度为 \(O(logk)\)
矩阵快速幂
先来复习一下矩阵的运算
先定义两个矩阵 \(A = \begin{bmatrix} a_{1,1}& a_{1,2} \\ a_{2,1}& a_{2,2} \end{bmatrix}, B = \begin{bmatrix} b_{1,1}& b_{1,2} \\ b_{2,1}& b_{2,2} \end{bmatrix}\)
\(1.\) 矩阵加减法
\[A \pm B = \begin{bmatrix} a_{1,1}\pm b_{1,1} & a_{1,2}\pm b_{1,2} \\ a_{2,1}\pm b_{2,1} & a_{2,2}\pm b_{2,2} \end{bmatrix} \]- 注意加减法必须为同型矩阵
\(2.\)矩阵乘法
\[A \times B = \begin{bmatrix} a_{1,1}\times b_{1,1}+ a_{1,2}\times b_{2,1} & a_{1,1}\times b_{1,2}+ a_{1,2}\times b_{2,2} \\ a_{2,1}\times b_{1,1}+ a_{2,2}\times b_{2,1} & a_{2,1}\times b_{1,2}+ a_{2,2}\times b_{2,2} \end{bmatrix} \]- 注意，矩阵乘法要求满足 \(A\) 矩阵的行数等于 \(B\) 矩阵的列数，而矩阵快速幂是自己乘自己，因此矩阵快速幂的矩阵要求为方阵。
- 任意矩阵乘上单位矩阵(除主对角线上为1其余都为0的方阵)还是原矩阵
- 还是不懂矩阵乘法的戳这里
学会了矩阵乘法和快速幂，将其进行组合就有了矩阵快速幂。什么PPAP
矩阵快速幂就是将快速幂中的乘法变成矩阵乘法就可以了。其中 \(res\) 初始化为单位矩阵
先定义一个矩阵：
struct matrix
{
int n, m; // 表示方阵的大小
ll e[N][N];
}a;
再写一个矩阵乘法：
matrix mul(matrix a, matrix b)
{
matrix ans;
ans.n = a.n;
ans.m = b.m;
for(int i = 1; i <= ans.n; i ++)
for(int j = 1; j <= ans.m; j ++)
ans.e[i][j] = 0; // 初始化
for(int i = 1; i <= a.n; i ++)
for(int j = 1; j <= b.m; j ++)
for(int k = 1; k <= a.m; k ++)
ans.e[i][j] = (ans.e[i][j] + (a.e[i][k] * b.e[k][j]) % p) % p;
return ans;
}
最后写上我们的快速幂：
matrix power(matrix a, ll k)
{
matrix ans;
ans.n = a.n;
ans.m = a.m;
for(int i = 1; i <= ans.n; i ++)
{
for(int j = 1; j <= ans.m; j ++)
{
if(i == j)ans.e[i][j] = 1;
else ans.e[i][j] = 0;
}
}
while(k)
{
if(k & 1)ans = mul(ans, a);
a = mul(a, a);
k >>= 1;
}
return ans;
}
至此我们的矩阵快速幂就介绍完毕，我们还可以对他进行封装一下并用vector动态保存
struct matrix
{
int length, width;
vector<vector<ll>> m;
matrix(int l, int w): length(l), width(w), m(vector<vector<ll>>(l, vector<ll>(w))) {}
vector<ll>& operator[](const int pos) {
return m[pos];
}
friend matrix operator*(matrix &lhs, matrix &rhs) {
int l = lhs.length, w = rhs.width;
matrix res(l, w);
for (int i = 0; i < l; i++) {
for (int j = 0; j < w; j++) {
for (int k = 0; k < w; k++) {
res[i][j] += lhs[i][k] * rhs[k][j] % mod;
res[i][j] %= mod;
}
}
}
return res;
}
};
matrix qpow(matrix a, ll k) {
int n = a.length;
matrix res(n, n);
for (int i = 0; i < n; i++)
res[i][i] = 1;
while (k) {
if (k & 1) {
res = res * a;
}
a = a * a;
k >>= 1;
}
return res;
}
实际运用
矩阵快速幂可以应用到很多递推题上。
比如说很经典的斐波那契数列
这题因为 \(n < 2^{63}\) 这恐怖的数量级让一般的递推做法成为不可能。
因此我们可以可以用矩阵去递推。
首先我们构造有两个项的初始矩阵：
\[A = \begin{bmatrix} F_{n - 1}&F_{n - 2} \end{bmatrix} \]而我们需要得到的目标矩阵为：
\[B = \begin{bmatrix} F_{n}&F_{n - 1} \end{bmatrix} \]递推的式子为：
\[AC = B \]于是我们的目标就成为了求出这个 \(C\)
我们讲上面的式子展开：
结果应等于 \(B\),因此：
- \(F_{n - 1}\times c_{1,1} +F_{n - 2}\times c_{2,1} = F_{n - 1} + F_{n - 2}\)，所以 \(c_{1,1}和c_{2,1}\) 都等于 \(1\)。
- \(F_{n - 1}\times c_{2,1} +F_{n - 2}\times c_{2,2} = F_{n - 1}\)，所以 \(c_{1,2} = 1, c_{2,2} = 0\)。
于是我们就得到了 \(C\)
\[C=\begin{bmatrix} 1 &1 \\ 1&0 \end{bmatrix} \]因为每次乘上一个 \(C\) 就能让 \(A\) 矩阵的 \(a_{1,1}\) 向后递推一位，因此我们只要定义一个矩阵 \(A=\begin{bmatrix} F_{1} &F_{2} \end{bmatrix}\),再乘上 \(C^{n-2}\) 之后，其 \(a_{1,1}\) 即为 \(F_{n}\)。
最后我们附上代码：
#include<iostream>
#include<cstring>
#define ll long long
using namespace std;
const int N = 110, p = 1e9 + 7;
struct matrix
{
int n, m;
ll e[N][N];
}a, b;
matrix mul(matrix a, matrix b)
{
matrix ans;
ans.n = a.n;
ans.m = b.m;
for(int i = 1; i <= ans.n; i ++)
for(int j = 1; j <= ans.m; j ++)
ans.e[i][j] = 0; // 初始化
for(int i = 1; i <= a.n; i ++)
for(int j = 1; j <= b.m; j ++)
for(int k = 1; k <= a.m; k ++)
ans.e[i][j] = (ans.e[i][j] + (a.e[i][k] * b.e[k][j]) % p) % p;
return ans;
}
matrix power(matrix a, ll k)
{
matrix ans;
ans.n = a.n;
ans.m = a.m;
for(int i = 1; i <= ans.n; i ++)
{
for(int j = 1; j <= ans.m; j ++)
{
if(i == j)ans.e[i][j] = 1;
else ans.e[i][j] = 0;
}
}
while(k)
{
if(k & 1)ans = mul(ans, a);
a = mul(a, a);
k >>= 1;
}
return ans;
}
int main()
{
ll n;
cin >> n;
if(n <= 2)
{
puts("1");
return 0;
}
a.n = a.m = 2;
a.e[1][1] = a.e[1][2] = a.e[2][1] = 1;
a = power(a, n - 2);
b.n = 1, b.m = 2;
b.e[1][1] = b.e[1][2] = 1;
b = mul(b, a);
cout << b.e[1][1] << endl;
return 0;
}