行列の構造と基本演算
行列は \(m\) 行 \(n\) 列の数値配置として定義され、各成分を \(a_{i,j}\) で表すとき、行列 \(\mathbf{A}\) は以下の形式で記述される。
\[ \mathbf{A} = \begin{pmatrix} a_{1,1} & \cdots & a_{1,n} \\ \vdots & \ddots & \vdots \\ a_{m,1} & \cdots & a_{m,n} \end{pmatrix} \]特殊な行列として「単位行列」が存在する。これは正方行列であり、主対角成分がすべて \(1\)、非対角成分が \(0\) の構造を持つ。スカラー算における \(1\) と等価であり、任意の行列 \(\mathbf{A}\) に対して \(\mathbf{I}\mathbf{A} = \mathbf{A}\mathbf{I} = \mathbf{A}\) を満たす恒等元として機能する。
行列の加算は同型（行・列数が一致）な行列間でのみ定義され、各成分を単純に加算する。一方、乗法 \(\mathbf{C} = \mathbf{A}\mathbf{B}\) は、\(\mathbf{A}\) の列数と \(\mathbf{B}\) の行数が一致する場合にのみ成立し、結果行列の成分 \(c_{i,j}\) は以下のドット積で求められる。
\[ c_{i,j} = \sum_{k} a_{i,k} b_{k,j} \]行列乗算は結合則 \((\mathbf{AB})\mathbf{C} = \mathbf{A}(\mathbf{BC})\) および分配則を満たすが、一般的に交換則 \(\mathbf{AB} \neq \mathbf{BA}\) は成り立たない。この非可換性がアルゴリズム設計時の重要な制約となる。
繰り返し二乗法（バイナリ法）の基礎
スカラー値の冪乗 \(x^y\) を線形計算ではなく対数時間で求解する手法が繰り返し二乗法である。指数 \(y\) を2進数展開し、最下位ビットから順に評価することで計算回数を大幅に削減する。
例えば \(3^5\) を計算する場合、\(5 = (101)_2\) であるため処理は以下のように展開される。
- 1ビット目が1：累積値に \(3^1\) を乗算
- 基底を自乗：\(3 \to 3^2 = 9\)
- 2ビット目が0：乗算をスキップ
- 基底を自乗：\(9 \to 9^2 = 81\)
- 3ビット目が1：累積値に \(3^4 = 81\) を乗算
最終結果は \(3 \times 81 = 243\) となる。このロジックは剰余演算と組み合わせることで、オーバーフローを防ぎつつ巨大な冪乗を効率的に処理できる。
行列への拡張：高速冪乗アルゴリズム
前述の繰り返し二乗法における乗算処理を行列演算に置換することで、正方行列の \(n\) 乗を \(O(k^3 \log n)\)（\(k\) は行列の次元）で計算可能となる。実装上の重要な点は、累積値の初期化を単位行列で行うこと、および正方行列のみが自己乗算の対象となることである。
以下に、キャッシュ効率を考慮したループ構成と0ベースインデックスを採用した実装例を示す。
#include <vector>
using ll = long long;
using Matrix = std::vector<std::vector<ll>>;
constexpr ll MOD = 1000000007;
Matrix create_identity(size_t dim) {
Matrix I(dim, std::vector<ll>(dim, 0));
for (size_t i = 0; i < dim; ++i) I[i][i] = 1;
return I;
}
Matrix multiply(const Matrix& lhs, const Matrix& rhs) {
size_t r1 = lhs.size();
size_t c2 = rhs[0].size();
size_t inner = rhs.size();
Matrix res(r1, std::vector<ll>(c2, 0));
for (size_t i = 0; i < r1; ++i) {
for (size_t k = 0; k < inner; ++k) {
if (lhs[i][k] == 0) continue;
for (size_t j = 0; j < c2; ++j) {
res[i][j] = (res[i][j] + lhs[i][k] * rhs[k][j]) % MOD;
}
}
}
return res;
}
Matrix power_matrix(Matrix base, ll exponent) {
size_t dim = base.size();
Matrix acc = create_identity(dim);
while (exponent > 0) {
if (exponent & 1) acc = multiply(acc, base);
base = multiply(base, base);
exponent >>= 1;
}
return acc;
}
線形漸化式の高速求解への適用
行列高速冪乗の最も一般的な用途は、フィボナッチ数列に代表される線形漸化式の第 \(N\) 項を対数時間で求めることである。漸化式 \(F_n = F_{n-1} + F_{n-2}\) を行列の線形変換として再定式化する。
状態ベクトル \(\begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix}\) を考えるとき、次の関係式が導かれる。
\[ \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_{n-1} \\ F_{n-2} \end{pmatrix} \]遷移行列 \(\mathbf{T} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}\) と定義すると、初期状態からの推移は \(\mathbf{T}^{n-1}\) を適用することで一発で計算できる。
\[ \begin{pmatrix} F_{n} \\ F_{n-1} \end{pmatrix} = \mathbf{T}^{n-1} \begin{pmatrix} F_{1} \\ F_{0} \end{pmatrix} \]これにより、\(O(N)\) の逐次計算を \(O(\log N)\) の行列冪乗に置き換えることが可能となる。実装時にはインデックスのずれや \(N=1,2\) の境界条件に注意する必要がある。
#include <iostream>
#include <vector>
using ll = long long;
using Matrix = std::vector<std::vector<ll>>;
constexpr ll MOD = 1000000007;
Matrix mat_mult_2x2(const Matrix& a, const Matrix& b) {
Matrix c(2, std::vector<ll>(2, 0));
for (int i = 0; i < 2; ++i) {
for (int k = 0; k < 2; ++k) {
ll val = a[i][k];
if (val == 0) continue;
for (int j = 0; j < 2; ++j) {
c[i][j] = (c[i][j] + val * b[k][j]) % MOD;
}
}
}
return c;
}
Matrix mat_pow_2x2(Matrix base, ll p) {
Matrix res = {{1, 0}, {0, 1}};
while (p) {
if (p & 1) res = mat_mult_2x2(res, base);
base = mat_mult_2x2(base, base);
p >>= 1;
}
return res;
}
ll solve_fibonacci(ll n) {
if (n <= 0) return 0;
if (n == 1 || n == 2) return 1;
Matrix trans = {{1, 1}, {1, 0}};
Matrix final_mat = mat_pow_2x2(trans, n - 1);
return final_mat[0][0];
}
int main() {
ll target_n;
std::cin >> target_n;
std::cout << solve_fibonacci(target_n) << std::endl;
return 0;
}