Recurrence & Convolution Laboratory

Intuition, Vocabulary and Foundations

A recurrence laboratory treats a sequence as both observations and algebra. From 0, 1, 1, 2, 3, 5, 8, Berlekamp–Massey can propose $(s_n=s_{n-1}+s_{n-2})$; held-out samples then test whether that proposal predicts rather than merely fits. The same rule has a companion matrix, characteristic polynomial $(x^2-x-1)$, and rational generating function.

Convolution is polynomial multiplication: [1,2,3] convolved with [4,5] is [4,13,22,15]. Kitamasa also multiplies polynomials and reduces them modulo a characteristic polynomial, so faster convolution can accelerate high-order recurrence evaluation. Formal power series connect both views through coefficient identities without requiring numerical convergence.

Every representation must describe the same sequence: samples, recurrence coefficients, matrix state, characteristic polynomial, and generating function. An inferred rule remains a hypothesis until it reproduces training and held-out terms over the stated field.

Concept	Example	Boundary or non-example	Question answered
Samples	`0,1,1,2,3,5,8`	A short prefix does not prove one infinite law	What observations constrain the model?
Linear recurrence	$(s_n=s_{n-1}+s_{n-2})$	$(s_n=s_{n-1}^2+1)$ is nonlinear	Do fixed coefficients generate later terms?
Minimal model	Fibonacci has order $(2)$	Padding with zero coefficients is not minimal	What is the shortest consistent state?
Kitamasa reduction	$(x^2\equiv x+1)$	Reducing by an unrelated polynomial	How does $(x^n)$ become a degree-$(\lt d)$ vector?
Matrix model	`[[1,1],[1,0]]`	Mismatched state ordering	How is one recurrence step linearized?
Convolution	`[1,2,3]*[4,5]` gives `[4,13,22,15]`	Pairwise multiplication only	What are product coefficients?
FFT	Complex roots plus rounding	Unbounded coefficients without error analysis	Can approximate transforms recover exact integers?
NTT	Exact arithmetic modulo $(998244353)$	Length not dividing (p-1)	Does the field contain the required roots?
Formal series	$(1/(1-x-x^2))$	Treating $(x)$ as a required numeric value	Which coefficient identities hold to a precision?

Step-by-Step Derivations and Algorithms

Berlekamp–Massey maintains a connection polynomial annihilating the processed prefix. A discrepancy of zero preserves it; a nonzero discrepancy shifts and scales an earlier correction using a field inverse. It returns a shortest recurrence for $(m)$ samples in $(O(m^2))$ field operations.

For recurrence coefficients $(c_1,\ldots,c_d)$, reduce $(x^n)$ modulo $(Q(x)=x^d-c_1x^{d-1}-\cdots-c_d)$. The remainder has degree below $(d)$, and its coefficients form weights for the initial terms. Binary polynomial exponentiation gives Kitamasa complexity $(O(d^2\log n))$ with naive multiplication. Companion-matrix exponentiation gives $(O(d^3\log n))$, usually slower for one scalar but valuable as an independent oracle. Naive convolution is quadratic; FFT and NTT reduce large products to $(O(L\log L))$, with rounding constraints for FFT and root-of-unity constraints for NTT.

Recurrence & Convolution Laboratory algorithm map

flowchart TD
    A[Collect samples and domain] --> B[Infer with Berlekamp Massey]
    B --> C{Held-out validation passes?}
    C -->|No| D[Reject or gather samples]
    C -->|Yes| E{Evaluation task}
    E -->|Reference or coupled state| F[Matrix exponentiation]
    E -->|Distant scalar term| G[Kitamasa reduction]
    G --> H{Polynomial size}
    H -->|Small| I[Naive convolution]
    H -->|Large friendly prime| J[NTT convolution]
    F --> K[Cross-check results]
    I --> K
    J --> K

Berlekamp–Massey is minimal only over a field and for the supplied prefix; held-out validation sets the empirical boundary. Kitamasa and matrices preserve the same recurrence state, while naive convolution costs $(O(d^2))$ and NTT costs $(O(L\log L))$; practical crossover points must be measured.

Implementation, Complexity and Engineering

The executable laboratory below works modulo $(998244353)$. It infers recurrences, evaluates them with Kitamasa, multiplies with naive or iterative NTT backends, and exposes the selection threshold. Deterministic randomized tests compare independent implementations; timing is reported rather than asserted because hardware and runtimes change crossover points.

from random import Random
from time import perf_counter

MOD = 998244353
ROOT = 3


def berlekamp_massey(s, mod=MOD):
    c, b = [1], [1]
    length, shift, last = 0, 1, 1
    for n in range(len(s)):
        discrepancy = s[n] % mod
        for i in range(1, length + 1):
            discrepancy = (discrepancy + c[i] * s[n - i]) % mod
        if discrepancy == 0:
            shift += 1
            continue
        old = c[:]
        scale = discrepancy * pow(last, mod - 2, mod) % mod
        if len(c) < len(b) + shift:
            c += [0] * (len(b) + shift - len(c))
        for i, value in enumerate(b):
            c[i + shift] = (c[i + shift] - scale * value) % mod
        if 2 * length <= n:
            length = n + 1 - length
            b, last, shift = old, discrepancy, 1
        else:
            shift += 1
    return [(-c[i]) % mod for i in range(1, length + 1)]


def combine(a, b, rec, mod=MOD):
    d = len(rec)
    product = [0] * (2 * d - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            product[i + j] = (product[i + j] + x * y) % mod
    for degree in range(2 * d - 2, d - 1, -1):
        for offset, coefficient in enumerate(rec, 1):
            product[degree - offset] = (product[degree - offset] + product[degree] * coefficient) % mod
    return product[:d]


def kitamasa(initial, rec, index, mod=MOD):
    if index < 0 or not rec or len(initial) < len(rec):
        raise ValueError("invalid recurrence state or index")
    if index < len(rec):
        return initial[index] % mod
    d = len(rec)
    result = [1] + [0] * (d - 1)
    base = [rec[0] % mod] if d == 1 else [0, 1] + [0] * (d - 2)
    while index:
        if index & 1:
            result = combine(result, base, rec, mod)
        base = combine(base, base, rec, mod)
        index >>= 1
    return sum(result[i] * initial[i] for i in range(d)) % mod


def naive_convolution(a, b, mod=MOD):
    if not a or not b:
        return []
    result = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            result[i + j] = (result[i + j] + x * y) % mod
    return result


def ntt(a, invert, mod=MOD):
    n = len(a)
    if not n or n & (n - 1) or (mod - 1) % n:
        raise ValueError("unsupported NTT length")
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j & bit:
            j ^= bit
            bit >>= 1
        j ^= bit
        if i < j:
            a[i], a[j] = a[j], a[i]
    length = 2
    while length <= n:
        root = pow(ROOT, (mod - 1) // length, mod)
        if invert:
            root = pow(root, mod - 2, mod)
        for start in range(0, n, length):
            factor = 1
            for offset in range(length // 2):
                even = a[start + offset]
                odd = a[start + offset + length // 2] * factor % mod
                a[start + offset] = (even + odd) % mod
                a[start + offset + length // 2] = (even - odd) % mod
                factor = factor * root % mod
        length *= 2
    if invert:
        inverse = pow(n, mod - 2, mod)
        for i in range(n):
            a[i] = a[i] * inverse % mod


def ntt_convolution(a, b, mod=MOD):
    if not a or not b:
        return []
    output = len(a) + len(b) - 1
    size = 1
    while size < output:
        size *= 2
    left = [x % mod for x in a] + [0] * (size - len(a))
    right = [x % mod for x in b] + [0] * (size - len(b))
    ntt(left, False, mod)
    ntt(right, False, mod)
    left = [x * y % mod for x, y in zip(left, right)]
    ntt(left, True, mod)
    return left[:output]


def convolution(a, b, threshold=64):
    if threshold < 1:
        raise ValueError("threshold must be positive")
    if min(len(a), len(b)) <= threshold:
        return naive_convolution(a, b), "naive"
    return ntt_convolution(a, b), "ntt"


fibonacci = [0, 1]
for _ in range(2, 80):
    fibonacci.append((fibonacci[-1] + fibonacci[-2]) % MOD)
recurrence = berlekamp_massey(fibonacci[:30])
assert recurrence == [1, 1]
for index, expected in enumerate(fibonacci):
    assert kitamasa([0, 1], recurrence, index) == expected

rng = Random(11)
for n in range(1, 24):
    for m in range(1, 24):
        a = [rng.randrange(MOD) for _ in range(n)]
        b = [rng.randrange(MOD) for _ in range(m)]
        assert ntt_convolution(a, b) == naive_convolution(a, b)

assert naive_convolution([1, 2, 3], [4, 5]) == [4, 13, 22, 15]
start = perf_counter()
trace = kitamasa([0, 1], recurrence, 1_000_000)
elapsed = perf_counter() - start
print("trace:", recurrence, trace, elapsed)
print("all tests passed")

Method	Complexity and domain	Best use	Failure mode and mitigation
Iteration	$(O(nd))$ time	Nearby terms or streams	Too slow for distant indices; route to exponentiation
Berlekamp–Massey	$(O(m^2))$ over a field	Inferring a short recurrence	Noninvertible discrepancy over a ring; require a field
Kitamasa	$(O(d^2\log n))$ with naive products	Distant scalar terms	Coefficient-order bugs; regenerate samples as an oracle
Companion matrix	$(O(d^3\log n))$	Small coupled states and cross-checking	State layout mismatch; test one-step transitions
Naive convolution	$(O(ab))$	Small inputs and oracle tests	Quadratic growth; expose a threshold
FFT	$(O(L\log L))$ with rounding	Bounded integer convolution	Precision loss; split coefficients or verify bounds
NTT	$(O(L\log L))$ exactly modulo a friendly prime	Large modular products	Unsupported length; verify divisibility by (p-1)
Multi-prime NTT	Several transforms plus CRT	Exact or arbitrary-modulus products	Insufficient modulus product; prove a coefficient bound

Applications, Case Studies and Exercises

An approachable experiment infers Fibonacci's recurrence, validates it on held-out samples, and compares distant terms with an independent matrix implementation. A production laboratory groups sampled sequences by field, preserves training and validation splits, benchmarks polynomial backends across orders, and records machine, seed, modulus, repetitions, warm-up policy, and raw timings.

Exercises should progress from translating a recurrence into a characteristic polynomial, through proving Kitamasa reduction, to adding arbitrary-modulus convolution. Debugging targets include inadequate samples, reversed coefficients, unsupported NTT lengths, composite-modulus inverses, and FFT rounding.

Deliverable

Apply Recurrence & Convolution Laboratory

Case study 1 — recover and predict. Generate 30 Fibonacci samples, train on 20, and validate on 10. Confirm coefficients [1, 1], then compare Kitamasa and matrix results at indices $(50)$, $(1000)$, and $(10^6)$. Success means all held-out and distant terms agree.

Case study 2 — measure crossover. For balanced lengths 8, 16, 32, 64, 128, 256, 512, generate coefficients with seed 7, warm both backends, and collect seven timing repetitions. Verify every coefficient and report the first size where NTT has a lower median. Do not hard-code the result from another machine.

Deliverable. Build infer, evaluate, convolve, and benchmark commands emitting modulus, seed, backend, threshold, input sizes, elapsed time, and verification status. Pass tests for zero sequences, order one, insufficient state, empty polynomials, negative coefficients, invalid thresholds, and unsupported transform lengths.

Use iteration nearby, matrices as a transparent oracle, Kitamasa for distant scalar terms, naive convolution below a measured threshold, and FFT or NTT only inside their numeric domains. A prefix fit is not a proven law, and an asymptotic advantage is not a benchmark result; preserve independent representations for trustworthy cross-checks.

Cookie Consent