Matrix Exponentiation

Modular exponentiation computed $a^n$ in $O(\log n)$ multiplications by reading the bits of $n$ and repeatedly squaring. That lesson noted that the doubling idea is not special to integers: replace multiply with matrix multiply and the same loop computes the $n$-th Fibonacci number in $O(\log n)$ steps. This lesson makes that precise and general. A linear recurrence — a rule that builds each term from a fixed linear combination of the previous few — is repeated application of one fixed matrix. Computing the $n$-th term is then computing the $n$-th power of that matrix, and exponentiation by squaring does it in logarithmic time.

A recurrence is a matrix

Start with Fibonacci: $F_0=0$, $F_1=1$, and $F_n=F_{n-1}+F_{n-2}$. The naive loop is $\Theta (n)$ additions, and the closed form (Binet's formula) involves irrational powers of the golden ratio, awkward to evaluate exactly. Instead, package two consecutive terms into a state vector and ask what one step does to it. The update $F_{n+1}=F_n+F_{n-1}$ together with the trivial $F_n=F_n$ is a linear map:

Call that $2\times 2$ matrix $M$. One multiplication by $M$ advances the state by one term. Two multiplications advance it by two. Applying $M$ a total of $n$ times to the initial state $\left[\begin{array}{c}{F}_1\\ F_0\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right]$ lands on $\left[\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right]$. So the whole question reduces to computing the matrix power $M^n$.

The top-right entry of $M^n$ is $F_n$ directly, so once we can raise $M$ to the $n$-th power we can read off $F_n$ with no vector multiply at all. It remains to compute $M^n$ fast.

Exponentiation by squaring, over matrices

Matrix multiplication is associative, which is the only property repeated squaring ever used. So the entire argument from the modular-exponentiation lesson carries over verbatim with the scalar product replaced by the matrix product and the scalar $1$ replaced by the identity matrix $I$. Write $n$ in binary as $n={\sum }_i{b}_i{2}^i$; then

where the powers $M^{2^i}$ are the repeated squares of $M$, each the square of the previous. Sweeping the bits of $n$ from low to high, we keep a running square and fold it into an accumulator whenever the bit is set.

There are $\lfloor {\log }_{2}n\rfloor +1$ squarings and at most that many accumulator multiplications, so $O(\log n)$ matrix multiplications in all. Each multiplication of two $k\times k$ matrices is $\Theta (k^3)$ by the schoolbook method, giving the headline cost.

For Fibonacci $k=2$, so this is a flat $O(\log n)$ — exponentially faster than the $\Theta (n)$ loop once $n$ is large, and exact (only integer arithmetic, no floating-point golden ratio). The schoolbook $\Theta (k^3)$ per product can be shaved to $\Theta (k^{2.807})$ with Strassen's algorithm, but for the small $k$ typical of recurrences the constant-factor simplicity of the cubic method wins.

from __future__ import annotations

from typing import Optional, Sequence

class Matrix:
  """
    A dense rectangular matrix of integers, stored row by row.\n
    The companion-matrix and walk-counting algorithms in this lesson all\n
    multiply and power square matrices through this one class.\n
  """

  def __init__(self, rows: Sequence[Sequence[int]]) -> None:
    """
      Build a matrix from a sequence of equal-length integer rows.\n
    """
    # copy the rows and record the shape.
    self.rows: list[list[int]] = [list(row) for row in rows]
    self.row_count: int = len(self.rows)
    self.column_count: int = len(self.rows[0]) if self.rows else 0

    # reject ragged input.
    for row in self.rows:
      if len(row) != self.column_count:
        raise ValueError("every row must have the same length")

  @classmethod
  def identity(cls, size: int) -> Matrix:
    """
      The size x size identity matrix — the unit for matrix multiply.\n
    """
    return cls([
      [1 if column == row else 0 for column in range(size)]
      for row in range(size)
    ])

  def __getitem__(self, position: tuple[int, int]) -> int:
    """
      The entry at (row, column).\n
    """
    row, column = position
    return self.rows[row][column]

  def __eq__(self, other: object) -> bool:
    if not isinstance(other, Matrix):
      return NotImplemented
    return self.rows == other.rows

  def __repr__(self) -> str:
    return f"Matrix({self.rows!r})"

  def multiply(self, other: Matrix, modulus: Optional[int] = None) -> Matrix:
    """
      The schoolbook product of this matrix with `other`, optionally\n
      reduced entrywise mod `modulus`. Costs O(rows * inner * columns).\n
    """
    if self.column_count != other.row_count:
      raise ValueError("inner dimensions must match for multiplication")

    # each output entry is a dot product of a left row with a right column.
    product: list[list[int]] = []
    for left_row in self.rows:
      result_row: list[int] = []
      for column in range(other.column_count):

        # dot the row against the column, reducing mod if asked.
        total: int = sum(
          left_row[inner] * other.rows[inner][column]
          for inner in range(self.column_count)
        )
        result_row.append(total % modulus if modulus is not None else total)

      product.append(result_row)
    return Matrix(product)

  def power(self, exponent: int, modulus: Optional[int] = None) -> Matrix:
    """
      This (square) matrix raised to a non-negative `exponent` by repeated\n
      squaring — the same low-to-high bit scan as scalar mod_pow, with the\n
      identity matrix as the accumulator's starting value. O(k^3 log n).\n
    """
    if self.row_count != self.column_count:
      raise ValueError("only square matrices can be raised to a power")
    if exponent < 0:
      raise ValueError("exponent must be non-negative")

    # accumulate from the identity; square repeatedly while scanning the bits.
    result: Matrix = Matrix.identity(self.row_count)
    square: Matrix = self
    remaining: int = exponent

    while remaining > 0:

      # fold the running square into the result when the low bit is set.
      if remaining & 1:
        result = result.multiply(square, modulus)
      square = square.multiply(square, modulus)
      remaining >>= 1

    # the seed identity is unreduced, so reduce it once at the end.
    if modulus is not None:
      result = Matrix([[entry % modulus for entry in row] for row in result.rows])
    return result

  def apply(
    self,
    vector: Sequence[int],
    modulus: Optional[int] = None,
  ) -> list[int]:
    """
      This matrix times a column `vector`, optionally reduced mod `modulus`.\n
    """
    if len(vector) != self.column_count:
      raise ValueError("vector length must match the column count")

    # each output entry dots one row against the vector.
    result: list[int] = []
    for row in self.rows:
      total: int = sum(entry * value for entry, value in zip(row, vector))
      result.append(total % modulus if modulus is not None else total)
    return result

The loop is byte-for-byte the iterative $\text{ModPow}$ from the previous lesson, with the scalar identity $1$ swapped for the matrix identity $I$ and $\cdot$ now meaning matrix multiply. To compute a recurrence modulo a prime — the usual demand when the true terms overflow — reduce every entry mod $p$ after each multiply, exactly as $\text{ModPow}$ reduces after each scalar product.

Worked example ($F_{10}$ by squaring). Take $M=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$ and compute $M^{10}$; the top-right entry is $F_{10}=55$. The exponent is $10=1010_2$, so the set bits are at positions $1$ and $3$. Build the repeated squares:

Each is the square of the previous — for instance $M^4=(M^2{)}^2$ has top-left $2\cdot 2+1\cdot 1=5$. The bits select $M^2$ and $M^8$, so

and reading the entries confirms $F_{11}=89$, $F_{10}=55$, $F_9=34$ — the Fibonacci identity $M^n=\left[\begin{array}{cc}{F}_{n+1}& F_n\\ F_n& F_{n-1}\end{array}\right]$ in a single squaring chain of four multiplies instead of ten additions. The gap only widens with $n$: reaching $F_{10^{18}}\mathrm{mod}p$ is $\approx 60$ multiplies, not $10^{18}$ additions.

The general $k$-term recurrence

Nothing above used $k=2$. Suppose a sequence obeys an order-$k$ linear recurrence with fixed coefficients $c_1,\dots ,c_k$,

Carry the last $k$ terms as the state vector $v_n=(a_n,a_{n-1},\dots ,a_{n-k+1}{)}^{\top }$. The new top entry is the recurrence itself; every other entry is just a copy shifted down by one. That compute the combination on top, shift the rest down rule is a single matrix, the companion matrix of the recurrence:

The top row holds the coefficients and produces $a_n$; the sub-diagonal of $1$s copies each old term down one slot, discarding the oldest. The figure makes the two roles visible: one row that mixes, and a staircase of $1$s that shifts.

Then $v_n=C^{n-k+1}{v}_{k-1}$ where $v_{k-1}$ holds the $k$ given initial terms, and $\text{Mat-Pow}$ delivers $C$'s power in $O(k^3\log n)$. The pattern generalizes: counting problems whose answer satisfies a fixed linear recurrence — tilings, constrained strings, walk counts — are all handled this way. Climbing Stairs is literally Fibonacci; Count Vowels Permutation and Student Attendance Record II track a constant number of states whose transitions are linear, so each is a small companion (or transfer) matrix raised to the $n$-th power.¹

from typing import Optional, Sequence

from matrix_power import Matrix

def companion_matrix(coefficients: Sequence[int]) -> Matrix:
  """
    The companion matrix of an order-k recurrence whose coefficients are\n
    (c1, c2, ..., ck) applied to (a(n-1), ..., a(n-k)). The top row holds\n
    the coefficients; the sub-diagonal of 1s copies each old term down one\n
    slot, discarding the oldest.\n
  """
  order: int = len(coefficients)
  if order == 0:
    raise ValueError("a recurrence needs at least one coefficient")

  # top row is the coefficients; each later row is a shifting 1.
  rows: list[list[int]] = [list(coefficients)]
  for shifted_row in range(1, order):
    row: list[int] = [0 for _ in range(order)]
    row[shifted_row - 1] = 1
    rows.append(row)

  return Matrix(rows)

def linear_recurrence(
  coefficients: Sequence[int],
  initial_terms: Sequence[int],
  index: int,
  modulus: Optional[int] = None,
) -> int:
  """
    The term a(index) of the order-k recurrence\n
    a(n) = c1*a(n-1) + ... + ck*a(n-k), given the k `initial_terms`\n
    a(0), a(1), ..., a(k-1) and the matching `coefficients` (c1, ..., ck).\n
    Computed in O(k^3 log index) matrix multiplications; pass `modulus` to\n
    reduce every entry after each product. The index must be non-negative.\n
  """
  order: int = len(coefficients)
  if len(initial_terms) != order:
    raise ValueError("need exactly one initial term per coefficient")
  if index < 0:
    raise ValueError("index must be non-negative")

  # the first k terms are given outright; no powering needed.
  if index < order:
    value: int = initial_terms[index]
    return value % modulus if modulus is not None else value

  # state holds (a(k-1), a(k-2), ..., a(0)) top-to-bottom; powering the
  # companion matrix (index - (k-1)) times advances it to (a(index), ...).
  state: list[int] = list(reversed(initial_terms))
  transition: Matrix = companion_matrix(coefficients)

  # advance the state to a(index) by powering the companion matrix.
  powered: Matrix = transition.power(index - (order - 1), modulus)
  return powered.apply(state, modulus)[0]

def affine_recurrence(
  coefficients: Sequence[int],
  constant: int,
  initial_terms: Sequence[int],
  index: int,
  modulus: Optional[int] = None,
) -> int:
  """
    The term a(index) of the inhomogeneous recurrence\n
    a(n) = c1*a(n-1) + ... + ck*a(n-k) + `constant`. The trick is to append\n
    a constant 1-track to the state and a row that carries `constant`, so\n
    the augmented (k+1) x (k+1) matrix is again fixed and the same power\n
    method applies. Initial terms and index are as in `linear_recurrence`.\n
  """
  order: int = len(coefficients)
  if len(initial_terms) != order:
    raise ValueError("need exactly one initial term per coefficient")
  if index < 0:
    raise ValueError("index must be non-negative")

  if index < order:
    value: int = initial_terms[index]
    return value % modulus if modulus is not None else value

  # augment the order-k state with a trailing constant 1; the extra column
  # in the top row injects `constant` into each new term, and a final
  # identity row keeps the constant track fixed at 1.
  # top row gains the `constant` column; older rows gain a trailing 0.
  base: Matrix = companion_matrix(coefficients)
  augmented_rows: list[list[int]] = [[*base.rows[0], constant]]
  augmented_rows.extend([*row, 0] for row in base.rows[1:])

  # a final identity row pins the constant track at 1.
  augmented_rows.append([0 for _ in range(order)] + [1])
  transition: Matrix = Matrix(augmented_rows)

  # carry the constant alongside the state, then power as before.
  state: list[int] = list(reversed(initial_terms)) + [1]
  powered: Matrix = transition.power(index - (order - 1), modulus)
  return powered.apply(state, modulus)[0]

Worked example (a $2\times 1$ tiling count). The number of ways to tile a $2\times n$ board with dominoes obeys $T_n=T_{n-1}+T_{n-2}$ with $T_0=T_1=1$ — a shifted Fibonacci. A richer case is the $3$-term tribonacci count $a_n=a_{n-1}+a_{n-2}+a_{n-3}$ with $a_0=0,a_1=a_2=1$, whose companion matrix is

The top row adds the last three terms; the two sub-diagonal $1$s shift the window. Starting from $v_2=(a_2,a_1,a_0{)}^{\top }=(1,1,0{)}^{\top }$, one multiply gives $C{v}_2=(1+1+0,1,1{)}^{\top }=(2,1,1{)}^{\top }=(a_3,a_2,a_1{)}^{\top }$, so $a_3=2$; another gives $a_4=2+1+1=4$, then $a_5=4+2+1=7$. To leap to $a_n$ directly, raise $C$ to the $n-2$ power by $\text{Mat-Pow}$: $a_{60}$ is $58$ matrix multiplies of a fixed $3\times 3$ matrix, independent of how large the terms grow.

Graph walks: the same matrix, a different reading

The companion matrix is one instance of a broader principle: whenever a process advances by a fixed linear rule, its $n$-step behaviour is a matrix power. The cleanest other instance is counting walks in a graph. Let $A$ be the adjacency matrix of a graph, so $A_{ij}=1$ when there is an edge $i\to j$. Then the matrix power counts paths:

Worked example (counting length-$2$ walks). Take the triangle-with-a-tail graph on vertices $\{1,2,3\}$ with edges $1-2$, $2-3$, $1-3$ (an undirected triangle), adjacency matrix

Read $(A^2{)}_{11}=2$: there are exactly two length-$2$ walks from vertex $1$ back to itself, namely $1\to 2\to 1$ and $1\to 3\to 1$. The off-diagonal $(A^2{)}_{12}=1$ counts the single length-$2$ walk $1\to 3\to 2$. Every entry is a sum over the intermediate vertex $\ell$, exactly ${\sum }_{\ell }{A}_{1\ell }{A}_{\ell 2}$, and raising $A$ to higher powers counts longer walks with no new idea — just more multiplies, each delivered in $O(\log n)$ of them by $\text{Mat-Pow}$.

So how many length-$n$ walks? is answered in $O(V^3\log n)$ by raising $A$ to the $n$-th power with $\text{Mat-Pow}$ — and the companion-matrix recurrence above is just this theorem applied to the small state-transition graph of the recurrence. The same $O(V^3\log n)$ power, taken over the boolean semiring (replace $+$ with OR, $\times$ with AND), computes reachability within $n$ steps; over the min-plus semiring (replace $+$ with $\min$, $\times$ with $+$) it computes shortest-path lengths, which underlies the all-pairs shortest-path algorithms.

from __future__ import annotations

from collections.abc import Hashable
from typing import Optional, TypeVar

from graph import Graph
from matrix_power import Matrix

Label = TypeVar("Label", bound=Hashable)

def adjacency_matrix(graph: Graph[Label]) -> tuple[Matrix, list[Label]]:
  """
    The 0/1 adjacency matrix of `graph` and the label order indexing it.\n
    Entry (i, j) is 1 when an edge runs from the i-th to the j-th label.\n
    Parallel edges between the same pair collapse to a single 1.\n
  """
  # fix a label order and an empty size x size grid of zeros.
  ordered_labels: list[Label] = graph.labels
  index_of: dict[Label, int] = {
    label: position for position, label in enumerate(ordered_labels)
  }
  size: int = len(ordered_labels)
  rows: list[list[int]] = [[0 for _ in range(size)] for _ in range(size)]

  # mark a 1 for every edge, collapsing parallel edges.
  for vertex in graph:
    source_index: int = index_of[vertex.label]
    for neighbor in vertex.neighbors():
      rows[source_index][index_of[neighbor.label]] = 1

  return Matrix(rows), ordered_labels

def count_walks(
  graph: Graph[Label],
  source: Label,
  target: Label,
  length: int,
  modulus: Optional[int] = None,
) -> int:
  """
    The number of walks of exactly `length` edges from `source` to `target`\n
    in `graph`, read off (A^length)[source][target]. A length of 0 counts\n
    the empty walk: 1 when source equals target, else 0. O(V^3 log length).\n
  """
  if length < 0:
    raise ValueError("length must be non-negative")

  # build A and recover the label-to-index map.
  matrix, ordered_labels = adjacency_matrix(graph)
  index_of: dict[Label, int] = {
    label: position for position, label in enumerate(ordered_labels)
  }

  # raise A to the length and read off the source-target entry.
  powered: Matrix = matrix.power(length, modulus)
  return powered[index_of[source], index_of[target]]

def total_walks(
  graph: Graph[Label],
  length: int,
  modulus: Optional[int] = None,
) -> int:
  """
    The total number of walks of exactly `length` edges over all ordered\n
    pairs of vertices — the sum of every entry of A^length.\n
  """
  if length < 0:
    raise ValueError("length must be non-negative")

  # power A and sum every entry — walks over all ordered pairs.
  matrix, _ = adjacency_matrix(graph)
  powered: Matrix = matrix.power(length, modulus)
  total: int = sum(entry for row in powered.rows for entry in row)
  return total % modulus if modulus is not None else total

import math
from typing import Sequence

# A walk of "no edges" has cost zero; an absent walk has infinite cost.
INFINITY: float = math.inf

def boolean_multiply(
  left: Sequence[Sequence[bool]],
  right: Sequence[Sequence[bool]],
) -> list[list[bool]]:
  """
    Boolean matrix product: entry (i, j) is True when some index k has both\n
    left[i][k] and right[k][j] — that is, OR over AND.\n
  """
  # OR over AND: (i, j) is True once any k links left[i][k] to right[k][j].
  size: int = len(left)
  product: list[list[bool]] = [[False for _ in range(size)] for _ in range(size)]
  for row in range(size):
    for column in range(size):
      product[row][column] = any(
        left[row][inner] and right[inner][column] for inner in range(size)
      )
  return product

def boolean_power(
  matrix: Sequence[Sequence[bool]],
  exponent: int,
) -> list[list[bool]]:
  """
    The `exponent`-th boolean power of a square matrix by repeated squaring.\n
    With an adjacency matrix, entry (i, j) of the result is True exactly\n
    when j is reachable from i by a walk of length exactly `exponent`. The\n
    identity is the boolean identity (True on the diagonal). O(V^3 log n).\n
  """
  if exponent < 0:
    raise ValueError("exponent must be non-negative")

  # seed with the boolean identity (True on the diagonal) and a copy of A.
  size: int = len(matrix)
  result: list[list[bool]] = [
    [row == column for column in range(size)] for row in range(size)
  ]
  square: list[list[bool]] = [list(row) for row in matrix]

  # square repeatedly, folding in when the low bit is set.
  remaining: int = exponent
  while remaining > 0:
    if remaining & 1:
      result = boolean_multiply(result, square)
    square = boolean_multiply(square, square)
    remaining >>= 1

  return result

def min_plus_multiply(
  left: Sequence[Sequence[float]],
  right: Sequence[Sequence[float]],
) -> list[list[float]]:
  """
    Min-plus (tropical) matrix product: entry (i, j) is the minimum over k\n
    of left[i][k] + right[k][j] — that is, min over +. INFINITY marks the\n
    absence of an edge or path.\n
  """
  # min over +: (i, j) is the cheapest two-hop path left[i][k] + right[k][j].
  size: int = len(left)
  product: list[list[float]] = [[INFINITY for _ in range(size)] for _ in range(size)]
  for row in range(size):
    for column in range(size):
      product[row][column] = min(
        (left[row][inner] + right[inner][column] for inner in range(size)),
        default=INFINITY,
      )
  return product

def min_plus_power(
  matrix: Sequence[Sequence[float]],
  exponent: int,
) -> list[list[float]]:
  """
    The `exponent`-th min-plus power of a weight matrix by repeated\n
    squaring. With a weighted adjacency matrix (INFINITY off the edges,\n
    0 on the diagonal), entry (i, j) is the lightest walk from i to j using\n
    at most `exponent` edges. The identity is 0 on the diagonal and\n
    INFINITY elsewhere — the min-plus unit. O(V^3 log n).\n
  """
  if exponent < 0:
    raise ValueError("exponent must be non-negative")

  # seed with the min-plus unit (0 on the diagonal, INFINITY off) and a copy.
  size: int = len(matrix)
  result: list[list[float]] = [
    [0.0 if row == column else INFINITY for column in range(size)]
    for row in range(size)
  ]
  square: list[list[float]] = [list(row) for row in matrix]

  # square repeatedly, folding in when the low bit is set.
  remaining: int = exponent
  while remaining > 0:
    if remaining & 1:
      result = min_plus_multiply(result, square)
    square = min_plus_multiply(square, square)
    remaining >>= 1

  return result

from collections.abc import Hashable, Iterator
from typing import Generic, Optional, TypeVar


Label = TypeVar("Label", bound=Hashable)


class Edge(Generic[Label]):
  """
    A directed connection from `source` to `target`, carrying a weight.\n
  """

  def __init__(
    self,
    source: Vertex[Label],
    target: Vertex[Label],
    weight: float = 1.0,
  ) -> None:
    self.source: Vertex[Label] = source
    self.target: Vertex[Label] = target
    self.weight: float = weight

  def __repr__(self) -> str:
    return f"Edge({self.source.label!r} -> {self.target.label!r}, w={self.weight})"


class Vertex(Generic[Label]):
  """
    A graph vertex: a label plus the list of edges leaving it.\n
  """

  def __init__(self, label: Label) -> None:
    self.label: Label = label
    self.outgoing: list[Edge[Label]] = []

  def neighbors(self) -> list[Vertex[Label]]:
    """
      The vertices reachable from this one by a single edge.\n
    """
    return [edge.target for edge in self.outgoing]

  def edge_to(self, label: Label) -> Optional[Edge[Label]]:
    """
      The outgoing edge to the vertex with `label`, or None.\n
    """
    for edge in self.outgoing:
      if edge.target.label == label:
        return edge
    return None

  def __repr__(self) -> str:
    return f"Vertex({self.label!r})"


class Graph(Generic[Label]):
  """
    A graph of Vertex objects linked by Edge objects.\n
    Pass `directed=True` for a digraph; otherwise each `add_edge` inserts\n
    the reverse edge too.\n
  """

  def __init__(self, directed: bool = False) -> None:
    self.directed: bool = directed
    self._vertices: dict[Label, Vertex[Label]] = {}

  def add_vertex(self, label: Label) -> Vertex[Label]:
    """
      Return the vertex for `label`, creating it if it is absent.\n
    """
    # reuse the existing vertex, or mint and register a fresh one.
    vertex = self._vertices.get(label)
    if vertex is None:
      vertex = Vertex(label)
      self._vertices[label] = vertex
    return vertex

  def add_edge(
    self,
    source_label: Label,
    target_label: Label,
    weight: float = 1.0,
  ) -> None:
    """
      Connect two labels (creating either vertex as needed).\n
      Adds the reverse edge as well when the graph is undirected.\n
    """
    source = self.add_vertex(source_label)
    target = self.add_vertex(target_label)

    # link source to target, and mirror it back when undirected.
    source.outgoing.append(Edge(source, target, weight))
    if not self.directed:
      target.outgoing.append(Edge(target, source, weight))

  def vertex(self, label: Label) -> Vertex[Label]:
    """
      The vertex carrying `label` (raises KeyError if absent).\n
    """
    return self._vertices[label]

  @property
  def vertices(self) -> list[Vertex[Label]]:
    """
      Every vertex, in insertion order.\n
    """
    return list(self._vertices.values())

  @property
  def labels(self) -> list[Label]:
    """
      Every vertex label, in insertion order.\n
    """
    return list(self._vertices)

  def edges(self) -> Iterator[Edge[Label]]:
    """
      Each edge once — an undirected edge is yielded a single time.\n
    """
    # track undirected endpoint pairs so each is emitted only once.
    seen: set[frozenset[Label]] = set()

    for vertex in self._vertices.values():
      for edge in vertex.outgoing:
        # skip an undirected edge already yielded from the other endpoint.
        if not self.directed:
          endpoints = frozenset((edge.source.label, edge.target.label))
          if endpoints in seen:
            continue
          seen.add(endpoints)

        yield edge

  def __contains__(self, label: Label) -> bool:
    return label in self._vertices

  def __iter__(self) -> Iterator[Vertex[Label]]:
    return iter(self._vertices.values())

  def __len__(self) -> int:
    return len(self._vertices)

Faster recurrence solvers and the semiring view

Matrix powering is the textbook route to a linear recurrence's $n$-th term, but it is neither the only route nor always the fastest.

Faster than cubic: Kitamasa and Fiduccia. The $O(k^3\log n)$ cost is dominated by full $k\times k$ multiplies. But the answer $a_n$ is a fixed linear combination of the first $k$ terms, and the coefficients of that combination are the residue of $x^n$ modulo the recurrence's characteristic polynomial$x^k-c_1{x}^{k-1}-\cdots -c_k$. Computing $x^n\mathrm{mod}$ that polynomial by repeated squaring of polynomials — the Fiduccia / Kitamasa method — costs $O(k^2\log n)$ with schoolbook polynomial multiplication, and $O(k\log k\log n)$ when the polynomial products use the FFT.² For a recurrence of large order $k$ this is a decisive improvement over the cubic matrix approach.

Eigen-decomposition and Binet's formula. If the companion matrix is diagonalizable, $C=PD{P}^{-1}$ with $D$ diagonal, then $C^n=P{D}^n{P}^{-1}$ and $D^n$ is just the eigenvalues raised to the $n$-th power. For Fibonacci the eigenvalues are the golden ratio $\phi =\frac{1+\sqrt{5}}{2}$ and its conjugate, which recovers Binet's formula $F_n=({\phi }^n-{\psi }^n)/\sqrt{5}$ exactly. The closed form is exact on paper but numerically fragile (the irrational powers must be carried to high precision), which is why the integer matrix power, exact by construction, is the algorithm of choice.³

One idea, three algorithms. The semiring view unifies several algorithms in this course. Over the ordinary ring, $A^n$ counts walks; over the boolean semiring $(\vee ,\wedge )$ it computes reachability, and the repeated-squaring version is essentially transitive closure; over the min-plus (tropical) semiring $(\min ,+)$ it computes shortest walks of bounded length, and squaring the distance matrix $\lceil \log n\rceil$ times gives the repeated-squaring all-pairs shortest paths algorithm that precedes Floyd–Warshall. The doubling trick on this page is the same one used in those graph algorithms.

Takeaways

• A linear recurrence of order $k$ is one fixed matrix applied repeatedly: pack the last $k$ terms into a state vector, and one step is a multiply by the companion matrix $C$ (top row = coefficients, sub-diagonal of $1$s = shift).
• Computing the $n$-th term is computing $C^n$, done in $O(\log n)$ matrix multiplications by exponentiation by squaring — the identical bit-scanning loop as $\text{ModPow}$, with scalar $\cdot$ replaced by matrix $\cdot$ and $1$ by $I$.
• Total cost is $O(k^3\log n)$ (cubic per multiply, logarithmically many multiplies); for Fibonacci $k=2$ this is a flat $O(\log n)$, exact and overflow-free under a modulus.
• The same matrix power counts length-$n$ walks in a graph via $(A^n{)}_{ij}$, and over other semirings (boolean, min-plus) computes reachability and shortest paths — one algebraic idea behind recurrences, walk-counting, and transitive closure.