3개월 정도 고민한 문제를 드디어 알았다.

머릿말

여름방학에 이 문제를 보고 너무 풀고 싶었고, 풀이도 찾았지만 그 당시에는 풀이가 이해가 안됬었다.

하지만, 3개월만에 이 문제에 대한 해답을 확실히 알게 되었으니 포스팅 해보자!!

문제 본문

BOJ 13716 : 피보나치 수열처럼 보이지만… Codeforces Round #392 (Div1) : C. Yet Another Number Sequence

피보나치 수열처럼 보이지만…

피보나치 수열은 다음과 같이 정의된다.

\[F_1 = 1, F_2 = 2, F_i = F_{i-1} + F_{i-2} (i > 2).\]

$A_i(k)$는 다음과 같이 정의 된다.

\[A_i(k) = F_i \times i^k (i \ge 1).\]

n과 k가 주어졌을 때, $A_1(k) + A_2(k) + … + A_n(k)$를 구하는 프로그램을 작성하시오.

입력

첫째 줄에 $n$과 $k$가 주어진다. ( $1 \le n \le 10^{17}, 1 \le k \le 40 $ )

출력

첫째 줄에 $A_1(k) + A_2(k) + … + A_n(k)$를 $10^9 + 7$로 나눈 나머지를 출력한다.

입출력 예시

예제 입력 1

1 1

예제 출력 1

1

예제 입력 2

4 1

예제 출력 2

34

예제 입력 3

5 2

예제 출력 3

316

예제 입력 4

7 4

예제 출력 4

73825

문제 풀이

우선 N의 범위를 보았을 때, N을 log수준으로 떨어뜨리거나 k를 통한 연산으로 반복해야한다.

또한 경험상 행렬곱 또는 DP를 이용하여 연산할 수 있다는 가정을 하였다.

1. 기본식 분해

우선 주어진 n,k에 대해 다음과 같은 정의를 내리자.

\[S_{n,k} = \sum_{i=1}^{n}A_i(k) = \sum_{i=1}^{n}F_i \times i^k\]

그렇다면 기본적인 점화식은 다음과 같다.

\[S_{n,k} = F_n \times n^k + \sum_{i=1}^{n-1}A_i(k)\]

만약 이걸 그대로 연산한다면 시간복잡도는 $O(n)$이므로 TLE다. 어떤식으로 추가적인 연산을 할 수 있을까.

2. 식 따로 생각하기

$(i+1)^k$은 $i^0$~$i^k$까지 다음과 같이 나타낼 수 있다.

\[(i+1)^k = C^0_ki^0 + C^1_ki^1 ... C^k_ki^k\]

그렇다면 유사-피보나치에 대입을 해보면 다음과 같다.

\[\begin{align} F_{q+1}\times (q+1)^k & = F_q\times (q+1)^k + F_{q-1}\times (q+1)^k \\ & = F_q \times (C^0_kq^0 + C^1_kq^1 ... C^k_kq^k) + F_{q-1} \times (C^0_kq^0 + C^1_kq^1 ... C^k_kq^k) \end{align}\]

여기서 어떻게 묶어야 행렬곱으로 표현하고, 연산을 효율적으로 할 수 있을까.

3. 행렬로 변환하기

피보나치와 거듭제곱의 결과를 $1 * k+1$ 행렬로 나타내보자.

\[FL_{i,j,k} = [F_ij^0 F_ij^1 F_ij^2 ... F_ij^k]\]

또한 $C_m^n$을 전처리하여 다음과 같은 행렬을 만들자.

\[C_k = \begin{bmatrix} C_0^0 & C_1^0 & C_2^0 & \cdots & C_k^0 \\ 0 & C_1^1 & C_2^1 & \cdots & C_k^1 \\ 0 & 0 & C_2^2 & \cdots & C_k^2 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & C_k^k \end{bmatrix}\]

4. 점화식 구체화하기

다음과 같이 정의한다면 아래와 같은 수식이 만족한다.

\[FL_{i,i,k} \times C_k = FL_{i,i+1,k}\]

이 식을 이용해서 우리가 구해야하는 점화식을 구해보면 다음과 같다.

\[S_{i,k} = S_{i-1,k} + FL_{i,i,k}[k]\] \[FL_{i,i,k} = FL_{i-2,i-1,k} \times C_k + FL_{i-1,i-1,k} \times C_k\]

위 식을 이용하여 행렬곱으로 나타내면 다음과 같다.

\[\begin{bmatrix} FL_{i-1,i,k} & FL_{i,i,k},S_{i-1,k}\end{bmatrix} \times \begin{bmatrix} 0 & C_k & 0 \\ C_k & C_k & I_k \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} FL_{i,i+1,k} & FL_{i+1,i+1,k},S_{i,k}\end{bmatrix}\]

단 여기서 $I_k$는 다음과 같다.

\[I_k = [0,0, ... 1 ] ^ t\]

위 행렬곱을 수정하면 다음과 같이 바꿀 수 있다.

\[\begin{bmatrix} FL_{0,1,k} & FL_{1,1,k},S_{0,k}\end{bmatrix} \times \begin{bmatrix} 0 & C_k & 0 \\ C_k & C_k & I_k \\ 0 & 0 & 1 \end{bmatrix} ^ n = \begin{bmatrix} FL_{n,n+1,k} & FL_{n+1,n+1,k},S_{n,k}\end{bmatrix}\]

초기 값은 정의에 따라 다음과 같다.

\[FL_{0,1,k} = FL_{1,1,k} = \begin{bmatrix}1 & 1 & \cdots & 1\end{bmatrix}\] \[S_{0,k} = 0\]

5. Source Code

구조체와 연산자 오버로딩으로 구현하였다. power(n)은 $O(log n)$으로 구현하였다.

구현할때, 배열크기를 잘못잡아 Runtime-Error를 한번 받았다.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9+7;

struct matrix{
    int N;
    ll A[82][82];
    matrix(int n){
        N = n;
        for(int i = 0 ; i < n ; i++){
            for(int j = 0 ; j < n ; j++){
                A[i][j] = 0;
            }
        }
    }

    matrix operator +(matrix t){
        matrix res(N);
        for(int i = 0 ; i < N ; i++){
            for(int j = 0 ; j < N ; j++){
                res.A[i][j] = (A[i][j] + t.A[i][j]) % mod;
            }
        }
        return res;
    }

    matrix operator *(matrix t){
        matrix res(N);
        for (int i=0;i<N;i++){
            for (int j=0;j<N;j++){
                for (int k=0;k<N;k++){
                    res.A[i][j] = (res.A[i][j] + A[i][k] * t.A[k][j]) % mod;
                }
            }
        }
        return res;
    }
};

ll bi[44][44];

int main(){
    ll N; int K; cin >> N >> K;
    for (int i = 0 ; i <= K ; i++){
        bi[i][0] = bi[i][i] = 1;
        for (int j = 1 ; j < i ;j++) bi[i][j] = (bi[i-1][j-1] + bi[i-1][j]) % mod;
    }

    matrix A(2*(K+1)+1);
    for (int i = 0 ; i <= K ; i++){
        for (int j = i ; j <= K ; j++){
            A.A[i][j+K+1] = A.A[i+K+1][j] = A.A[i+K+1][j+K+1] = bi[j][i];
        }
    }
    A.A[2*K+2][2*K+2] = A.A[2*K+1][2*K+2] = 1;

    matrix R(2*(K+1)+1);
    for(int i = 0 ; i <= 2*K+2 ; i++) R.A[i][i] = 1;

    while (N){
        if (N & 1) R = R * A;
        A = A * A;
        N >>= 1;

    }
    ll ans = 0;
    for (int i = 0 ; i < 2*K+2 ; i++) ans = (ans + R.A[i][2*K+2]) % mod;
    cout << ans;
}

맺음말

처음 풀이는 코드포스에서 크리님 코드를 통해 봤다. 연산자 오버로딩 말고는 왜 이렇게 되는지 몰라서 넘어갔었다. 연산자 오버로딩은 그대로 가져왔기에 어떻게 보면 카피같기도 하고…

크리님 풀이는 2의 power를 이용하여 행렬을 다른 방식으로 곱하여 더하는데 어떤 방식인지 이해하지 못했다.

13716 풀이

출처 : BOJ 문제에 대한 모든 권리는 BOJ(acmicpc.net, startlink)에 있음

Leave a Comment