"""Part 7: find the rule behind a sequence of numbers.

Given the first terms of a sequence, Berlekamp-Massey finds the shortest linear recurrence
    s[n] = c1 s[n-1] + c2 s[n-2] + ... + cL s[n-L]
that produces all of them. It reads the terms one at a time and keeps the best recurrence so far.
When the next term disagrees with what the recurrence predicts, it corrects the recurrence using
the last one that failed, scaled so the mistake cancels, and only makes the recurrence longer
when it has to. With exact fractions there is no rounding anywhere. Given 2L terms of a sequence
that really does follow a recurrence of length L, the answer is that recurrence.

Elwyn Berlekamp invented it for decoding error-correcting codes (1968), and James Massey saw that
it finds the shortest linear feedback shift register for any sequence (1969).
"""
from fractions import Fraction


def berlekamp_massey(seq):
    """Returns [c1, ..., cL], as Fractions, for the shortest recurrence the terms satisfy."""
    s = [Fraction(x) for x in seq]
    c, b = [Fraction(1)], [Fraction(1)]        # connection polynomials, current and last
    L, m, last = 0, 1, Fraction(1)
    for n in range(len(s)):
        d = s[n] + sum(c[i] * s[n - i] for i in range(1, L + 1))   # how wrong c is here
        if d == 0:
            m += 1
            continue
        coef = d / last
        t = list(c)
        c = c + [Fraction(0)] * (len(b) + m - len(c))
        for i, bi in enumerate(b):
            c[i + m] -= coef * bi
        if 2 * L <= n:
            L, b, last, m = n + 1 - L, t, d, 1
        else:
            m += 1
    return [-x for x in c[1:L + 1]]


def extend(seq, rec, n):
    """The sequence continued to n terms with the recurrence."""
    out = list(seq)
    while len(out) < n:
        out.append(sum(c * out[-1 - i] for i, c in enumerate(rec)))
    return out