Matrix Exponentiation

Why Matrix Exponentiation Matters

A process repeatedly applies the same linear update: $x_{t+1}=Ax_t$. After $p$ steps, its state is $x_p=A^p x_0$. Computing the update one step at a time takes work proportional to $p$, which becomes expensive when only a distant state is needed.

Matrix exponentiation groups those repeated updates by powers of two. Instead of forming every intermediate power, compute $A$, $A^2$, $A^4$, and so on, then combine the powers selected by the binary digits of $p$. The idea supports linear recurrences, fixed-length graph-walk counts, repeated geometric transformations, and discrete-time transition models.

This lesson implements integer matrix powers with plain nested arrays. Nonnegative powers require multiplication and an identity matrix. Negative powers additionally require an invertible matrix; the programs compute its inverse with partial pivoting and a floating-point tolerance, then apply the same repeated-squaring loop. They print numerical checks for both paths. The operation here is the integer power $A^p$, not the matrix exponential $e^A$.

Core Concepts and Terminology

For a square $d\times d$ matrix, $A^p$ means repeated matrix multiplication when $p\gt0$, $A^0=I_d$, and $A^{-p}=(A^{-1})^p$ when $p\gt0$ and the inverse exists. The identity convention also applies to the zero matrix's zeroth power. A rectangular matrix is outside this implementation's contract.

Matrix multiplication is associative: $(AB)C=A(BC)$. It is generally not commutative: $AB$ may differ from $BA$. Exponentiation uses associativity and products of powers of the same base; it does not license rearranging unrelated matrices. Ordinary multiplication computes $(AB)_{ij}=\sum_{k=0}^{d-1}A_{ik}B_{kj}$, rather than raising each entry independently.

The running variables are an accumulated result $R$, a current base $B$, and a remaining nonnegative exponent $q$. For target base $C$ and magnitude $m$, begin with $R=I$, $B=C$, and $q=m$. The invariant is

$$ R B^q=C^m. $$

If the requested exponent is negative, choose $C=A^{-1}$; otherwise choose $C=A$. Halving $q$ consumes one binary digit per iteration. NIST's repeated-squaring reference describes this accumulate-and-square pattern.

Our worked base is

$$ A=\begin{pmatrix}0&1\\1&1\end{pmatrix},\qquad x_0=\begin{pmatrix}0\\1\end{pmatrix}. $$

Multiplication maps a state $(u,v)^{\mathsf T}$ to $(v,u+v)^{\mathsf T}$. Thus $A^p x_0=(F_p,F_{p+1})^{\mathsf T}$ for the Fibonacci sequence $F_0=0$, $F_1=1$. This connects matrix powers to a familiar recurrence while keeping every sample entry small.

How Matrix Exponentiation Works, Step by Step

First validate a nonempty square matrix with finite entries. The exponent must be an integer. Copy the base so the caller's matrix stays unchanged.

1. If the exponent is negative, compute an inverse by Gauss–Jordan elimination on $[A\mid I]$. At each column, select the largest absolute pivot among the active rows, swap the whole augmented rows, normalize the pivot row, and eliminate that column from every other row.
2. During inversion, use scale $\alpha=\max_{ij}\lvert A_{ij}\rvert$ and reject a pivot whose magnitude is at most $\tau\alpha$, with $0\lt\tau\lt1$. This is a numerical acceptance rule, not an exact rank certificate.
3. Set $R=I$, let $B$ be the copied matrix or its computed inverse, and let $q=\lvert p\rvert$.
4. While $q\gt0$, multiply $R\leftarrow RB$ if $q$ is odd.
5. Halve $q$ using integer division. If it remains positive, square $B\leftarrow BB$ for the next iteration.
6. When $q=0$, return $R$. Skipping the final unused square saves one matrix multiplication and avoids generating an unnecessary larger power.

Every multiplication writes into a fresh matrix. Reusing an input matrix's entries as output before all dot products finish would corrupt later terms. For positive powers, pivoting is unnecessary; it belongs only to the inversion path for negative powers.

For the sample exponent $5$, the first odd iteration selects $A$, then prepares $A^2$ and leaves $q=2$. The next even iteration leaves $R$ unchanged, prepares $A^4$, and leaves $q=1$. The final iteration selects $A^4$, obtaining $A^5$.

flowchart TD
    A["Read square matrix A and integer p"] --> B{"Valid finite matrix?"}
    B -- No --> C["Reject input"]
    B -- Yes --> D{"p negative?"}
    D -- Yes --> E["Compute pivoted inverse with relative tolerance"]
    E --> F{"Every pivot accepted?"}
    F -- No --> G["Reject negative power"]
    F -- Yes --> H["Set base to inverse"]
    D -- No --> I["Set base to copied A"]
    H --> J["Set result = identity; q = magnitude of p"]
    I --> J
    J --> K{"q positive?"}
    K -- No --> R["Return result"]
    K -- Yes --> L{"q odd?"}
    L -- Yes --> M["Multiply result by base"]
    L -- No --> N["Halve q"]
    M --> N
    N --> O{"q still positive?"}
    O -- Yes --> P["Square base"]
    P --> K
    O -- No --> R

Worked Example with a Full Trace

Use $A=\begin{pmatrix}0&1\\1&1\end{pmatrix}$ and exponent $p=5$. Initialize $R=I_2$, $B=A$, and $q=5$.

First iteration. The exponent is odd, so $R=IA=A$. Halve $q$ to $2$ and square the base:

$$ B=A^2=\begin{pmatrix}0\cdot0+1\cdot1&0\cdot1+1\cdot1\\1\cdot0+1\cdot1&1\cdot1+1\cdot1\end{pmatrix} =\begin{pmatrix}1&1\\1&2\end{pmatrix}. $$

The state is $R=A$, $B=A^2$, $q=2$. The invariant reads $A(A^2)^2=A^5$.

Second iteration. The exponent is even, so keep $R=A$. Halve $q$ to $1$ and square the base:

$$ B=A^4=\begin{pmatrix}1\cdot1+1\cdot1&1\cdot1+1\cdot2\\1\cdot1+2\cdot1&1\cdot1+2\cdot2\end{pmatrix} =\begin{pmatrix}2&3\\3&5\end{pmatrix}. $$

The state is $R=A$, $B=A^4$, $q=1$.

Third iteration. Select the current base:

$$ R=AA^4=\begin{pmatrix}0\cdot2+1\cdot3&0\cdot3+1\cdot5\\1\cdot2+1\cdot3&1\cdot3+1\cdot5\end{pmatrix} =\begin{pmatrix}3&5\\5&8\end{pmatrix}. $$

Halve $q$ to zero and stop. There is no final square. The algorithm used two accumulator products and two squarings. Applying the result to $x_0=(0,1)^{\mathsf T}$ gives $(5,8)^{\mathsf T}$.

For the negative-power path, inversion starts at

$$ [A\mid I]=\left[\begin{array}{cc|cc}0&1&1&0\\1&1&0&1\end{array}\right]. $$

The first column selects row $1$ and swaps the augmented rows. Both pivots are then $1$, so normalization leaves them unchanged. Eliminating the upper entry of column $1$ gives

$$ [I\mid A^{-1}]=\left[\begin{array}{cc|cc}1&0&-1&1\\0&1&1&0\end{array}\right]. $$

For exponent $-5$, the same binary loop uses base $D=A^{-1}$. Its squares are $D^2=\begin{pmatrix}2&-1\\-1&1\end{pmatrix}$ and $D^4=\begin{pmatrix}5&-3\\-3&2\end{pmatrix}$. The selected product is $D D^4=\begin{pmatrix}-8&5\\5&-3\end{pmatrix}=A^{-5}$. Multiplying $A^5A^{-5}$ gives $I_2$, providing an additional hand check.

Implementation in Python, C++, and Java

Each program implements square-matrix multiplication, pivoted inversion, and signed integer powering using plain nested arrays. Nonnegative powers also work for singular matrices; negative powers require inversion to pass the tolerance test. Exponent zero returns a fresh identity after matrix validation. The tolerance is used when inversion is required.

The portable exponent interface is signed $64$-bit in C++ and Java; Python also supports larger integers. The C++ magnitude calculation avoids negating the minimum signed value. Java treats the magnitude's bit pattern as unsigned and shifts it with >>>, which also handles Long.MIN_VALUE. These details prevent an exponent-boundary bug without changing the mathematical loop.

All programs use double precision. They print A^0 = [[1.000000, 0.000000], [0.000000, 1.000000]], A^5 = [[3.000000, 5.000000], [5.000000, 8.000000]], A^-1 = [[-1.000000, 1.000000], [1.000000, 0.000000]], A^-5 = [[-8.000000, 5.000000], [5.000000, -3.000000]], inverse residual = 0.000e+00, power error = 0.000e+00, and check = PASS, one per line. The inverse residual is the maximum absolute entry of $AA^{-1}-I$; the power error compares $A^5$ with five successive multiplications. Save the Java panel as Main.java.

The examples remain exactly representable, but general powers accumulate roundoff and may overflow or underflow. The code targets moderate, sensibly scaled matrices and does not estimate conditioning. A small inverse residual alone does not ensure accurate large negative powers.

from math import isfinite


def checked_copy(a: list[list[float]]) -> list[list[float]]:
    n = len(a)
    if n == 0 or any(len(row) != n for row in a):
        raise ValueError("require a nonempty square matrix")
    if any(not isfinite(v) for row in a for v in row):
        raise ValueError("matrix entries must be finite")
    return [row[:] for row in a]


def identity(n: int) -> list[list[float]]:
    return [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]


def multiply(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:
    n = len(a)
    out = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            for j in range(n):
                out[i][j] += a[i][k] * b[k][j]
    return out


def inverse(a: list[list[float]], tol: float = 1e-12) -> list[list[float]]:
    a = checked_copy(a)
    if not isfinite(tol) or not 0 < tol < 1:
        raise ValueError("require 0 < tol < 1")
    n = len(a)
    eye = identity(n)
    aug = [a[i] + eye[i] for i in range(n)]
    scale = max(abs(v) for row in a for v in row)
    for k in range(n):
        pivot = max(range(k, n), key=lambda i: abs(aug[i][k]))
        if abs(aug[pivot][k]) <= tol * scale:
            raise ValueError("inverse rejected by pivot tolerance")
        aug[k], aug[pivot] = aug[pivot], aug[k]
        divisor = aug[k][k]
        for j in range(2 * n):
            aug[k][j] /= divisor
        for i in range(n):
            if i != k:
                factor = aug[i][k]
                for j in range(2 * n):
                    aug[i][j] -= factor * aug[k][j]
                aug[i][k] = 0.0
    return [row[n:] for row in aug]


def matrix_power(a: list[list[float]], exponent: int,
                 tol: float = 1e-12) -> list[list[float]]:
    base = checked_copy(a)
    if not isinstance(exponent, int):
        raise ValueError("exponent must be an integer")
    if exponent < 0:
        base = inverse(base, tol)
    remaining = abs(exponent)
    result = identity(len(base))
    while remaining:
        if remaining & 1:
            result = multiply(result, base)
        remaining //= 2
        if remaining:
            base = multiply(base, base)
    return result


def max_error(a: list[list[float]], b: list[list[float]]) -> float:
    return max(abs(a[i][j] - b[i][j]) for i in range(len(a)) for j in range(len(a)))


def show(a: list[list[float]]) -> str:
    return "[" + ", ".join("[" + ", ".join(f"{v:.6f}" for v in row) + "]" for row in a) + "]"


def main() -> None:
    a = [[0.0, 1.0], [1.0, 1.0]]
    powers = [matrix_power(a, p) for p in [0, 5, -1, -5]]
    for p, value in zip([0, 5, -1, -5], powers):
        print(f"A^{p} = {show(value)}")
    naive = identity(2)
    for _ in range(5):
        naive = multiply(naive, a)
    residual = max_error(multiply(a, powers[2]), identity(2))
    error = max_error(powers[1], naive)
    print(f"inverse residual = {residual:.3e}")
    print(f"power error = {error:.3e}")
    print("check = " + ("PASS" if residual <= 1e-10 and error <= 1e-10 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

Matrix checked_copy(const Matrix& a) {
    if (a.empty()) throw std::invalid_argument("matrix must be nonempty");
    for (const auto& row : a) {
        if (row.size() != a.size()) throw std::invalid_argument("matrix must be square");
        for (double v : row) if (!std::isfinite(v)) throw std::invalid_argument("entries must be finite");
    }
    return a;
}

Matrix identity(std::size_t n) {
    Matrix out(n, std::vector<double>(n, 0));
    for (std::size_t i = 0; i < n; ++i) out[i][i] = 1;
    return out;
}

Matrix multiply(const Matrix& a, const Matrix& b) {
    const std::size_t n = a.size();
    Matrix out(n, std::vector<double>(n, 0));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t k = 0; k < n; ++k)
            for (std::size_t j = 0; j < n; ++j) out[i][j] += a[i][k] * b[k][j];
    return out;
}

Matrix inverse(const Matrix& input, double tol = 1e-12) {
    Matrix a = checked_copy(input);
    if (!std::isfinite(tol) || tol <= 0 || tol >= 1) throw std::invalid_argument("invalid tolerance");
    const std::size_t n = a.size();
    Matrix aug(n, std::vector<double>(2 * n, 0));
    double scale = 0;
    for (std::size_t i = 0; i < n; ++i) {
        for (std::size_t j = 0; j < n; ++j) {
            aug[i][j] = a[i][j];
            scale = std::max(scale, std::abs(a[i][j]));
        }
        aug[i][n + i] = 1;
    }
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        for (std::size_t i = k + 1; i < n; ++i)
            if (std::abs(aug[i][k]) > std::abs(aug[pivot][k])) pivot = i;
        if (std::abs(aug[pivot][k]) <= tol * scale) throw std::runtime_error("inverse rejected by pivot tolerance");
        std::swap(aug[k], aug[pivot]);
        const double divisor = aug[k][k];
        for (std::size_t j = 0; j < 2 * n; ++j) aug[k][j] /= divisor;
        for (std::size_t i = 0; i < n; ++i) if (i != k) {
            const double factor = aug[i][k];
            for (std::size_t j = 0; j < 2 * n; ++j) aug[i][j] -= factor * aug[k][j];
            aug[i][k] = 0;
        }
    }
    Matrix out(n, std::vector<double>(n));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) out[i][j] = aug[i][n + j];
    return out;
}

Matrix matrix_power(const Matrix& a, long long exponent, double tol = 1e-12) {
    Matrix base = checked_copy(a);
    unsigned long long remaining;
    if (exponent < 0) {
        base = inverse(base, tol);
        remaining = static_cast<unsigned long long>(-(exponent + 1)) + 1ULL;
    } else {
        remaining = static_cast<unsigned long long>(exponent);
    }
    Matrix result = identity(base.size());
    while (remaining != 0) {
        if ((remaining & 1ULL) != 0) result = multiply(result, base);
        remaining /= 2;
        if (remaining != 0) base = multiply(base, base);
    }
    return result;
}

double max_error(const Matrix& a, const Matrix& b) {
    double result = 0;
    for (std::size_t i = 0; i < a.size(); ++i)
        for (std::size_t j = 0; j < a.size(); ++j) result = std::max(result, std::abs(a[i][j] - b[i][j]));
    return result;
}

void show(long long exponent, const Matrix& a) {
    std::cout << "A^" << exponent << " = [" << std::fixed << std::setprecision(6);
    for (std::size_t i = 0; i < a.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << '[';
        for (std::size_t j = 0; j < a.size(); ++j) {
            if (j != 0) std::cout << ", ";
            std::cout << a[i][j];
        }
        std::cout << ']';
    }
    std::cout << "]\n";
}

int main() {
    const Matrix a = { {0, 1}, {1, 1}};
    const Matrix zero = matrix_power(a, 0), positive = matrix_power(a, 5);
    const Matrix inv = matrix_power(a, -1), negative = matrix_power(a, -5);
    show(0, zero); show(5, positive); show(-1, inv); show(-5, negative);
    Matrix naive = identity(2);
    for (int i = 0; i < 5; ++i) naive = multiply(naive, a);
    const double residual = max_error(multiply(a, inv), identity(2));
    const double error = max_error(positive, naive);
    std::cout << std::scientific << std::setprecision(3)
              << "inverse residual = " << residual << '\n' << "power error = " << error << '\n';
    std::cout << "check = " << (residual <= 1e-10 && error <= 1e-10 ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.Locale;

public class Main {
    private static double[][] checkedCopy(double[][] a) {
        if (a.length == 0) throw new IllegalArgumentException("matrix must be nonempty");
        double[][] out = new double[a.length][];
        for (int i = 0; i < a.length; ++i) {
            if (a[i].length != a.length) throw new IllegalArgumentException("matrix must be square");
            out[i] = a[i].clone();
            for (double v : a[i]) if (!Double.isFinite(v)) throw new IllegalArgumentException("entries must be finite");
        }
        return out;
    }

    public static double[][] identity(int n) {
        double[][] out = new double[n][n];
        for (int i = 0; i < n; ++i) out[i][i] = 1;
        return out;
    }

    public static double[][] multiply(double[][] a, double[][] b) {
        int n = a.length;
        double[][] out = new double[n][n];
        for (int i = 0; i < n; ++i)
            for (int k = 0; k < n; ++k)
                for (int j = 0; j < n; ++j) out[i][j] += a[i][k] * b[k][j];
        return out;
    }

    public static double[][] inverse(double[][] input, double tol) {
        double[][] a = checkedCopy(input);
        if (!Double.isFinite(tol) || tol <= 0 || tol >= 1) throw new IllegalArgumentException("invalid tolerance");
        int n = a.length;
        double[][] aug = new double[n][2 * n];
        double scale = 0;
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                aug[i][j] = a[i][j];
                scale = Math.max(scale, Math.abs(a[i][j]));
            }
            aug[i][n + i] = 1;
        }
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            for (int i = k + 1; i < n; ++i)
                if (Math.abs(aug[i][k]) > Math.abs(aug[pivot][k])) pivot = i;
            if (Math.abs(aug[pivot][k]) <= tol * scale) throw new IllegalArgumentException("inverse rejected by pivot tolerance");
            double[] row = aug[k]; aug[k] = aug[pivot]; aug[pivot] = row;
            double divisor = aug[k][k];
            for (int j = 0; j < 2 * n; ++j) aug[k][j] /= divisor;
            for (int i = 0; i < n; ++i) if (i != k) {
                double factor = aug[i][k];
                for (int j = 0; j < 2 * n; ++j) aug[i][j] -= factor * aug[k][j];
                aug[i][k] = 0;
            }
        }
        double[][] out = new double[n][n];
        for (int i = 0; i < n; ++i) System.arraycopy(aug[i], n, out[i], 0, n);
        return out;
    }

    public static double[][] matrixPower(double[][] a, long exponent, double tol) {
        double[][] base = checkedCopy(a);
        if (exponent < 0) base = inverse(base, tol);
        long remaining = exponent < 0 ? -exponent : exponent;
        double[][] result = identity(base.length);
        while (remaining != 0) {
            if ((remaining & 1L) != 0) result = multiply(result, base);
            remaining >>>= 1;
            if (remaining != 0) base = multiply(base, base);
        }
        return result;
    }

    public static double maxError(double[][] a, double[][] b) {
        double result = 0;
        for (int i = 0; i < a.length; ++i)
            for (int j = 0; j < a.length; ++j) result = Math.max(result, Math.abs(a[i][j] - b[i][j]));
        return result;
    }

    private static void show(long exponent, double[][] a) {
        System.out.print("A^" + exponent + " = [");
        for (int i = 0; i < a.length; ++i) {
            if (i != 0) System.out.print(", ");
            System.out.print("[");
            for (int j = 0; j < a.length; ++j) {
                if (j != 0) System.out.print(", ");
                System.out.printf(Locale.ROOT, "%.6f", a[i][j]);
            }
            System.out.print("]");
        }
        System.out.println("]");
    }

    public static void main(String[] args) {
        double[][] a = { {0, 1}, {1, 1}};
        double[][] zero = matrixPower(a, 0, 1e-12), positive = matrixPower(a, 5, 1e-12);
        double[][] inv = matrixPower(a, -1, 1e-12), negative = matrixPower(a, -5, 1e-12);
        show(0, zero); show(5, positive); show(-1, inv); show(-5, negative);
        double[][] naive = identity(2);
        for (int i = 0; i < 5; ++i) naive = multiply(naive, a);
        double residual = maxError(multiply(a, inv), identity(2));
        double error = maxError(positive, naive);
        System.out.printf(Locale.ROOT, "inverse residual = %.3e%npower error = %.3e%n", residual, error);
        System.out.println("check = " + (residual <= 1e-10 && error <= 1e-10 ? "PASS" : "FAIL"));
    }
}

Complexity, Correctness, and Invariants

With ordinary multiplication, each dense $d\times d$ product costs $\Theta(d^3)$. For magnitude $m=\lvert p\rvert\ge1$, let $\ell=\lfloor\log_2 m\rfloor+1$ and let $s$ be the number of set bits in $m$. The shown loop performs exactly $\ell-1$ squarings and $s$ accumulator products. Thus nonzero positive powers take $\Theta(d^3\log(m+1))$ time. Negative powers add a $\Theta(d^3)$ inversion step, preserving that overall bound on successful inputs.

Exponent zero takes $\Theta(d^2)$ time to validate, copy, and return the identity. For fixed dimension and a fixed positive magnitude, matrix values do not alter the multiplication loop counts, so best, average, and worst cases have the same asymptotic bound. Among magnitudes of a fixed bit length, the set-bit count changes the number of accumulator products but not the logarithmic order. Negative inputs may be rejected early by inversion; an average including failures requires an input distribution.

The code retains only a constant number of dense matrices at a time, giving $\Theta(d^2)$ live workspace. Recursive exponentiation would also use recursion stack space. Floating-point arithmetic costs are treated as constant per operation; exact big-integer entries introduce bit costs that grow with the powers.

For the loop proof, assume $RB^q=C^m$. If $q=2t$, keeping $R$ and squaring $B$ gives $R(B^2)^t=C^m$. If $q=2t+1$, multiplying the accumulator by $B$ gives $(RB)(B^2)^t=C^m$. Associativity suffices in both cases. When $q=0$, the invariant yields $R=C^m$. This proves the nonnegative powering loop in exact arithmetic. For a negative exponent, row operations on $[A\mid I]$ first establish the inverse when all required pivots are nonzero; numerical tolerance and rounding qualify the floating-point result.

Variants, Pitfalls, and Optimizations

Matrix power versus elementwise power. For the sample, $A^2=\begin{pmatrix}1&1\\1&2\end{pmatrix}$, while squaring its entries individually gives $\begin{pmatrix}0&1\\1&1\end{pmatrix}$. They represent different operations. Likewise, the matrix exponential $e^A=\sum_{k\ge0}A^k/k!$ is a separate matrix function with different algorithms.

Zero and negative exponents. Initialize the accumulator to the identity, not the zero matrix. A singular matrix still has nonnegative powers, but no negative integer power. For negative powers, a global pivot threshold can reject an invertible but poorly scaled matrix; it is not a condition-number estimate or exact singularity classification.

Exact recurrence and walk counts. When entries are counts, use exact integers or modular multiplication rather than floating-point arithmetic if exact answers are required. An adjacency matrix with $A_{ij}$ equal to the number of edges from $i$ to $j$ has $(A^p)_{ij}$ equal to the number of length-$p$ walks. These counts can grow quickly. Choose a type or modulus that safely handles intermediate products, not just final outputs.

State orientation and affine updates. The lesson uses column states $x_{t+1}=Ax_t$. A row-state formulation multiplies on the other side. To model $x_{t+1}=Ax_t+b$, append a constant coordinate and power the block matrix $\begin{pmatrix}A&b\\0&1\end{pmatrix}$. Time-varying transitions require their ordered product, not a power of one fixed matrix.

Sparsity and repeated queries. Squaring a sparse matrix can create many new nonzeros; sparse input alone does not imply sparse powers. If only $A^p x$ is needed, applying selected powers directly to a vector can avoid accumulating a full result matrix, although forming the squared bases can still dominate. Cache powers of two when many queries share the same base and account for the extra storage.

Floating-point growth. Repeated squaring reduces the number of products, but it does not guarantee less numerical error for every matrix. Powers can magnify perturbations, negative powers can amplify inversion error, and large entries can overflow. Compare suitably scaled errors for general applications. The small exact samples provide checks of the implementation's indexing, bit logic, and pivoting rather than a universal accuracy guarantee.

Practice Problems and Real-World Use

Choose the state vector, define one update, and verify that the same matrix applies at every step. Derive a tiny case by hand before selecting an exponentiation routine. The examples below separate algebraic correctness, exact counting, and numerical behavior.

Verify positive powers against sequential multiplication on small exponents, compare negative powers with independently solved inverse systems, and include zero, one, singular matrices, and row-swap cases. Check that inputs remain unchanged. Exponent-boundary tests on identity or permutation matrices exercise the integer loop without introducing numerical growth.