[백준 C++] 피사노 주기와 행렬 제곱을 활용한 피보나치의 수 구하기

피보나치 수 3

피보나치 수 6

백준 블로그를 참고했다.

N 번째의 피보나치 수를 구해야 할 때 N 의 수가 충분히 작다면 재귀함수나 반복문, 메모이제이션을 활용해서 문제를 풀면 된다. 하지만 '피보나치 수 3' 나 '피보나치 수 6' 처럼 N의 범위가 1,000,000,000,000,000,000 까지라면 위에서 언급한 방법으로는 시간을 맞출 수 없다. 그래서 두가지 방법을 활용할 수 있는데 하나는 '파사노 주기이며 다른 하나는 행렬의 제곱이다.

피보나치 수 3

피보나치 수 3은 N의 범위가 1,000,000,000,000,000,000 까지고 결과값을 1,000,000으로 나눈 값을 출력해야 한다. 피보나치 수를 M 으로 나눈 수는 항상 주기를 가지게 되는데, 이를 피사노 주기(Pisano Period)라고 한다. 예를 들면 다음과 같다.

N	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19
Fibo	0	1	1	2	3	5	8	13	21	34	55	89	144	233	377	610	987	1597	2584	4181
M=2	0	1	1	0	1	1	0	1	1	0	1	1	0	1	1	0	1	1	0	1
M=3	0	1	1	2	0	2	2	1	0	1	1	2	0	2	2	1	0	1	1	2
M=4	0	1	1	2	3	1	0	1	1	2	3	1	0	1	1	2	3	1	0	1

피보나치 수를 M 으로 나눈 나머지는 빨간색 숫자들이 반복된다. 이 원리를 활용하면 연산을 적게 하면서도 나머지 값을 구할 수 있다. 이 문제의 경우 피사노 주기를 상대적으로 쉽게 구할 수 있는데, M 의 값이 $10^k(k>2)$라면 피사노 주기는 항상 $15*10^{k-1}$이다. 만약 M 이 다른 값이라면 코드를 작성해서 피사노 주기를 구할 수 있다. 피사노 주기를 구하는 코드는 다음과 같다.

long long get_pisano_period(long long m) {
  long long c, a = 0, b = 1;
  for (int i = 0; i < m * m; i++) {
    c = (a + b) % m;
    a = b;
    b = c;
    if (a == 0 && b == 1) return i + 1;
  }
}

아래는 피보나치 수 3 소스코드다. 반복문은 n 을 피사노 주기로 나눈 값 만큼만 돌리면 된다. 계산을 마치고 MOD 로 나눠주는 것을 잊으면 안된다.

#include <iostream>
using namespace std;
#define MOD 1000000
const int pisano = MOD / 10 * 15;
int arr[pisano] = {0, 1};
int main() {
  long long n;
  cin >> n;
  long long j = n % pisano;
  for (int i = 2; i <= j; i++) {
    arr[i] = arr[i - 1] + arr[i - 2];
    arr[i] %= MOD;
  }
  cout << arr[j] << '\n';
  return 0;
}

피보나치 수 6

피보나치 수 3 은 피사노 주기를 활용해서 풀었다면 피보나치 수 6 에서는 그것이 불가능하다. 3 에서는 M 의 값이 상대적으로 작은 1,000,000이었고, 피사노 주기도 1,500,000 에 불과했기 때문에 빠른 시간 안에 답을 도출해낼 수 있었다. 하지만 6에서는 결과값을 1,000,000,007 로 나눈 나머지를 찾아야 하는데, 이는 피사노 주기를 찾더라도 그 값이 너무 커서 반복문으로 피사노 주기까지의 값을 계산하면 시간이 오래 걸린다. 그래서 이 문제는 행렬의 제곱을 활용해서 풀어야 한다.

피보나치 수를 행렬로 나타내면 다음과 같다.

$$\binom{F_{n+2}}{F_{n+1}} = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix} + \binom{F_{n+1}}{F_n}$$

이 식을 정리해보면 다음과 같아진다.

$$\begin{pmatrix} F_{n+1} & F_n\\ F_n & F_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}^n$$

#include <iostream>
#include <vector>
using namespace std;
 
typedef long long ll;
typedef vector<vector<ll>> matrix;
const ll mod = 1000000007LL;
 
matrix operator * (const matrix &a, const matrix &b) {
  int n = a.size();
  matrix c(n, vector<ll>(n));
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      for (int k = 0; k < n; k++) {
        c[i][j] += a[i][k] * b[k][j];
      }
      c[i][j] %= mod;
    }
  }
  return c;
}
 
 
int main() {
  cin.tie(NULL);
  ios::sync_with_stdio(false);
 
  ll n;
  cin >> n;
 
  if (n < 2) {
    cout << n << '\n';
    return 0;
  }
 
  matrix ans = {{1, 0}, {0, 1}};
  matrix a = {{1, 1}, {1, 0}};

  while (n > 0) {
    if (n % 2) ans = ans * a;
    a = a * a;
    n /= 2;
  }
  cout << ans[0][1] << '\n';
 
  return 0;
}

백준 블로그에는 코드만 덩그라니 있어서 다른 블로그를 참고했다. matrix c(n, vector<ll>(n)); 이런 코드는 굉장히 낯설었는데 행렬 c 는 vector<long long>(n) 를 원소로 갖는 벡터 n 개라는 뜻이라고 한다.

아래 코드는 제곱 수를 계산하는 코드이다. 이항계수에서 제곱을 구할 때 사용했던 것과 원리가 비슷하다. a 는 {{1, 0}, {0, 1}} 인데, 이 값은 일반 숫자로 치면 1에 해당하는 것으로 곱셈을 할 때 초기값으로 주는 것이다. a 에는 자기 자신을 곱해주고, n 은 2로 나눠준다. n 을 2 로 나눈 나머지가 1 일 때는 정답(ans)에 a 를 곱해야 한다.

while (n > 0) {
  if (n % 2) ans = ans * a;
  a = a * a;
  n /= 2;
}