Berlekamp-Massey 알고리즘은 특정한 DP의 점화식을 찾아주는 알고리즘이다. $10^{18}$ 번째 피보나치 수를 찾기 위해서 행렬 곱셈을 짜고, 타일 채우기 문제를 풀기 위해서 수많은 점화식과 씨름하던 옛 시간은 이젠 안녕. 이제는 백트래킹 짜고 하드코딩해서 넣으면 끝난다.

이 글은 알고리즘의 구현법, 동작 원리나 증명에 대해서 거의 설명하지 않는다. 그 이유는 내가 구현법과 동작 원리, 증명을 모르기 때문이다. 알고리즘 구현은 여기에서 복붙해서 사용하면 된다. 이론적 배경지식이 상당히 깊지만, 그 활용도가 매우 높기 때문에, 일단 이해하지 말고 작동법부터 제대로 깨우친 후, 나중에 다시 돌아와서 방법을 이해하는 것을 추천한다.

1967년 이 알고리즘을 개발한 수학자 Elwyn Berlekamp가 최근 (2019년 4월 9일) 사망했습니다. 고인의 명복을 빕니다. A final game with Elwyn Berlekamp

링크에 주어진 코드에는 크게 네 가지 부분이 있다.

• berlekamp_massey 함수. $n$개의 DP 결과값이 주어졌을 때, Berlekamp-Massey 알고리즘을 사용하여 해당 DP의 점화식을 찾아주는 루틴. $O(n^2)$ 시간 복잡도에 작동한다 (정확히는, 답이 되는 점화식의 크기가 $k$ 일때 $O(nk+n\log mod)$)
• get_nth 함수. 주어진 DP 점화식을 사용하여 $n$번째 항을 구하는 루틴.
• get_min_poly 함수. 3번 챕터에서 설명한다. 지금은 몰라도 상관없다.
• det 함수. 4번 챕터에서 설명한다. 지금은 몰라도 상관없다.

berlekamp_massey 함수에 구현된 Berlekamp-Massey는 입력 수열을 받아서, 위 수열을 만족하는 가장 짧은 점화식 $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 을 찾는다. 이 때 DP 점화식의 형태는 정확히 저러한 형태여야 하며, 이러한 점화식을 선형 점화식 (Linear recurrence) 라고 부른다. 예를 들어, $F_n = F_{n-1} + F_{n-2}$ 라는 점화식을 가지는 피보나치 수열, $P_n = P_{n-1} + P_{n-5}$ 라는 점화식을 가지는 파도반 수열 등이 선형 점화식의 예이다. $C_n = \sum_{i = 0}^{n-1}{C_i \times C_{n-i-1}}$ 형태의 점화식을 가지는 카탈란 수열은 선형 점화식의 예가 아니기 때문에 이 알고리즘을 통해서 구할 수 없다. 어떠한 패턴의 점화식이 계산되는지는 아래에서 계속 설명한다.

get_nth 함수는 Berlekamp-Massey 알고리즘과 전혀 상관이 없는 코드로, 점화식이 주어졌을 때 이 점화식에서 $DP_n$ 이 무엇인지를 계산해 주는 함수이다. 이 함수의 시간 복잡도는 $O(k^2 \log n)$ 이다. 일반적인 DP 점화식을 계산하는 데 드는 시간 복잡도는 $O(nk)$ 이고, 실질적으로 $k \le n$ 이니, 대부분의 경우 충분히 빠른 알고리즘이다.

get_nth 함수는 흔히 키타마사법이라고 불리는 방법으로 구현되어 있다. 흔히 저러한 DP 점화식을 구하는 데는 행렬을 사용하여 $O(k^3 \log n)$ 에 구하는 방법이 잘 알려져 있으나, 키타마사법을 사용하면 더 빠른 시간 복잡도인 $O(k^2 \log n)$ 에 이를 구할 수 있다. 키타마사법에 대해서는 이 동영상 이 블로그 에 잘 설명되어 있으나, 여기서도 한번 더 욕심을 내서 설명해 보겠다. (관심이 없으면 넘겨도 이후 글을 읽는데 상관이 없다.)

키타마사법은 생성함수와 비슷한 형태로 이해를 하면 상당히 간단하다. $DP_n$ 을 구하는 과정을 생각해 보자. $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 으로 대체(substitution)해주면, 더 작은 DP값들인 $DP_{n-1}, DP_{n-2}, \cdots, DP_{n-k}$ 에 대한 식으로 나타내어 줄 수 있다. $DP_{n-1}$ 을 또 대체해주면 $DP_{n-2} \cdots DP_{n-k-1}$ 가 되고, $DP_{n-2}, DP_{n-3}$ 등을 점화식을 사용하여 계속 똑같은 변환 방법으로 대체해주면 결국 언젠가는 $C_0 \times DP_0 + C_1 \times DP_1 + \cdots + C_{k-1} \times DP_{k-1}$ 형태가 되어서 문제를 풀 수 있을 것이다. 흥미로운 관찰은, 여기서 최대항을 대체한다는 것은 일반적으로 수학에서 다항식에 나머지 연산을 취하는 것과 동일하다는 것이다. $DP_n$ 을 $x^n$ 이라고 보자. $x^n$ 을 $\sum_{i=1}^{k}{C_i x^{n-i}}$ 로 대체하는 것을 반복하는 연산은, 사실 $x^n$ 을 $(x^k - \sum_{i=0}^{k-1}{C_{i+1} \times x^i})$ 으로 나눈 나머지를 구하는 것과 동일하다.

$x^n$ 을 단순히 $(x^k - \sum_{i=0}^{k-1}{C_{i+1} \times x^i})$ 으로 나눈 나머지를 구하는 것은 물론 효율적이지 않다. 이를 단순히 구현하면 $O(nk)$ 의 시간 복잡도가 들어서 일반적인 DP 계산과 차이가 없다 (아마 코드도 일반적인 DP 계산과 별로 다르지 않을 것이다). 그런데, 일반적인 수에 대해서 $x^n \mod p$ 을 $O(\log n)$ 시간에 구하는 다음과 같은 빠른 알고리즘이 잘 알려져 있다:

• $x^{2k+1} = x^{2k} \times x \mod p$
• $x^{2k} = x^k \times x^k \mod p$

여기서 $p$ 를 단순히 다항식으로 치환해주면, $O(\log n)$ 번의 나머지 연산과 곱셈만으로 문제를 해결할 수 있다. 다항식에 mod가 되어 있으니 항의 크기는 항상 $k$ 이하임이 보장되고, 이 경우 두 다항식의 곱셈과 나눗셈은 모두 $O(k^2)$ 에 구현할 수 있다. 고로 $O(k^2 \log n)$ 에 문제를 풀 수 있고, 다항식 곱셈과 나눗셈을 $O(k \log k$) 에 구현하면 복잡도가 $O(k\log k\log n)$ 으로 줄어든다. (곱셈은 FFT를 사용하면 매우 간단하고, 나눗셈 역시 FFT를 사용하면 되나 간단하지는 않으니 관심이 있으면 찾아보도록 하자.)

모든 것을 알았다면 Berlekamp-Massey 알고리즘을 사용하는 것은 간단하다. 예를 들어, 점화식이 $F_n = F_{n-1} + F_{n-2}$ 인 피보나치 함수의 점화식 {1, 1, 2, 3, 5, 8, 13, 21} 을 넣었다면, 결과물은 {1, 1} 이다. 파도반 수열의 경우에는 $P_n = P_{n-2} + P_{n-3}$ 과 같은 동치의 점화식이 존재하기 때문에, 이 경우 {0, 1, 1} 이 반환된다.

조심해야 할 것은 두 가지다. 첫 번째로, 초기항이 충분히 많아야 한다는 점이다. 만약에 점화식이 $DP_n = \sum_{i = 1}^{k}{DP_{n - i} \times C_i}$ 의 형태라면, 항을 최소한 $3k$ 개보다 많이 넣어줘야 한다. Berlekamp-Massey는 가장 짧은(shortest), 즉 $k$ 를 최소로 하는 선형 점화식을 반환한다. 예를 들어, 만약 항을 $x < 2k$ 개 넣었다면, 항상 $x/2$ 정도 크기의 선형 점화식을 찾을 수 있다 (방정식을 풀어주면 해가 나올 것이니..). 해당 점화식에 대해서 원하는 결과를 얻기 위해서는 최소 $3k$ 개의 항을 제공해야 한다. $2k$ 가 아니라 $3k$ 인지, $3k$ 를 넘어가면 자명하지 않은 선형 점화식이 항상 유일한지는 사실 잘 모르겠다. 꽤 중요한 질문인데 (예를 들어 선형 점화식이 유일하지 않으면, 조건을 만족하는 가장 짧은 점화식이 우리가 원하는 점화식이 아니라 틀릴 수 있으니) 이 부분에 대해서 질문한 사람도 없고 답을 구하지도 못했다. 좋은 의견이 있다면 댓글로 알려주세요.

두 번째로는, 모듈러가 소수여야 한다는 점이다. 만약에 BigInteger 형태로 문제를 해결하고 싶다면, 적당한 소수로 해싱한 후 Berlekamp-Massey를 사용하면 될 것이다 (물론 이 경우 get_nth 함수는 쓸모가 없을 것이다). 모듈러가 소수가 아니라면 더 복잡한 Reeds-Sloane 알고리즘을 사용해야 한다.

사실 여기까지 들어보면 상당히 쓸데 없는 것 같다. 왜냐하면 대부분의 복잡한 문제들은 점화식이 1차원으로 끝나지 않고, 2차원이거나 여러 선형 점화식의 연립으로 표현되기 때문이다. 하지만…

일반적으로 동적 계획법이 저렇게 깔끔한 선형 점화식 하나로 표현되는 경우는 많지 않다. 예를 들어서, 2차원 점화식의 형태이거나 ($dp_{i, j} = dp_{i-1, j} + dp_{i-1, j-1}$). 혹은 여러 개의 선형 점화식이 서로간의 dependency를 가지는 식 ($F_i = G_{i-1} - H{i-1}, G_i = H_{i-1} \times 2$..) 이 등장하는 경우들이 있다. 이러한 경우는 위의 경우처럼 선형 점화식으로 모델링하는 것이 불가능해 보이지만, 사실은 가능하다.

이러한 식들을 행렬로 나타내 보자. 예를 들어 선형 점화식인 $DP_{i} = DP_{i-1} + DP_{i-2} + DP_{i-5}$ 는

$ \begin{bmatrix} DP_i \\ DP_{i-1} \\ DP_{i-2} \\ DP_{i-3} \\ DP_{i-4} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} DP_{i-1} \\ DP_{i-2} \\ DP_{i-3} \\ DP_{i-4} \\ DP_{i-5} \end{bmatrix}$

와 같은 행렬로 나타낼 수 있다. 그 외에도, $DP_{i, j} = 2\times DP_{i-1, j} + 5 \times DP_{i-1, j-1}$ 과 같은 식은 ($0 \le j \le 4$ 라 했을 때)

$ \begin{bmatrix} DP_{i, 0} \\ DP_{i, 1} \\ DP_{i, 2} \\ DP_{i, 3} \\ DP_{i, 4} \end{bmatrix} = \begin{bmatrix} 2 & 0 & 0 & 0 & 0 \\ 5 & 2 & 0 & 0 & 0 \\ 0 & 5 & 2 & 0 & 0 \\ 0 & 0 & 5 & 2 & 0 \\ 0 & 0 & 0 & 5 & 2 \end{bmatrix} \times \begin{bmatrix} DP_{i-1, 0} \\ DP_{i-1, 1} \\ DP_{i-1, 2} \\ DP_{i-1, 3} \\ DP_{i-1, 4} \end{bmatrix}$

와 같은 행렬로 나타낼 수 있고. $F_i = 2 \times G_{i-1}, G_i = F_{i-1} + 3 \times H_{i-1}, H_i = F_{i-1} + G_{i-1}$ 와 같은 식은,

$\begin{bmatrix} F_i \\ G_i \\ H_i \end{bmatrix} = \begin{bmatrix} 0 & 2 & 0 \\ 1 & 0 & 3 \\ 1 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} F_{i-1} \\ G_{i-1} \\ H_{i-1} \end{bmatrix}$

로 나타낼 수 있다.

중요한 예시들을 대부분 적었으나 이 예시가 전부는 아니다. 행렬식으로 바꿀 수 없는 점화식이 뭔지, 직관적으로 행렬식으로 바꿀 수 없어 보이지만 사실 가능한 것들은 뭔지에 대해서 생각해 보는 것은 좋을 것이다. 예를 들어서, $F_i = F_{i-1} + G_{i-1} + 10, G_i = 2 \times F_{i} - G_{i-1} - 6$ 과 같은 식은 행렬 식으로 바꿀 수 없어 보이지만.

$\begin{bmatrix} F_i \\ G_i \\ 1 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 10 \\ 2 & -1 & -6 \\ 0 & 0 & 1 \end{bmatrix} \times \begin{bmatrix} F_{i-1} \\ G_{i-1} \\ 1 \end{bmatrix}$

과 같은 형태로 행렬식으로 쓸 수 있으며, 한편 $C_i = \sum_{0 \le j \le n-1}{C_j \times C_{i-j-1}}$ 과 같은 카탈란 점화식은 이러한 트릭을 쓸 수 없다.

이제 중요한 결론을 제시한다: DP 상태 전이가 행렬로 표현 가능하다면, 해당 DP에 대한 선형 점화식이 존재한다. 정확히는, 위 DP 식에서의 $DP_i, DP_{i ,0}, F_i$ (각각) 에 대해서 선형 점화식이 존재하며, 이 점화식의 크기는 행렬의 크기가 $n \times n$ 일때 최대 $n$이다. 고로 3n개의 DP 값을 직접 계산한 후, Berlekamp-Massey를 사용하면 $O(n^2)$ 시간에 점화식을 찾을 수 있으며, 이 점화식의 $k$ 번째 항은 $O(n^2 \log k)$ 에 계산할 수 있다.

이러한 방식의 장점은 다음과 같다:

• 점화식이 저렇게 생겼을 거 같으면 찍어보면 된다. 어떠한 경우에는, $DP_{i, state}$ 에 대한 점화식을 찾았는데 정확한 상태 전이를 찾기 귀찮거나, 아예 그것을 찾는 것이 문제일 수 있다. 어떠한 경우에는, 점화식을 찾지 못했으나 점화식의 형태가 위와 같을 거라는 직관만이 있을 수도 있다. 어떠한 경우에는 아는 게 아무것도 없고, 그냥 이것이 유일한 희망일 수도 있다. 이러한 경우에, 단순히 Berlekamp-Massey를 로컬에서 돌려본 후 작은 점화식이 있는지를 확인해 볼 수 있다. 만약에 충분히 큰 $n$에 대해서 $\frac{n}{2}$ 보다 유의미하게 작은 점화식을 꾸준히 찾을 수 있다면, 문제를 거의 다 풀었다는 좋은 청신호이다.
• 시간 복잡도가 빠르다. 선형 점화식 형태의 전이는 위의 키타마사법이라는 방법으로 빠르게 계산하는 방법이 잘 알려져 있으나, 행렬 형태의 전이는 보통 행렬 곱셈을 동반하여 이러한 최적화가 불가능하다. Berlekamp-Massey는 후자의 전이를 전자로 변환해 주는 알고리즘이 된다.
• 커팅을 동반한다. Berlekamp-Massey는 최소 크기의 점화식을 찾는다. 고로 꼭 행렬의 크기가 $n \times n$ 이라고 점화식의 크기가 $n$ 이 나오는 것이 아니고 더 작은 점화식이 나올 수 있다. 예를 들어서 $DP_{i, state}$ 와 같은 점화식을 찾았는데 상태 전이가 작고 (sparse matrix) 특정 상태가 도달 불가능한 상황이 있을 수 있다. Berlekamp-Massey를 돌리면 그러한 경우를 우리가 따져 줄 필요가 없이, 알아서 최적의 점화식을 찾아준다. 이는 매우 중요한 요소이다. 어떠한 문제의 경우에는 이러한 커팅을 찾는 것이 (요컨데, 방문할 수 없는 상태를 찾아 시간 복잡도나 상수를 유의미하게 줄이는 것이) 문제의 어려운 부분일 수 있다. Berlekamp-Massey 를 사용하면, 이러한 문제를 내가 해결할 필요가 없이 알고리즘단에서 이러한 커팅을 최적으로 사용할 수 있다.

Berlekamp-Massey 알고리즘을 사용한 다양한 연습 문제와 풀이, 몇가지 사실의 증명, 그리고 get_min_poly와 det 함수에 대한 이야기는 구사과 블로그 에서 계속됩니다.