행렬 거듭제곱

무엇을 가속하나

피보나치 \(f_n = f_{n-1} + f_{n-2}\) 같은 선형 점화식의 \(N\)번째 항을 단순
반복으로 구하면 \(O(N)\)이 든다. \(N\)이 \(10^{18}\)처럼 크면 불가능하다. 행렬
거듭제곱은 이를 \(O(K^3 \log N)\)으로 줄인다(\(K\)는 점화식의 차수).

점화식을 행렬로 옮기기

상태 벡터를 점화식이 참조하는 최근 항들로 잡는다. 피보나치라면
\(\begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix}\)를 한 단계 진행시키는 변환은

$$ \begin{pmatrix} f_{n+1} \\ f_n \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} f_n \\ f_{n-1} \end{pmatrix} $$

이다. 이 변환 행렬을 \(A\)라 하면, 초기 벡터에 \(A\)를 \(n\)번 곱한 것이 답이다.
즉 \(\begin{pmatrix} f_{n+1} \\ f_n \end{pmatrix} = A^n \begin{pmatrix} f_1 \\ f_0 \end{pmatrix}\).

일반화: \(f_n = c_1 f_{n-1} + \dots + c_k f_{n-k}\)이면 변환 행렬은 첫 행이
\((c_1, \dots, c_k)\)이고 그 아래에 항등 부분 행렬을 둔 동반 행렬이다.

거듭제곱을 빠르게: 분할 정복

행렬 거듭제곱 \(A^n\)은 정수 거듭제곱과 똑같이 분할 정복(이진 거듭제곱) 으로
계산한다.

$$ A^n = \begin{cases} (A^{n/2})^2 & n \text{ 짝수} \\ A \cdot A^{n-1} & n \text{ 홀수} \end{cases} $$

곱셈 횟수가 \(O(\log N)\)번이고, 한 번의 \(K \times K\) 행렬 곱이 \(O(K^3)\)이므로
전체 \(O(K^3 \log N)\)이다.

왜 옳은가

행렬 곱은 결합법칙을 만족하므로, \(A\)를 \(n\)번 곱하는 것은 어떤 순서로 묶어
곱해도 결과가 같다. 따라서 지수를 절반씩 나누는 분할 정복이 그대로 성립한다.
이는 정수 빠른 거듭제곱이 옳은 이유와 정확히 같은 원리다.

적용 범위

선형 점화식뿐 아니라 그래프에서 길이 \(L\) 경로의 수(인접 행렬의 \(L\)제곱),
상태 전이가 일정한 동적 계획법, 타일링·문자열 카운팅 등 "같은 선형
변환을 여러 번 반복"하는 모든 문제에 쓸 수 있다. 답이 크면 보통 모듈러
\(10^9 + 7\)로 계산한다.

모듈러 행렬 거듭제곱 레퍼런스

크기 \(K\)의 정사각 행렬을 vector<vector<ll>>로 두고, 곱셈과 이진 거듭제곱을
구현한다. 모듈러를 곱셈마다 적용하지 않으면 오버플로가 난다.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll MOD = 1000000007;

typedef vector<vector<ll>> Mat;
int K;

Mat mul(const Mat& A, const Mat& B) {
    Mat C(K, vector<ll>(K, 0));
    for (int i = 0; i < K; i++)
        for (int k = 0; k < K; k++) {
            if (A[i][k] == 0) continue;        // 0이면 건너뛰어 가속
            for (int j = 0; j < K; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
        }
    return C;
}

Mat identity() {
    Mat I(K, vector<ll>(K, 0));
    for (int i = 0; i < K; i++) I[i][i] = 1;
    return I;
}

Mat power(Mat A, ll n) {
    Mat result = identity();
    while (n > 0) {
        if (n & 1) result = mul(result, A);    // 홀수 비트면 곱한다
        A = mul(A, A);                         // 제곱
        n >>= 1;
    }
    return result;
}

int main() {
    ll n; scanf("%lld", &n);
    K = 2;
    Mat A = {{1, 1}, {1, 0}};
    Mat An = power(A, n);
    // f_n = (A^n)[0][1] * f_1 + (A^n)[1][1] * f_0, 여기선 f_0=0, f_1=1
    printf("%lld\n", An[0][1] % MOD);
    return 0;
}

흔한 실수

• 모듈러 누락 — A[i][k] * B[k][j]는 두 수가 각각 \(10^9\) 근처면 \(10^{18}\)가
  되어 ll에 겨우 들어간다. 누적합에 반드시 % MOD를 적용해야 한다.
• 항등 행렬 초기값 — 거듭제곱의 시작은 단위 행렬이다. 영행렬로 시작하면
  결과가 모두 0이 된다.
• 초기 벡터 정렬 — 상태 벡터의 항 순서와 변환 행렬의 행렬 정의가 일치해야
  한다. \(f_n\)이 위인지 아래인지 한 번 더 확인하자.
• 음수 모듈러 — 점화식에 음수 계수가 있으면 ((x % MOD) + MOD) % MOD로
  보정한다.

응용 패턴

상수항이 붙은 점화식 \(f_n = a\ f_{n-1} + b\) 같은 경우, 상태 벡터에 항상 1을
실어 두는 칸을 하나 더 추가하면 affine 변환을 선형 행렬로 표현할 수 있다.
차수가 더 큰 점화식이나 \(N\)이 매우 클 때 행렬 곱 자체가 부담되면, 이후 배울
키타마사법으로 \(O(K^2 \log N)\)까지 줄인다.