Berlekamp-Massey & Kitamasa

Intuition, Vocabulary and Foundations

A recurrence is a compact predictor. Fibonacci satisfies $F_n=F_{n-1}+F_{n-2}$, so two coefficients replace an arbitrarily long table. Berlekamp–Massey reverses that idea: from a sample prefix it infers the shortest linear recurrence consistent with every observed term. Kitamasa then uses the recurrence to jump to a distant index.

For $s_n=c_1s_{n-1}+\cdots+c_Ls_{n-L}$, the characteristic relation is $X^L=c_1X^{L-1}+\cdots+c_L$. Reducing $X^n$ by this relation expresses $s_n$ as a linear combination of the first $L$ terms. Berlekamp–Massey needs a field because each nonzero discrepancy must be invertible; this is automatic modulo a prime but false in general modulo a composite.

Berlekamp–Massey compresses a trusted prefix into its shortest linear predictor, and Kitamasa reduces $X^n$ modulo that predictor's characteristic polynomial. Field division and independent sample validation are correctness requirements, not optional optimizations.

Concept	Exact example	Non-example or boundary	Question answered
Linear complexity	Fibonacci has complexity $2$ with coefficients $(1,1)$	A list of terms alone is not a recurrence proof	What is the shortest matching rule?
Connection polynomial	$1-X-X^2$ annihilates Fibonacci shifts	$1-X$ fails on the prefix	How are prediction errors represented?
Discrepancy	$d_n=s_n+\sum_{i=1}^L C_i s_{n-i}$	Zero discrepancy on a short prefix proves nothing about unseen data	Does the current rule predict sample $n$?
Characteristic reduction	$X^2\equiv X+1$ for Fibonacci	Ordinary numeric division is unrelated	How is a distant term reduced to initial terms?
Field requirement	Every nonzero value modulo prime $p$ is invertible	$2$ has no inverse modulo $6$	When is every update defined?
Sample sufficiency	At least $2L$ exact samples normally identify order $L$	A short accidental fit may fail later	How much evidence supports inference?

Step-by-Step Derivations and Algorithms

Berlekamp–Massey maintains a connection polynomial $C$, its length $L$, and a saved polynomial $B$ from the last length increase. At index $n$, compute

$$d=s_n+\sum_{i=1}^{L}C_i s_{n-i}.$$

If $d=0$, the current rule already predicts the sample. Otherwise subtract a shifted copy of $B$ scaled by $d/b$, where $b$ is the saved nonzero discrepancy. This cancels the current error. When $2L\le n$, minimality forces the new feasible length $n+1-L$ and the old $C$ becomes the backup.

Kitamasa computes $X^n$ in the quotient by $X^L-c_1X^{L-1}-\cdots-c_L$. Binary exponentiation repeatedly multiplies coefficient vectors, then reduces degrees from $2L-2$ down to $L$. If $X^n\equiv\sum_{i=0}^{L-1}a_iX^i$, then $s_n=\sum a_i s_i$.

Berlekamp-Massey & Kitamasa algorithm map

flowchart TD
    A[Normalize field-valued samples] --> B[Compute discrepancies]
    B --> C{Discrepancy zero?}
    C -->|Yes| D[Keep connection polynomial]
    C -->|No| E[Apply scaled shifted correction]
    D --> F[Validate inferred recurrence]
    E --> F
    F --> G[Exponentiate X modulo characteristic polynomial]
    G --> H[Combine with initial terms]

After processing index $n$, the connection polynomial annihilates the observed prefix and has minimal feasible length; every nonzero update cancels the current discrepancy. Classical inference costs $O(NL)$, at worst $O(N^2)$, while Kitamasa evaluates index $n$ in $O(L^2\log n)$ time and $O(L)$ working space.

Implementation, Complexity and Engineering

The implementation requires a prime modulus, normalizes samples, validates the inferred rule against the complete prefix, handles the all-zero sequence, and rejects invalid indices. A prefix fit is only evidence about that prefix; when samples come from measurements or an unknown program, reserve a disjoint suffix for validation.

from math import isqrt


def is_prime(n):
    if not isinstance(n, int) or isinstance(n, bool) or n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    d = 3
    while d <= isqrt(n):
        if n % d == 0:
            return False
        d += 2
    return True


def validate_recurrence(sequence, recurrence, modulus):
    order = len(recurrence)
    if order == 0:
        if any(value % modulus for value in sequence):
            raise ValueError("empty recurrence requires zero samples")
        return
    for index in range(order, len(sequence)):
        predicted = sum(
            recurrence[offset - 1] * sequence[index - offset]
            for offset in range(1, order + 1)
        ) % modulus
        if predicted != sequence[index] % modulus:
            raise ValueError("recurrence does not reproduce samples")


def berlekamp_massey(samples, modulus):
    if not is_prime(modulus):
        raise ValueError("modulus must be prime")
    sequence = [value % modulus for value in samples]
    if not sequence:
        return []
    connection, backup = [1], [1]
    length, shift, saved = 0, 1, 1
    for index in range(len(sequence)):
        discrepancy = sequence[index]
        for offset in range(1, length + 1):
            discrepancy += connection[offset] * sequence[index - offset]
        discrepancy %= modulus
        if discrepancy == 0:
            shift += 1
            continue
        previous = connection[:]
        scale = discrepancy * pow(saved, -1, modulus) % modulus
        needed = len(backup) + shift
        connection.extend([0] * max(0, needed - len(connection)))
        for offset, coefficient in enumerate(backup):
            position = offset + shift
            connection[position] = (
                connection[position] - scale * coefficient
            ) % modulus
        if 2 * length <= index:
            length = index + 1 - length
            backup, saved, shift = previous, discrepancy, 1
        else:
            shift += 1
    recurrence = [(-connection[i]) % modulus for i in range(1, length + 1)]
    validate_recurrence(sequence, recurrence, modulus)
    return recurrence


def combine(left, right, recurrence, modulus):
    order = len(recurrence)
    product = [0] * (2 * order - 1)
    for i, a_value in enumerate(left):
        for j, b_value in enumerate(right):
            product[i + j] = (product[i + j] + a_value * b_value) % modulus
    for degree in range(2 * order - 2, order - 1, -1):
        for offset, coefficient in enumerate(recurrence, start=1):
            product[degree - offset] = (
                product[degree - offset] + product[degree] * coefficient
            ) % modulus
    return product[:order]


def kitamasa(initial, recurrence, index, modulus):
    if not is_prime(modulus):
        raise ValueError("modulus must be prime")
    if not isinstance(index, int) or isinstance(index, bool) or index < 0:
        raise ValueError("index must be nonnegative")
    recurrence = [value % modulus for value in recurrence]
    order = len(recurrence)
    if order == 0:
        if any(value % modulus for value in initial):
            raise ValueError("empty recurrence requires zero initial values")
        return 0
    if len(initial) < order:
        raise ValueError("not enough initial values")
    initial = [value % modulus for value in initial[:order]]
    if index < order:
        return initial[index]
    result = [1] + [0] * (order - 1)
    base = [recurrence[0]] if order == 1 else [0, 1] + [0] * (order - 2)
    exponent = index
    while exponent:
        if exponent & 1:
            result = combine(result, base, recurrence, modulus)
        base = combine(base, base, recurrence, modulus)
        exponent >>= 1
    return sum(a * b for a, b in zip(result, initial)) % modulus


MOD = 1_000_000_007
fib = [0, 1]
for _ in range(18):
    fib.append((fib[-1] + fib[-2]) % MOD)
recurrence = berlekamp_massey(fib, MOD)
assert recurrence == [1, 1]
for index, expected in enumerate(fib):
    assert kitamasa(fib, recurrence, index, MOD) == expected
assert kitamasa(fib, recurrence, 50, MOD) == 586268941

geometric = [pow(3, index, 101) for index in range(20)]
assert berlekamp_massey(geometric, 101) == [3]
assert kitamasa(geometric, [3], 1000, 101) == pow(3, 1000, 101)
assert berlekamp_massey([0] * 12, 101) == []
assert kitamasa([0], [], 10**6, 101) == 0

for operation in (
    lambda: berlekamp_massey([1, 2, 3], 8),
    lambda: kitamasa([0, 1], [1, 1], -1, MOD),
):
    try:
        operation()
    except ValueError:
        pass
    else:
        raise AssertionError("invalid input was accepted")

print("all tests passed")

Method	Cost	Failure mode	Mitigation
Direct iteration	$O(Ln)$ time, $O(L)$ memory	Huge indices are impractical	Use Kitamasa
Berlekamp–Massey	$O(NL)$ time, $O(L)$ state	Too few or corrupted samples yield a false rule	Use at least $2L$ trusted samples and held-out validation
Classical Kitamasa	$O(L^2\log n)$ time, $O(L)$ state	Quadratic multiplication dominates for large $L$	Switch to fast polynomial arithmetic when justified
Matrix power	$O(L^3\log n)$ classically	Ignores scalar recurrence structure	Prefer Kitamasa unless the full state is needed
Prime-field arithmetic	All nonzero discrepancies invert	Composite moduli contain nonunits	Reject rings or use a ring-aware method
Cached model	Amortizes inference across queries	Wrong modulus or initial vector may be reused	Store recurrence, modulus, and initial state together

Applications, Case Studies and Exercises

From twenty Fibonacci samples modulo $10^9+7$, Berlekamp–Massey returns $(1,1)$; Kitamasa then evaluates $F_{50}\equiv586268941\pmod{10^9+7}$ without generating the intervening terms. In a production pipeline, infer from exact training samples, validate on a disjoint suffix, and only then answer indices as large as $10^{18}$.

Practice

Apply Berlekamp-Massey & Kitamasa

Case study 1. Generate twenty Fibonacci terms modulo $101$, trace every nonzero discrepancy, and recover $(1,1)$. Evaluate indices $20$, $50$, and $1000$ and compare with fast doubling. Success means every result agrees.

Case study 2. Hide an order-$8$ recurrence modulo $1{,}000{,}000{,}007$. Generate $32$ training and $16$ validation samples, recover the rule, validate the held-out suffix, and answer $100$ indices up to $10^{18}$. Success means every answer matches the original generator.

Exercises. Define discrepancy and linear complexity; prove that the shifted correction cancels the current error; add batched queries; demonstrate a noninvertible discrepancy modulo a composite; and benchmark an NTT-backed polynomial multiplier against the quadratic version.

Use iteration for nearby indices, matrix powers when full state transitions matter, and Berlekamp–Massey plus Kitamasa when trustworthy field-valued samples reveal a short scalar recurrence. Never extrapolate an unvalidated prefix or perform field division in a composite ring; these boundaries lead directly to fast polynomial and generating-function techniques.

Cookie Consent