BAEKJOON 12728: n제곱 계산

 

2008년 Google code jam round A에 나왔던 문제이다. 이 문제에서는 매우 정밀한 수의 연산을 요구하기 때문에 double이나 float, Decimal 등을 이용한 연산이 옳은 결과를 내지 못하도록 설계되어 있다.

 

Solution

먼저, \( (3 + \sqrt{5})^{n} \)에 대해서 생각해 보자. \(\sqrt{5}\)가 살짝 거슬린다. 이를 해결하기 위해 \( (3 - \sqrt{5})^{n}\)을 더해준다. 그러면 \(\sqrt{5}\)가 존재하는 모든 항을 없앨 수 있다. 조금 더 자세히 설명하자면,

\( (3+\sqrt{5})^{n}+(3-\sqrt{5})^{n}=\sum_{i=0}^{n} (3)^{n-i}(\sqrt{5})^{i}+\sum_{i=0}^{n} (3)^{n-i}(\sqrt{5})^{i}(-1)^{i}\)이고, 이는 \( \sum_{i=0}^{\lfloor \frac{n}{2} \rfloor} (3)^{n-2i}(5)^i\)로 표현할 수 있다.

 

이제, \((3+\sqrt{5})^n\)을 정리하였을 때를 생각해 보자. \(\sqrt{5}\)가 붙은 부분과 그렇지 않은 부분으로 나뉠 것이다. 그렇다면 유리수 부분을 \(a_n\), 무리수 부분을 \(b_{n} \sqrt{5}\)라 하자.

 

이 포스팅에서 처음 도출한 식인 \( \sum_{i=0}^{\lfloor \frac{n}{2} \rfloor} (3)^{n-2i}(5)^i\)으로 돌아가자. 이는 \(2a_n\)과 같다. 그리고 \((3-\sqrt{5})^n\)는 항상 (0,1]의 범위에 들어간다. 그러므로, 우리가 구하고자 하는 \((3-\sqrt{5})^n\)의 값은 \(2a_{n}-1\)과 같다. 이를 이용하여 답을 구해줄 수 있다.

 

이제 \(a_n\)을 구해보자.

\((3+\sqrt{5})^{n+1}=(3+\sqrt{5})(a_{n}+b_{n} \sqrt{5})=(3a_{n}+5b_{n})+(a_{n}+3b_{n})\sqrt{5}\)이다.

그러므로, \(a_{n+1}=3a_{n}+5b_{n}\), \(b_{n+1}=a_{n}+3b_{n}\) 가 된다.

 

이를 행렬로 표현하여 나타냄으로서 빠르게 \(a_n\)을 구할 수 있다. 위 점화식을 행렬로 나타내면

\( \begin{pmatrix} a_n \\ b_n  \end{pmatrix} = \begin{pmatrix}3&5 \\ 1&3 \end{pmatrix}^n \begin{pmatrix} a_0 \\b_0 \end{pmatrix}\) 이다. 여기서 \((3+\sqrt{5})^0=1\) 이므로, \(a_0=1\), \(b_0=0\) 이다.

 

행렬의 거듭제곱은 분할 정복을 통해 \(O(\log(n))\)에 수행하여 문제를 해결할 수 있다.

 

정답 코드

#include <bits/stdc++.h>
using namespace std;
using ll=long long;
const int mod=1000;
void solve(){
    int n;cin>>n;
    int A[2][2]={{3,5},{1,3}},B[2][2]={{3,5},{1,3}};n--;
    while(n){
        int a,b,c,d;
        if(n&1){
            a=(B[0][0]*A[0][0]+B[0][1]*A[1][0])%mod;
            b=(B[0][0]*A[0][1]+B[0][1]*A[1][1])%mod;
            c=(B[1][0]*A[0][0]+B[1][1]*A[1][0])%mod;
            d=(B[1][0]*A[0][1]+B[1][1]*A[1][1])%mod;
            B[0][0]=a;B[0][1]=b;B[1][0]=c;B[1][1]=d;
        }
        a=(A[0][0]*A[0][0]+A[0][1]*A[1][0])%mod;
        b=(A[0][0]*A[0][1]+A[0][1]*A[1][1])%mod;
        c=(A[1][0]*A[0][0]+A[1][1]*A[1][0])%mod;
        d=(A[1][0]*A[0][1]+A[1][1]*A[1][1])%mod;
        A[0][0]=a;A[0][1]=b;A[1][0]=c;A[1][1]=d;
        n>>=1;
    }
    ll x=(B[0][0]*2+999)%mod;
    if(x<10)cout<<"00";
    else if(x<100)cout<<"0";
    cout<<x;
}
int main(){
    ios_base::sync_with_stdio(false);cin.tie(NULL);
    int T;cin>>T;
    for(int i=1;i<=T;i++){
        cout<<"Case #"<<i<<": ";
        solve();
        cout<<"\n";
    }
    return 0;
}