[알고리즘] Berlekamp-Massey 알고리즘 톺아보기

Berlekamp-Massey 알고리즘은 내가 가장 좋아하는 알고리즘 중 하나이다. 지금까지 이 알고리즘을 블랙박스 형태로 내부 동작을 모르는 채 사용했지만 이번 기회에 한 번 알아보는 것도 좋을 것 같아서 찾아보게 되었다. (몰라도 쓰는데 문제는 없다.) Shift-Register Synthesis and BCH Decoding 이라는 1969년에 나온 오래되고 클래식한 논문에 있는 알고리즘이었다. 정말 정말 좋은 리소스가 있어서 쉽게 접할 수 있었다. 엄밀한 증명은 하지 않고 어떻게 동작하는지만 남긴다.

 

선형 점화식

DP에도 많이 나오는 선형 점화식은 1차원 수열상에서 이전 몇 개의 항들의 선형 결합을 통해 다음 항을 만드는 수열이 된다.

$a_{i} = \displaystyle\sum_{j=1}^n c_j a_{i-j}$

위와 같은 형식을 선형 점화식이라고 한다. 수식이 무시무시해보이지만 가장 널리 알려진 피보나치 수열도 위와 같은 선형 점화식으로 나타낼 수 있다.

$F_1 = 1, F_2 = 1, F_{i} = F_{i-1} + F_{i-2} = \displaystyle\sum_{j=1}^2 F_{i-j} (i \ge 3)$

위 선형 점화식을 통해 5개의 항을 계산하면 다음과 같다.

$F_3 = 2, F_4 = 3, F_5 = 5, F_6 = 8, F_7 = 13$

하지만 아래와 같은 점화식은 선형 점화식이 아니다.

$a_i = \displaystyle\sum_{j=1}^n a_{j} a_{n-j+1} $

$a$에 대한 선형성을 가지지 않기 때문이다. 위 수열은 카탈란 수의 점화식이며 이는 선형적이지 않다. 이런 경우 Berlekamp-Massey 알고리즘을 사용할 수 없다.

또한 선형점화식은 다음과 같은 형태일 수도 있다.

$a_i = p + \displaystyle\sum_{j=1}^n c_j a_{i-j}$

여기서 $p$는 상수이다. 이런 경우도 선형 점화식이며 non-homogeneous 선형 점화식이라고 불린다. $p = 0$이면 homogeneous 선형 점화식이다. non-homogeneous 인 경우에도 Berlekamp-Massey 알고리즘을 사용할 수 있다. non-homogeneous 형식은 homogeneous 형태로 바꿀 수도 있다.

$a_i = p + \displaystyle\sum_{j=1}^n c_j a_{i-j}$

$a_{i-1} = p + \displaystyle\sum_{j=1}^n c_j a_{i-1-j}$

$a_i - a_{i-1} = \displaystyle\sum_{j=1}^n c_j (a_{i-j} - a_{i-1-j})$

$d_i = a_i - a_{i-1}$

$d_i$ 는 homogeneous 선형 점화식이 된다.

Berlekamp-Massey 알고리즘

Berlekamp-Massey 알고리즘은 수열의 몇 항이 주어졌을 때 해당 수열의 최소 선형점화식을 찾아주는 신기한 알고리즘이다. 예를 들어 Berlekamp-Massey의 input으로 피보나치 수열 $1,1,2,3,5,8,13,21$을 넣으면 $1,1$을 반환해준다. 이는 피보나치 수열의 점화식의 계수에 해당한다. 

여기서 중요한 점은 최소 선형점화식이라는 것이다. 최소가 아니라면 어떤 점화식이든 올 수 있기 때문이다. 대표적으로 아래와 같은 사진의 경우가 될 수 있다.

 

위 사진은 4차 다항식으로 구해서 위와 같은 결과가 나왔지만 Berlekamp-Massey 알고리즘은 1차 다항식인 $f(x) = 2x-1$을 구해줄 것이다. 여기에서 Berlekamp-Massey 알고리즘으로 다항식을 구할 수 있다는 것도 흥미롭다.

사전 정의

Berlekamp-Massey알고리즘의 동작을 알아보기 위해 몇 가지만 정의하고 넘어가자.

• 계수 표현식이라는 것을 도입하자. 계수 표현식은 선형 점화식을 $c = (c_1, c_2, ..., c_n)$으로 표기한다. (원본 블로그 글과 표기가 다르다!!) 피보나치 선형 점화식은 $c = \{1, 1\}$이다.
• 선형 점화식의 길이는 $n = |c|$ 로 나타낸다. 피보나치 선형 점화식의 길이는 2이다.
• $c(i) = \displaystyle\sum_{i=1}^n c_i s_{i-j}$를 인덱스 $i$에서 계산한다고 하며 $c(i) = s_i$이면 맞았다라고 한다.
• 선형 점화식은 선형성을 가진다. 즉, 두 선형 점화식 $c = (c_1, c_2, ..., c_n), d = (d_1, d_2, ..., d_n)$이 있을 때 $k (c + d) = (k(c_1 + d_1), k(c_2 + d_2), ..., k(c_n + d_n))$이다. $k$는 상수이다.
• 빈 선형 점화식 $c = ()$이 있다면 모든 $i$에서 $c(i) = 0$이다.

알고리즘 Flow

기본적인 아이디어는 다음과 같다.

현재 가지고 있는 $c$로 $i < t$에서 $c(i) = s_i$이지만 $c(t) \neq s_t$라면 $d$를 잡아서 $i < t$일 때 $d(i) = 0$, $d(t) + c(t) = s_t$로 만들어서 $c$를 $d+c$로 대체하는 형식이다.

 

예를 들어 수열 $s = \{ 1, 2, 4, 8, 13, 20, 28, 215, 757, 2186 \}$가 있다고 하자. 이는 계수가 아닌 값을 나타낸 것이다.

초기 $c = ()$에서 시작한다.

$i = 1 ... |s|$까지 진행하며 $c$에 대해 계산한다.

 

$i = 1$

$s_1 = 1$이고 $c$는 비어있기 때문에 정의상 $c(1) = 0$이다. $s_1 \neq c(1)$이므로 $c$의 교체가 필요하다. $c$를 초깃값으로 설정하는 것은 아무거나 해도 된다. $c = (0)$이어도 되고 $c = (12374)$여도 된다. 어차피 이후에 틀리면 다시 대체할 것이니 적당히 1로 초기화하자.

 

$i = 2$

$ s_2 = 2 $이고 $c = \{1\}$이므로 계산하면 $c(2) = c_1 s_1 = 1$이다. $s_2 \neq c(2)$이므로 $c$의 교체가 필요하다. 지금은 처음 교체하는 것이기 때문에 $c = (2)$로 교체하면 된다. 이 때, 현재 교체했다는 사실을 기억한다.

기억해야 할 것은 $i = 2$랑 교체하기 전 $c$인 $(1)$이다.

 

$i = 3$

$ s_3 = 4$이고 $c(3) = c_1 s_2 = 4$이므로 맞았다. 맞았으므로 통과한다.

 

$i = 4$

$ s_4 = 8$이고 $c(4) = c_1 s_3 = 8$이므로 맞았다. 맞았으므로 통과한다. 

 

$i = 5$

$ s_5 = 13$이고 $c(5) = c_1 s_4 = 16$이므로 틀렸다. $s_5 \neq c(5)$이므로 $c$의 교체가 필요하다.

기본적인 아이디어를 되돌아보면 현재 가지고 있던 $c$로 $i<5$에서 $c(i) = s_i$이지만 $c(5) \neq s_5$이므로 $d$를 하나 잡아서$i<5$일 때 $d(i) = 0$, $d(5) + c(5) = s_5$를 만족하게 하여 $c$를 $c+d$로 대체해야 한다.

 

우리는 $d(5) = s_5 - c(5) = -3$이 되고 싶다. 이 값을 $\Delta = -3$이라고 하자.

위 조건을 만족하는 $ d $를 찾는 법은 다음과 같다.

1. $i = 2$에서 교체하기 전 $c = (1)$을 $d$로 둔다. -> $d = (1)$
2. 선택한 $d$의 원소들에 $-1$을 곱한다. -> $d = (-1)$
3. $d$의 왼쪽에 $1$을 추가한다. -> $d = (1, -1)$
4. 직전에 교체한 인덱스를 $f$로 둔다 -> $f = 2$
5. $d$의 모든 원소에 $\frac{\Delta}{d(f)}$를 곱한다. -> $d(f) = d(2) = d_1 s_1 = 1$, $d = (-3, 3)$
6. $d$의 왼쪽에 $i-f-1 = 5-2-1 = 2$개의 $0$을 추가한다. -> $d = (0, 0, -3, 3)$

이제 $d$는 위 조건을 만족한다. 현재 $c = \{2\}$였으므로 $c + d = (2, 0, -3, 3)$이 되고 이것을 $c$로 대체한다. 그렇게 하면 $c = (2, 0, -3, 3)$이 된다.

그리고 교체했다는 사실을 기억한다. $i = 5$, 교체하기 전 $c$인 $(2)$

 

이제 이 과정을 $s$의 모든 원소를 돌 때까지 진행하면 최소 선형 점화식을 구할 수 있다. 이 과정을 쉽게 말하면 $s = \{1, 2, 4, 8, 13, 20\}$을 구하기 위해 $a = \{1, 2, 4, 8, 16, 32\}$를 구하고 $b = \{0, 0, 0, 0, -3, -12\}$인 $b$를 구하여 최종적으로 $s = a + b$하여 구하는 형식이다.

 

여기까지만 알아도 코드 구현하는데 문제는 없다.

 

핵심은 $d$를 구하는 과정인 것 같다. 저 과정으로 어떻게 조건에 맞는 $d$를 구할 수 있는 것일까?

 

$d$ 구하기

위의 6단계에서 단계별로 분석하자.

1단계를 보면 이전에 교체한 인덱스를 $f$라고 두면 해당 위치에서 교체하기 전 점화식인 $c = (c_1, c_2, ..., c_n)$을 $d$로 뒀다. 그렇다면 $d$는 $ i < f$에서 $c(i) = d(i) = s_i$를 만족한다.

 

3단계까지 진행하면 $d = (1, -c_1, -c_2, ..., -c_n)$이 된다. 이 점화식을 계산하면

$d(i+1) = s_{i} + \displaystyle\sum_{j=1}^{n} (-c_j)s_{i-j} = s_{i} - c(i)$

1단계에서 $i < f$에서 $c(i) = s_i$라고 했으므로 $i < f + 1$에서 $d(i) = 0$, $i \ge f + 1$에서 $d(i) = s_{i-1} - c(i-1)$이 되고 $d(f+1) \neq 0$이 된다. $d(f + 1) = 0$이라면 $s_{f} = c(f)$가 되고 이는 $f$의 정의에 따르면 모순이다.

 

그렇다면 5단계에서 곱한 $\frac{\Delta}{d(f)}$는 무엇을 의미할까? 이는 $i = 2$에서 $c = (1)$을 $c = (2)$로 바꾼것과 유사하다. 

 

TBD

 

코드

위 설명을 코드로 옮기면 된다.

mod = 10**9+7


def berlekamp_massey(x):
    if len(x) <= 1:
        return []
    a, b = x[:2]
    x = [i % mod for i in x]
    f, cur, d = 1, [1], [0]
    if a != b:
        cur, d = [b * pow(a, mod - 2, mod)], [1]

    def get(c, ix):
        res = 0
        for i in range(len(c)):
            res = (res + c[i] * x[ix - i]) % mod
        return res

    for i in range(2, len(x)):
        t = get(cur, i - 1)
        if t == x[i]:
            continue
        delta = (x[i] - t) % mod
        d = [1] + [mod - i for i in d]
        mul = delta * pow(get(d, f), mod - 2, mod) % mod
        d = [0] * (i - f - 1) + [(j * mul) % mod for j in d]
        for j in range(len(cur)):
            d[j] = (d[j] + cur[j]) % mod
        cur, d, f = d, cur, i

    return cur

위 Berlekamp-Massey 의 시간 복잡도는 $O(NK)$이다. 여기서 $N$은 입력으로 넣은 수열의 길이이고 $K$는 해당 수열의 선형 점화식 길이이다. $K$는 출력의 길이인데 시간 복잡도에 들어가는 것이 어색할 수 있다. 그래서 $O(N^2)$으로 봐도 무방하다.

보통 수가 너무 커질 것을 생각해서 mod 로 나눈 값을 계산하므로 위 코드에 해당 내용이 반영되어있다.

 

 

Berlekamp-Massey 알고리즘은 선형 점화식만 구해주고 $N$번째 항을 구해주진 않는다. 선형 점화식을 구했으니 직접 계산하면 $O(NK)$, 행렬곱을 사용하면 $O(K^3 \log N)$, 키타마사를 사용하면 $O(K^2 \log N)$, 키타마사 + FFT 를 사용하면 $O(K \log K \log N)$에 구할 수 있다. ($K$는 선형 점화식의 길이) 보통 $N$의 범위를 $10^{18}$로 주기 때문에 엄청난 효율을 보여준다. 하지만 이는 Berlekamp-Massey에 포함된 내용은 아니므로 여기서 다루지 않는다.