백준 13246, 행렬 제곱의 합

개요

문제 링크
골드 2, 분할정복
아주 큰 지수에 대해 행렬의 제곱의 합 구하기

---

접근

1. 분할 정복을 이용한 거듭제곱은 선형 대수학의 존재 이유인 선형 동차방정식을 좀 더 컴퓨터에 적합한 과정에서 접근할 수 있게 해주는 방법론이라고 생각하면 될 것 같다. 예를들면 피보나치 수열의 일반항도 선형 조합이기 때문에 점화식을 잘 해결해주면 행렬표현이 가능한데, 행렬 제곱의 계산에서 분할정복을 이용하면 log 수준에서 해결할 수 있다.
   개인적으로 공학 수학의 근간은 선형 대수학이라 생각하는데, 대부분의 공학적인 방법을 선형 다항식으로 퉁쳐버리고 행렬이라는 무지막지한 놈을 이용해 N차 방정식으로 끙끙매는 공대생들에게 사칙연산급의 프리패스 알사탕을 던져주기 때문이다. 여기에 C++의 미친 장점 operator가 추가된다면? 컴퓨터로 기초적인 공학의 정복이 가능하다.

2. 선형대수학을 좋아해 여담이 길었는데, 문제로 돌아와서, 제곱의 합을 분할정복으로 접근하려면 우선 식을 절반으로 나누어야 한다. 정수를 예시로 들어보자.

   2 + 2^2 + 2^3 + 2^4 + 2^5 + 2^6 + 2^7 + 2^8
   = (2 + 2^2 + 2^3 + 2+4) + 2^4 * (2 + 2^2 + 2^3 + 2^4)

   겹치는 부분을 발견할 수 있다. 그러면 연산을 절반어치만큼만 해주면 된다. 이걸 재귀적으로 이용하면? logN번의 연산으로 문제를 해결할 수 있다. B가 10^11인데 log(10^11)을 하면 40번 안에 해결이 된다. 만약 지수가 홀수라면? 2+2×(half + 2^4 × half)를 해주면 된다. 비슷하게 작동만 하면 아무렇게나 써도 된다.

3. 이때 위 식에서 오른쪽 덩어리에 곱해주는 2^4는 따로 계산해주어야 하는데, 이것도 분할정복을 이용한다. 거듭제곱의 분할정복은 매우 유명한데, 아래와 같이 구현한다.

   power(n,k) {
       v = 1
       while(k)
           if (k%2)
               v*=n
           n*=n, k/=2
       return v
   }

   지수인 k를 절반씩 나누면서 진행하는데, k가 홀수일 때만 지금까지 곱해온 덩어리를 출력값에 곱해준다. 원리는 비트마스킹을 생각하면 쉬운데, k를 이진수로 나타내고 bit가 있는 부분마다 곱을 해주는 것이다.

4. 나머지는 행렬곱의 개념을 알고 있으면 크게 어렵지 않다. 앞서 말했듯이 C++의 operator를 이용해 코드를 단순화할 수 있다.

---

Pseudo code

power(a, k) {
    v = a, k--
    while(k)
        if (k%2)
            v *= a
        a *= a, k /= 2
    return v
}

sum(a, b) {
    if (b = 1) return a
    half = sum(a, b/2)
    if (b%2)
        return a+a*(half+power(a,b/2)*half)
    return half+power(a,b/2)*half
}

---

Source code

#include <iostream>
using namespace std;

using ld = long long;
int R = 1000;
ld n,b;
struct mat {
    int a[5][5];
    mat operator + (mat b) {
        mat c;
        for (int i=0; i<n; i++)
            for (int j=0; j<n; j++) {
                c.a[i][j] = a[i][j]+b.a[i][j];
                c.a[i][j] %= R;
            }
        return c;
    }
    mat operator * (mat b) {
        mat c;
        for (int i=0; i<n; i++)
            for (int j=0; j<n; j++) {
                c.a[i][j] = 0;
                for (int k=0; k<n; k++)
                    c.a[i][j] += a[i][k]*b.a[k][j];
                c.a[i][j] %= R;
            }
        return c;
    }
    mat pow(ld k) {
        mat a = *this;
        mat v = a;
        k--;
        while (k) {
            if (k%2)
                v = v*a;
            a = a*a, k/=2;
        }
        return v;
    }
};

mat sum(mat a, ld b) {
    if (b==1) return a;
    mat c = sum(a,b/2);
    if (b%2)
        return a+a*(c+a.pow(b/2)*c);
    else
        return c+a.pow(b/2)*c;
}

mat a;
int main() {
    cin>>n>>b;
    for (int i=0; i<n; i++)
        for (int j=0; j<n; j++) {
            cin>>a.a[i][j];
            a.a[i][j] %= R;
        }
    mat c = sum(a,b);
    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++)
            cout << c.a[i][j] << " ";
        cout << "\n";
    }
}