Back to Number Theory Dsa

Linear Recurrences, Matrices & Generating Functions

September 22, 2026 Wasil Zafar 6 min read

Accelerate recurrence evaluation with algebraic representations.

Table of Contents

  1. Intuition, Vocabulary and Foundations
  2. Step-by-Step Derivations and Algorithms
  3. Implementation, Complexity and Engineering
  4. Applications, Case Studies and Exercises

Intuition, Vocabulary and Foundations

A recurrence is a machine with limited memory. Fibonacci begins with $F_0=0$, $F_1=1$ and uses $F_n=F_{n-1}+F_{n-2}$. In general, an order-$k$ homogeneous linear recurrence is

$$a_n=c_0a_{n-1}+c_1a_{n-2}+\cdots+c_{k-1}a_{n-k}.$$

Its contract must specify coefficient order, base values, arithmetic domain, and first valid index. A companion matrix converts one step into a linear transformation; binary exponentiation then jumps to a distant index. The characteristic polynomial

$$P(x)=x^k-c_0x^{k-1}-\cdots-c_{k-1}$$

lets Kitamasa reduce $x^n$ to degree below $k$. An ordinary generating function $A(x)=\sum_{n\ge0}a_nx^n$ turns shifts into multiplication by $x$; Fibonacci has $A(x)=x/(1-x-x^2)$. An exponential generating function uses $a_nx^n/n!$ and naturally models labeled structures.

A recurrence algorithm is correct only relative to an explicit contract: coefficient orientation, base values, valid indices, and arithmetic domain. Matrices, characteristic polynomials, and generating functions must all preserve that same contract.

ConceptDefinitionExampleBoundaryQuestion answered
Linear recurrence$a_n=\sum_{i=1}^kc_{i-1}a_{n-i}$Fibonacci coefficients are $(1,1)$$a_n=a_{n-1}^2+1$ is nonlinearHow does the next term depend on history?
Companion matrixMatrix that shifts the recurrence stateFibonacci uses $\begin{bmatrix}1&1\\1&0\end{bmatrix}$Reversed state reverses coefficientsHow can steps become powers?
Characteristic polynomial$x^k-\sum c_jx^{k-1-j}$$x^2-x-1$Sign conventions must matchHow are high powers reduced?
OGF$A(x)=\sum a_nx^n$$x/(1-x-x^2)$Coefficients are not factorial-scaledHow do shifts become algebra?
EGF$E(x)=\sum a_nx^n/n!$$e^x$ encodes all onesIt is not an OGFHow are labeled objects encoded?
Berlekamp-MasseyShortest recurrence over a fieldInfers Fibonacci from samplesComposite moduli may lack inversesWhat recurrence explains samples?

Step-by-Step Derivations and Algorithms

With state $S_n=[a_n,a_{n-1},\ldots,a_{n-k+1}]^{\mathsf T}$ and companion matrix $M$, we have $S_n=M^{n-k+1}S_{k-1}$. Binary exponentiation gives $O(k^3\log n)$ time using dense multiplication.

Kitamasa computes

$$x^n\bmod P(x)=w_0+w_1x+\cdots+w_{k-1}x^{k-1},$$

then returns $a_n=\sum_iw_ia_i$. Dense polynomial combine-and-reduce costs $O(k^2)$, giving $O(k^2\log n)$ evaluation.

Berlekamp-Massey processes samples left to right. A nonzero discrepancy is canceled by a shifted, scaled correction polynomial. This requires field division and costs $O(N^2)$ for $N$ samples. A rational OGF $P(x)/Q(x)$ induces a recurrence from the denominator after any exceptional prefix.

Linear Recurrences, Matrices & Generating Functions algorithm map
flowchart TD
    A["Start with samples or recurrence"] --> B{"Recurrence known?"}
    B -->|"No"| C["Berlekamp-Massey"]
    C --> D["Validate on held-out terms"]
    B -->|"Yes"| D
    D --> E{"Choose representation"}
    E -->|"Coupled state"| F["Matrix exponentiation"]
    E -->|"One distant term"| G["Kitamasa reduction"]
    E -->|"Symbolic model"| H["Generating-function algebra"]
    H --> G

Matrix exponentiation preserves the state transition and costs $O(k^3\log n)$; Kitamasa preserves $x^n\bmod P(x)$ and costs $O(k^2\log n)$. Berlekamp-Massey requires field inverses and must be validated on samples not used for inference.

Implementation, Complexity and Engineering

The reference supplies an iterative oracle and Kitamasa evaluation under any positive modulus. It treats $n\lt k$ as a base lookup and uses the coefficient convention stated above.

def validate(initial, coefficients, modulus):
    if not initial or len(initial) != len(coefficients):
        raise ValueError("equal positive lengths required")
    if modulus <= 0:
        raise ValueError("positive modulus required")

def iterative_term(initial, coefficients, n, modulus):
    validate(initial, coefficients, modulus)
    if n < 0: raise ValueError("nonnegative index required")
    k = len(initial)
    values = [x % modulus for x in initial]
    if n < k: return values[n]
    for index in range(k, n + 1):
        values.append(sum(coefficients[j] * values[index - 1 - j]
                          for j in range(k)) % modulus)
    return values[n]

def combine(left, right, coefficients, modulus):
    k = len(coefficients)
    product = [0] * (2 * k - 1)
    for i, a in enumerate(left):
        for j, b in enumerate(right):
            product[i + j] = (product[i + j] + a * b) % modulus
    for degree in range(2 * k - 2, k - 1, -1):
        value = product[degree]
        for offset, coefficient in enumerate(coefficients):
            product[degree - 1 - offset] = (
                product[degree - 1 - offset] + value * coefficient
            ) % modulus
    return product[:k]

def kitamasa_term(initial, coefficients, n, modulus):
    validate(initial, coefficients, modulus)
    if n < 0: raise ValueError("nonnegative index required")
    k = len(initial)
    if n < k: return initial[n] % modulus
    result = [1] + [0] * (k - 1)
    base = [coefficients[0] % modulus] if k == 1 else [0, 1] + [0] * (k - 2)
    exponent = n
    while exponent:
        if exponent & 1:
            result = combine(result, base, coefficients, modulus)
        base = combine(base, base, coefficients, modulus)
        exponent >>= 1
    return sum(w * a for w, a in zip(result, initial)) % modulus

if __name__ == "__main__":
    mod = 1_000_000_007
    initial, coefficients = [0, 1], [1, 1]
    expected = [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
    for n, value in enumerate(expected):
        assert iterative_term(initial, coefficients, n, mod) == value
        assert kitamasa_term(initial, coefficients, n, mod) == value
    tri_initial, tri_coefficients = [0, 0, 1], [1, 1, 1]
    for n in range(40):
        assert kitamasa_term(tri_initial, tri_coefficients, n, 10007) == \
               iterative_term(tri_initial, tri_coefficients, n, 10007)
    assert kitamasa_term([7], [3], 4, 1000) == 567
    print("all tests passed")
MethodPreparationQuery costSpaceBest useRisk
Direct iterationNone$O(nk)$$O(k)$ with ring bufferSmall $n$ and oracle testsToo slow for distant indices
Dense matricesCompanion matrix$O(k^3\log n)$$O(k^2)$Coupled state and many outputsState-orientation errors
KitamasaCoefficients$O(k^2\log n)$$O(k)$One distant scalar termWrong reduction convention
Berlekamp-Massey$N$ samples$O(N^2)$ inference$O(N)$Discovery over a fieldNoninvertible discrepancy or overfitting
Rational OGFSymbolic numerator and denominatorExtraction-dependentRepresentation-dependentCounting derivationsMissing prefix corrections

Applications, Case Studies and Exercises

Recurrence acceleration applies to automaton walks, tilings, population models, linear feedback systems, and dynamic programs with fixed transitions. Matrices naturally support coupled state, while Kitamasa is attractive for one scalar term. Inference should train on one prefix and validate on a disjoint suffix.

Practice

Apply Linear Recurrences, Matrices & Generating Functions

Distant Fibonacci: compute $F_{10^{18}}\bmod1{,}000{,}000{,}007$ using $(F_0,F_1)=(0,1)$ and coefficients $(1,1)$. Confirm small terms and compare two logarithmic methods.

Inferred event counts: over the prime field modulo $998{,}244{,}353$, infer a recurrence from a training prefix with Berlekamp-Massey, validate on held-out samples, then evaluate index $10^{12}$ with Kitamasa. Report order, coefficients, and validation length.

Hands-on engine: accept base values, coefficients, modulus, index, and method. Differentially test iteration and Kitamasa for random orders $1$ through $12$. Success requires $10{,}000$ agreeing cases, correct $n\lt k$ behavior, and rejected malformed inputs.

Exercises: build a Tribonacci companion matrix; derive the Fibonacci OGF; prove the binary-power invariant; reduce $x^{10}$ modulo $x^2-x-1$; implement Berlekamp-Massey; demonstrate overfitting from a short prefix; extend a matrix state for affine recurrences.

Use iteration for small indices, matrices for coupled states, Kitamasa for one distant term, and Berlekamp-Massey only over a field with independent validation. Most failures come from ambiguous base indexing, reversed coefficients, or unjustified modular division.

Next in the Series

In Polynomial Algorithms, FFT & NTT, we'll explore Turn polynomial operations into fast convolution primitives..

Continue the Series

Modular Combinatorics

Evaluate combinatorial quantities under prime and composite moduli.

Read Article