¶ 矩阵快速幂加速递推
可以简单把矩阵看作一个二维数组。
一个 行 列的矩阵 一共包含 个元素。第 行，第 列的元素为 。
¶ 矩阵乘法
矩阵乘法的定义为：
即一个 行 列的矩阵 ，与一个 行 列的矩阵 相乘，可以得到一个 行 列的矩阵 。
其中， 的第 行第 列的值为：
怎么理解这个式子呢？我们可以把 与 看作是两张价目表：
- 表示从城市 到城市 乘坐火车的价格
- 表示从城市 到城市 乘坐飞机的价格。
如上图所述， 表示的就是先乘坐火车从 到达一个中转站，然后由中转站乘坐飞机到达 。显然中转站只能是 或 （允许火车坐回到原地）。 的值就是把“每个方案的两趟票价的乘积”相加。
实际上， 就是把 “ 的第 行”与 “ 的第 列”依次相乘，然后把所有乘积累加。
为什么是“所有乘积相加”？可以考虑：
- 表示从城市 到城市 乘坐火车有几种方案。
- 表示从城市 到城市 乘坐飞机有几种方案。
那么乘积中的 就表示从城市 到城市 ，先火车再飞机有几种方案。根据数学中的加法原理和乘法原理，自然就是把“所有乘积相加”了。实际上这也是矩阵在图论中的一个常见应用，可以很方便计算在图中走 步、每个点到其他点有多少种走法。
¶ 矩阵乘法的性质
- 乘法的两个矩阵必须满足第一个矩阵的列数等于第二个矩阵的行数（中转站数量一致）
- 不满足交换律，即不满足
- 满足结合律，即满足 =
- 行 列的单位矩阵 ，为除了对角线()上为 ，其他位置都是 的矩阵。
- .
- 因为满足结合律，所以我们可以使用快速幂来快速计算矩阵的幂。
¶ 矩阵乘法完成一次递推
我们经常需要做类似于：
- .
- .
这样的递推式。
我们可以用一个 行 列的矩阵 作为初始值矩阵， 行 列的矩阵 作为递推矩阵，并希望能计算出初始矩阵递推一次后的矩阵 。我们可以采用这样的矩阵乘法的方式，完成一次递推。
¶ 例子 1
比如对于递推式 来说，因为由两项就可以推出来下一项。
- 我们可以构建初始矩阵 为 .
- 希望递推一次后得到的 为 .
那么递推矩阵是一个两行两列的，根据矩阵乘法的性质，我们可以很容易得出递推矩阵：
- 我们希望 为 ，而 ，因此递推矩阵的第一列从上往下就分别为
- 我们希望 为 ，而 ，因此递推矩阵的第一列从上往下就分别为
只需要抓住矩阵乘法的特点：“乘法结果的第 行、第 列的值，就是初始矩阵第 行分别与递推矩阵的第 列相乘后相加”。
用来优化递推时，初始矩阵一般都是一行的，所以可以进一步简单记住特点：“乘法结果的第 个数，就是初始矩阵依次与递推矩阵的第 列相乘后相加”。
¶ 例子 2
如果递推式为：
这里首先 与前三项都有关，并且有一个常数项，因此我们可以如下构建矩阵：
- 初始矩阵为 第四项的 用来达成最后的
- 希望一次乘法后的矩阵为 .
- 一项项看结果矩阵的每个数
- ，所以递推矩阵的第一列为
- ，所以递推矩阵的第二列为
- ，所以递推矩阵的第三列为
- ，所以递推矩阵的第四列为
¶ 矩阵乘法加速递推
既然 可以得到一次递推后的结果，那么 就可以得到递推两次后的结果。
即递推 次后的结果。而 可以使用矩阵快速幂，在 快速算完。
¶ 模板
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int MAXR = 3; // 矩阵最大行数
const int MAXC = 3; // 矩阵最大列数
int MOD;
struct Mat
{
int r, c; // 行数列数
int m[MAXR + 1][MAXC + 1];
Mat() { memset(m, 0, sizeof(m)); }
};
Mat operator*(const Mat &a, const Mat &b)
{
assert(a.c == b.r);
Mat res;
res.r = a.r, res.c = b.c;
for (int i = 1; i <= a.r; i++)
for (int k = 1; k <= a.c; k++)
{
int temp = a.m[i][k];
for (int j = 1; j <= b.c; j++)
res.m[i][j] = (res.m[i][j] + temp * b.m[k][j]) % MOD;
}
return res;
}
Mat quick_pow(Mat a, int b)
{
assert(a.r == a.c);
Mat res;
res.r = res.c = a.r;
for (int i = 1; i <= res.r; i++)
res.m[i][i] = 1;
while (b > 0)
{
if (b % 2)
res = res * a;
a = a * a;
b = b / 2;
}
return res;
}
signed main()
{
ios::sync_with_stdio(false);
cin.tie(0);
int x;
cin >> x >> MOD;
// 构造矩阵
Mat now;
now.r = 1, now.c = 3;
now.m[1][1] = 1, now.m[1][2] = a + 1, now.m[1][3] = 1;
Mat ans;
ans.r = ans.c = 3;
ans.m[1][1] = 0, ans.m[1][2] = 0, ans.m[1][3] = 0;
ans.m[2][1] = 1, ans.m[2][2] = a, ans.m[2][3] = 0;
ans.m[3][1] = 0, ans.m[3][2] = 1, ans.m[3][3] = 1;
// 计算答案
ans = quick_pow(ans, x - 1);
now = now * ans;
cout << now.m[1][1] % MOD << "\n";
return 0;
}
¶ 其他
- 如果不会重载运算符，可以使用
Mat mul(const Mat &a, const Mat &b){...}的函数封装乘法。 - 如果矩阵非常大（行数列数很多），矩阵乘法的三重循环需要注意顺序，尽可能让访问数组时为同一行的连续访问，可以加快速度。如果只是非常小的矩阵的话，就无所谓了。
- 以斐波那契数列为例，因为 ，矩阵乘法最终就是在找这个 。显然算出这个 后，就有 。
- 矩阵快速幂本质上可以说是在用更好的结构（矩阵描述），来倍增地处理递推。
¶ 高斯消元模板
¶ 只需要有解时给出解
//------------高斯消元模板-Start--------------
const double eps = 1e-12;
double a[105][105]; // 高斯消元的方程
//如果无解/多解：返回 false
//有解：返回 true
bool gauss(int x)
{
for (int i = 1; i <= x; i++)
{
// 当前方程是第 i 个方程
// 试图把当前方程第 i 个变量系数调整为 1
// 找到下面第 i 个变量系数最大的方程
int row = i;
for (int j = i; j <= x; j++)
if (fabs(a[row][i]) < fabs(a[j][i]))
row = j;
// 用系数最大的那个作为当前方程
if (row != i)
std::swap(a[row], a[i]);
// 如果最大的系数都是 0 了，下面所有方程的就都是 0 了
double div1 = a[i][i];
if (fabs(div1) < eps)
return false;
// 调整当前方程
for (int j = 1; j <= x + 1; j++)
a[i][j] /= div1;
// 用当前方程调整下面的方程，把下面方程的第 i 个系数小调
for (int j = i + 1; j <= x; j++)
{
double div2 = a[j][i];
for (int k = 1; k <= x + 1; ++k)
a[j][k] -= a[i][k] * div2;
}
}
// 每个方程向上把对应系数清零
for (int i = x; i >= 1; i--)
for (int j = i - 1; j >= 1; j--)
{
a[j][x + 1] -= a[j][i] * a[i][x + 1];
a[j][i] = 0;
}
return true;
}
//------------高斯消元模板-End--------------
¶ 需要判断到底是无解还是多解
//------------高斯消元模板-Start--------------
const double eps = 1e-12;
double a[105][105]; // 高斯消元的方程
// 无解：返回 -1
// 多解：返回 0
// 有解：返回 1
int gauss(int x)
{
bool flag = false; // 多解flag
int line = 1; // 当前方程的编号
for (int i = 1; i <= x; i++) // 主元编号
{
// 当前方程是第 line 个方程
// 试图把当前方程主元变量(第i个变量)系数调整为 1
// 找到下面主元系数最大的方程
int row = line;
for (int j = line + 1; j <= x; j++)
if (fabs(a[row][i]) < fabs(a[j][i]))
row = j;
// 用系数最大的那个作为当前方程
if (row != line)
std::swap(a[row], a[line]);
// 如果最大的系数都是 0 了，下面所有方程的就都是 0 了，不用往下调整了
double div1 = a[line][i];
if (fabs(div1) < eps)
continue;
// 调整当前方程
for (int j = 1; j <= x + 1; j++)
a[line][j] /= div1;
// 用当前方程调整下面的方程，把下面方程的第 i 个系数小调
for (int j = line + 1; j <= x; j++)
{
double div2 = a[j][i];
for (int k = 1; k <= x + 1; ++k)
a[j][k] -= a[line][k] * div2;
}
line++;
}
// 判无解和多解
if (line <= x)
{
for (int i = line; i <= x; i++)
if (a[i][x + 1] != 0)
return -1;
return 0;
}
// 有唯一解时依次处理
for (int i = x; i >= 1; i--)
{
for (int j = i - 1; j >= 1; j--)
{
a[j][x + 1] -= a[j][i] * a[i][x + 1];
a[j][i] = 0;
}
}
return 1;
}
//------------高斯消元模板-End--------------
¶ 线性基
https://oi.men.ci/linear-basis-notes/
const int MAXL = 60;
struct LinearBasis
{
long long a[MAXL + 1];
LinearBasis()
{
std::fill(a, a + MAXL + 1, 0);
}
void insert(long long t)
{
// 逆序枚举二进制位
for (int j = MAXL; j >= 0; j--)
{
// 如果 t 的第 j 位为 0，则跳过
if (!(t & (1ll << j))) continue;
// 如果 a[j] != 0，则用 a[j] 消去 t 的第 j 位上的 1
if (a[j]) t ^= a[j];
else
{
// 找到可以插入 a[j] 的位置
// 用 a[0...j - 1] 消去 t 的第 [0, j) 位上的 1
// 如果某一个 a[k] = 0 也无须担心，因为这时候第 k 位不存在于线性基中，不需要保证 t 的第 k 位为 0
for (int k = 0; k < j; k++) if (t & (1ll << k)) t ^= a[k];
// 用 t 消去 a[j + 1...L] 的第 j 位上的 1
for (int k = j + 1; k <= MAXL; k++) if (a[k] & (1ll << j)) a[k] ^= t;
// 插入到 a[j] 的位置上
a[j] = t;
// 不要忘记，结束插入过程
return;
}
// 此时 t 的第 j 位为 0，继续寻找其最高位上的 1
}
// 如果没有插入到任何一个位置上，则表明 t 可以由 a 中若干个元素的异或和表示出，即 t 在 span(a) 中
}
};