这篇 Blog 记录了一些线性代数的板子,以便不时之需。

矩阵快速幂

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
const int maxN = 105;
struct matrix{
    long long a[maxN][maxN];
}mat, m1;
long long n, p;
matrix operator*(matrix x, matrix y) {
    matrix z;
    for (long long i = 1; i <= n; i++)
        for (long long j = 1; j <= n; j++) {
            z.a[i][j] = 0;
            for (long long k = 1; k <= n; k++)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j]) % MOD;
        }
    return z;
}
matrix ident(long long k) {
    matrix id;
    for (long long i = 1; i <= k; i++)
        for (long long j = 1; j <= k; j++)
            if (i == j) id.a[i][j] = 1;
            else id.a[i][j] = 0;
    return id;
}
matrix qpow(matrix s, long long k) {
    matrix res = ident(n);
    for ( ; k; k >>= 1, s = s * s)
        if (k & 1) res = res * s;
    return res;
}
int main() {
    cin >> n >> p;
    for (long long i = 1; i <= n; i++)
        for (long long j = 1; j <= n; j++)
            cin >> mat.a[i][j];
    matrix ans = qpow(mat, p);
    for (long long i = 1; i <= n; i++) {
        for (long long j = 1; j <= n; j++)
            cout << ans.a[i][j] << ' ';
        cout << endl;
    }
    return 0;
}

矩阵加速的斐波那契

$$\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix}$$

#include <bits/stdc++.h>
using namespace std;
const int MOD = 1e9 + 7;
const int maxN = 5;
const int N = 2;
struct matrix{
    long long a[maxN][maxN];
}mat;
long long p;
matrix operator*(matrix x, matrix y) {
    matrix z;
    for (long long i = 1; i <= N; i++)
        for (long long j = 1; j <= N; j++) {
            z.a[i][j] = 0;
            for (long long k = 1; k <= N; k++)
                z.a[i][j] = (z.a[i][j] + x.a[i][k] * y.a[k][j]) % MOD;
        }
    return z;
}
matrix ident() {
    matrix id;
    id.a[1][1] = id.a[2][2] = 1;
    id.a[2][1] = id.a[1][2] = 0;
    return id;
}
matrix qpow(matrix s, long long k) {
    matrix res = ident();
    for ( ; k; k >>= 1, s = s * s)
        if (k & 1) res = res * s;
    return res;
}
int main() {
    cin >> p;
    mat.a[1][1] = mat.a[1][2] = mat.a[2][1] = 1;
    mat.a[2][2] = 0;
    matrix ans = qpow(mat, p);
    cout << ans.a[2][1] << endl;
    return 0;
}