線形漸化式（例：フィボナッチ数列）を単純に計算すると、nが10⁸を超えると時間がかかりすぎます。これを解決するのが「行列高速累乗」です。この手法は、二進展開と行列乗算を組み合わせ、計算量をO(log n)まで削減します。

基礎知識：二進高速累乗

整数のべき乗を高速に計算するアルゴリズムです。指数を二進数で分解し、繰り返し二乗しながら結果を構築します。

long long binpow(long long base, long long exp, long long mod) {
    long long result = 1;
    while (exp > 0) {
        if (exp & 1) result = (result * base) % mod;
        base = (base * base) % mod;
        exp >>= 1;
    }
    return result;
}

基礎知識：行列乗算

行列同士の積は、行と列の内積で定義されます。サイズが一致している必要があります。

using Matrix = vector<vector<long long>>;

Matrix matmul(const Matrix& A, const Matrix& B, long long mod) {
    int n = A.size(), m = B[0].size(), p = A[0].size();
    Matrix C(n, vector<long long>(m, 0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            for (int k = 0; k < p; ++k)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
    return C;
}

行列高速累乗の実装

単位行列から始めて、二進高速累乗の要領で行列を累乗します。

Matrix matpow(Matrix base, long long exp, long long mod) {
    int n = base.size();
    Matrix result(n, vector<long long>(n, 0));
    for (int i = 0; i < n; ++i) result[i][i] = 1;

    while (exp > 0) {
        if (exp & 1) result = matmul(result, base, mod);
        base = matmul(base, base, mod);
        exp >>= 1;
    }
    return result;
}

応用例1：フィボナッチ数列

状態ベクトル [Fₙ, Fₙ₋₁] を遷移行列 [[1,1],[1,0]] で更新します。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const ll MOD = 1e9+7;

Matrix matmul(const Matrix& A, const Matrix& B, ll mod) {
    int n = A.size(), m = B[0].size(), p = A[0].size();
    Matrix C(n, vector<ll>(m, 0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            for (int k = 0; k < p; ++k)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod;
    return C;
}

Matrix matpow(Matrix base, ll exp, ll mod) {
    int n = base.size();
    Matrix res(n, vector<ll>(n, 0));
    for (int i = 0; i < n; ++i) res[i][i] = 1;
    while (exp) {
        if (exp & 1) res = matmul(res, base, mod);
        base = matmul(base, base, mod);
        exp >>= 1;
    }
    return res;
}

void solve_fib(ll n) {
    if (n <= 2) { cout << 1 << endl; return; }
    Matrix T = {{1,1}, {1,0}};
    Matrix Fn = matpow(T, n-1, MOD);
    cout << Fn[0][0] << endl;
}

応用例2：一般化線形漸化

aₙ = p·aₙ₋₁ + q·aₙ₋₂ の形にも拡張可能です。

void solve_general(ll p, ll q, ll a1, ll a2, ll n, ll mod) {
    if (n == 1) { cout << a1 % mod << endl; return; }
    if (n == 2) { cout << a2 % mod << endl; return; }

    Matrix T = {{p % mod, q % mod}, {1, 0}};
    Matrix V = {{a2 % mod}, {a1 % mod}};
    T = matpow(T, n-2, mod);
    Matrix res = matmul(T, V, mod);
    cout << res[0][0] << endl;
}

応用例3：三項間漸化式

aₙ₊₁ = aₙ + aₙ₋₂ のような複雑なケースでも、3×3行列で表現できます。

void solve_tribo(ll n) {
    if (n <= 3) { cout << 1 << endl; return; }
    Matrix T = {{1,0,1}, {1,0,0}, {0,1,0}};
    Matrix Fn = matpow(T, n, MOD);
    cout << Fn[1][0] << endl;
}

応用例4：累積和の反復計算

下三角行列を使って、多次元前缀和を行列累乗で高速化できます。

void solve_prefix_sum(int n, int m, vector<ll>& arr) {
    Matrix A(n, vector<ll>(1));
    for (int i = 0; i < n; ++i) A[i][0] = arr[i];

    Matrix T(n, vector<ll>(n, 0));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j <= i; ++j)
            T[i][j] = 1;

    T = matpow(T, m, MOD);
    Matrix res = matmul(T, A, MOD);
    for (int i = 0; i < n; ++i)
        cout << res[i][0] << " ";
    cout << endl;
}

応用例5：線形合同法の高速シミュレーション

Xₙ₊₁ = (a·Xₙ + c) mod m も行列で表現可能。

void solve_lcg(ll m, ll a, ll c, ll x0, ll n, ll g) {
    Matrix V = {{x0 % m}, {1}};
    Matrix T = {{a, c % m}, {0, 1}};
    T = matpow(T, n, m);
    Matrix res = matmul(T, V, m);
    cout << (long long)(res[0][0] % g) << endl;
}