[Baekjoon,10830] 행렬 제곱

일단 행렬의 곱셈을 모른다면 winner blog : 행렬의 연산(곱셈)을 참고하자.

문제에서 행렬을 제곱해야 하는데 이 제곱하는 양이 엄청 크다. 따라서 일반적인 방법으로 하나씩 곱해나가면 주어진 시간에 해결할 수 없다.

// 정방 행렬을 표현하는 클래스가 있다고 가정하자.
class Matrix;
// ret은 항등 행렬으로 초기화 한다. 
Matrix A,ret;
// 1000억번의 연산을 감당할 수 없다.
for(long long i=0;i<100000000000;++i)
  ret*=A;

---

분할 정복

이를 해결하는 방법은 제곱을 분할하여 계산하는 것이다. 정방 행렬 $A$를 $p$번 제곱한다고 했을 때 $p$를 절반으로 나누어 보자.

• $ A^p = A^{p/2} \times A^{p/2} $

이렇게 진행하면 1000억번의 연산의 크기가 대략 36번으로 확 줄어든다.

이를 재귀적으로 함수를 구현할 수 있다. 기저 사례를 제곱근이 0이 되었을 때로 설정하였다. 이때 기저 사례로 반환하는 값은 아래와 같은 행렬을 0 제곱한 값인 항등 행렬(identity matrix)이다.

• $ \begin{bmatrix} 1 & 0 & 0 \\
  0 & 1 & 0 \\
  0 & 0 & 1 \end{bmatrix} $

왜 행렬을 $0$ 제곱한 값이 항등 행렬인지 간단하게 설명하면 수에서 $1$은 아무리 거듭 제곱해도 $1$이다. 즉 수의 곱셈에 대한 항등원은 $1$이다. 마찬가지로 행렬인 자신을 아무리 거듭 제곱해도 자기 자신인, 즉 행렬의 곱셈에 대한 항등원은 항등 행렬이다. 이에 대한 자세한 설명은 수학방에서 볼 수 있다.

// 행렬의 거듭 제곱
Mat pow(const Mat& mat, long long p)
{
    // 기저 사례: 항등 행렬 반환
    if(p==0)
    {
        Mat temp(n,vector<int>(n));
        for(int i=0;i<n;++i)
            temp[i][i]=1;
        return temp;
    }
    if(p&1)
        return mat*pow(mat,p-1);
    Mat half=pow(mat,p/2);
    return half*half;
}
// 코드 참고: 구종만, 프로그래밍 대회에서 배우는 알고리즘 문제해결전략, 179p, 인사이트 

---

반복문으로 분할 정복

직관적이지는 않지만 반복문으로도 구현이 가능하다.

Mat pow(Mat mat, long long p)
{
    Mat ret(n,vector<int>(n));
    for(int i=0;i<n;++i)
        ret[i][i]=1;
    while(p)
    {
        if(p&1)
            ret*=mat;
        p>>=1;
        mat*=mat;
    }
    return ret;
}

---

구현

이 문제를 해결하는 모든 것을 구현하면 아래와 같다. 따로 클래스를 만들지 않고 이중 vector를 이용했으며, 연산자 *를 오버로딩하여 직관적으로 코드를 구현했다.

#include <iostream>
#include <vector>
using namespace std;

typedef vector<vector<int> > Mat;

const int MOD=1000;

int n;

Mat operator*(const Mat& a, const Mat& b)
{
    Mat ret(n,vector<int>(n,0));
    for(int i=0;i<n;++i)
        for(int j=0;j<n;++j)
        {
            for(int k=0;k<n;++k)
                ret[i][j]+=a[i][k]*b[k][j];
            ret[i][j]%=MOD;
        }
    return ret;
}

// 행렬의 거듭 제곱
Mat pow(const Mat& mat, long long p)
{
    // 기저 사례: 항등 행렬 반환
    if(p==0)
    {
        Mat temp(n,vector<int>(n));
        for(int i=0;i<n;++i)
            temp[i][i]=1;
        return temp;
    }
    if(p&1)
        return mat*pow(mat,p-1);
    Mat half=pow(mat,p/2);
    return half*half;
}

int main()
{
    long long p;
    cin>>n>>p;
    
    Mat mat(n,vector<int>(n));
    for(int i=0;i<n;++i)
        for(int j=0;j<n;++j)
            cin>>mat[i][j];
    
    Mat ret=pow(mat,p);
    for(int i=0;i<n;++i)
    {
        for(int j=0;j<n;++j)
            cout<<ret[i][j]<<' ';
        cout<<endl;
    }

    return 0;
}