diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py
index 268ce7b..a897108 100644
--- a/src/lemke/bimatrix.py
+++ b/src/lemke/bimatrix.py
@@ -1,4 +1,4 @@
-# bimatrix class
+"""Bimatrix game representation and Nash equilibrium computation."""
 
 import fractions
 import random  # random.seed
@@ -11,20 +11,15 @@
 from .randomstart import MAX_ACCURACY
 from .utils import MAXDECIMALS
 
-# file format:
-# <m> <n>
-# m*n entries of A, separated by blanks / newlines
-# m*n entries of B, separated by blanks / newlines
-#
-# blank lines or lines starting with "#" are ignored
-
 # defaults
 # MAXDIM = 2000 # largest allowed value for m and n; not used yet
 
 
-# list generated from string s such as "1-3,10,4-7", all not
-# larger than endrange (50 is arbitrary default) and at least 1
 def rangesplit(s, endrange=50):
+    """Parse a comma-separated range string
+    such as "1-3,10,4-7" into a list of integers.
+    Open-ended ranges like "5-" extend to `endrange`.
+    """
     result = []
     for part in s.split(","):
         if part != "":
@@ -41,9 +36,46 @@ def rangesplit(s, endrange=50):
     return result
 
 
-# used for both A and B
 class payoffmatrix:
-    # create matrix from any numerical matrix
+    """A single player's payoff matrix, stored as exact fractions.
+
+    Wraps a numerical matrix as an array of `fractions.Fraction` entries.
+
+    Parameters
+    ----------
+    A : array_like
+        Numerical payoff matrix of shape ``(m, n)``.
+        Entries may be int, float, or Fraction;
+        they are converted to `Fraction`.
+
+    Attributes
+    ----------
+    numrows : int
+        Number of rows, `m`.
+    numcolumns : int
+        Number of columns, `n`.
+    matrix : numpy.ndarray of fractions.Fraction
+        The payoff matrix as exact fractions.
+    max : fractions.Fraction
+        Running maximum entry of `matrix`.
+    min : fractions.Fraction
+        Running minimum entry of `matrix`.
+    negshift : int
+        ``int(max) + 1``, the shift used to make all entries positive.
+    negmatrix : numpy.ndarray
+        `matrix` shifted by `negshift` so all entries are positive
+        (used when constructing the LCP).
+
+    Examples
+    --------
+    >>> from lemke.bimatrix import payoffmatrix
+    >>> payoffs = payoffmatrix([[3, 0], [0, 2]])
+    >>> payoffs.numrows
+    2
+    >>> payoffs.max
+    3
+    """
+
     def __init__(self, A):
         AA = np.array(A)
         m, n = AA.shape
@@ -66,6 +98,19 @@ def __str__(self):
         return out
 
     def updatemaxmin(self, fromrow, fromcol):
+        """Update `max`, `min`, and `negmatrix` over a submatrix range.
+
+        Recomputes `max`/`min` by scanning entries from
+        ``(fromrow, fromcol)`` to the end of the matrix,
+        then recomputes `negshift` and `negmatrix` from the result.
+
+        Parameters
+        ----------
+        fromrow : int
+            First row index to include in the scan.
+        fromcol : int
+            First column index to include in the scan.
+        """
         m = self.numrows
         n = self.numcolumns
         for i in range(fromrow, m):
@@ -77,31 +122,114 @@ def updatemaxmin(self, fromrow, fromcol):
         self.negmatrix = np.full((m, n), self.negshift, dtype=int) - self.matrix
 
     def fullmaxmin(self):
+        """Recompute `max` and `min` over the entire matrix."""
         self.max = self.matrix[0][0]
         self.min = self.matrix[0][0]
         self.updatemaxmin(0, 0)
 
-    # add full row, row must be of size n
     def addrow(self, row):
+        """Append a row to the matrix and update max/min.
+
+        Parameters
+        ----------
+        row : array_like
+            Row to append; must have length `numcolumns`.
+        """
         self.matrix = np.vstack([self.matrix, row])
         self.numrows += 1
         self.updatemaxmin(self.numrows - 1, 0)
 
-    # add full column, col must be of size m
     def addcolumn(self, col):
+        """Append a column to the matrix and update max/min.
+
+        Parameters
+        ----------
+        col : array_like
+            Column to append; must have length `numrows`.
+        """
         self.matrix = np.column_stack([self.matrix, col])
         self.numcolumns += 1
         self.updatemaxmin(0, self.numcolumns - 1)
 
 
 class bimatrix:
+    """A two-player bimatrix game and methods to find its equilibria.
+
+    Parameters
+    ----------
+    A : payoffmatrix
+        Row player's payoff matrix.
+    B : payoffmatrix
+        Column player's payoff matrix.
+
+    Attributes
+    ----------
+    A : payoffmatrix
+        Row player's payoff matrix.
+    B : payoffmatrix
+        Column player's payoff matrix.
+
+    Examples
+    --------
+    >>> from lemke.bimatrix import bimatrix, payoffmatrix
+    >>> game = bimatrix(
+    ...     A=payoffmatrix([[3, 0], [0, 2]]),
+    ...     B=payoffmatrix([[2, 0], [0, 7]]),
+    ... )
+    >>> game.A.numrows
+    2
+    """
+
     def __init__(self, A, B):
         self.A = A
         self.B = B
 
-    # create A,B from file
     @classmethod
     def from_file(cls, filename):
+        """Create a bimatrix game from a text file.
+
+        Expects the file format::
+
+            <m> <n>
+            m*n entries of A
+            m*n entries of B
+
+        Blank lines and lines starting with ``#`` are ignored.
+
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Path to the game file.
+
+        Returns
+        -------
+        bimatrix
+            The game described by the file.
+
+        Raises
+        ------
+        SystemExit
+            If the number of payoff values in the file is incorrect.
+
+        Examples
+        --------
+        Given a file ``game.txt`` containing::
+            # m,n=
+            2 2
+
+            # A=
+            3 0
+            0 2
+
+            # B=
+            2 0
+            0 7
+
+        >>> from lemke.bimatrix import bimatrix
+        >>> game = bimatrix.from_file("game.txt")
+        >>> print(game.A.numrows)
+        2
+        """
         lines = utils.stripcomments(filename)
         # flatten into words
         words = utils.towords(lines)
@@ -128,6 +256,14 @@ def __str__(self):
         return out
 
     def createLCP(self):
+        """Build the linear complementarity problem for this game.
+
+        Returns
+        -------
+        lcp
+            The LCP corresponding to this game.
+        """
+
         m = self.A.numrows
         n = self.A.numcolumns
         lcpdim = m + n + 2
@@ -155,6 +291,23 @@ def createLCP(self):
         return lemke.lcp(M, q, d)
 
     def runLH(self, droppedlabel):
+        """Run Lemke-Howson from a given dropped label.
+
+        Parameters
+        ----------
+        droppedlabel : int
+            The label to drop as the starting label for the Lemke-Howson path.
+
+        Returns
+        -------
+        tuple of fractions.Fraction
+            The equilibrium strategy pair found, as a tuple of length ``m + n``.
+
+        Raises
+        ------
+        RuntimeError
+            If the underlying Lemke's algorithm unexpectedly fails to find a solution.
+        """
         lcp = self.createLCP()
         lcp.d[droppedlabel - 1] = 0  # subsidize this label
         # tabl.runlemke(verbose=True, lexstats=True, z0=gz0)
@@ -167,6 +320,30 @@ def runLH(self, droppedlabel):
         return tuple(equilibrium)
 
     def LH(self, LHstring):
+        """Find equilibria by running Lemke-Howson from several labels.
+
+        Parameters
+        ----------
+        LHstring : str
+            Range specification of labels to try, e.g. ``"1-3,7"``.
+
+        Returns
+        -------
+        dict
+            Mapping from each distinct equilibrium found (as a tuple)
+            to the list of labels that produced it.
+            Also printed to stdout as equilibria are found.
+
+        Examples
+        --------
+        Run Lemke-Howson for all labels. This prints progress as each label
+        is tried, then a summary of all distinct equilibria found with
+        the list of labels that produced it:
+
+        >>> from lemke.bimatrix import bimatrix
+        >>> game = bimatrix.from_file("game.txt")
+        >>> game.LH("1-")
+        """
         m = self.A.numrows
         n = self.A.numcolumns
         lhset = {}  # dict of equilibria and list by which label found
@@ -184,6 +361,25 @@ def LH(self, LHstring):
         return lhset
 
     def runtrace(self, xprior, yprior):
+        """Run the tracing procedure from a given prior.
+
+        Parameters
+        ----------
+        xprior : array_like
+            Row player's prior strategy (probabilities summing to 1).
+        yprior : array_like
+            Column player's prior strategy (probabilities summing to 1).
+
+        Returns
+        -------
+        tuple of fractions.Fraction
+            The equilibrium strategy pair found, as a tuple of length ``m + n``.
+
+        Raises
+        ------
+        RuntimeError
+            If the underlying Lemke's algorithm unexpectedly fails to find a solution.
+        """
         lcp = self.createLCP()
         Ay = self.A.negmatrix @ yprior
         xB = xprior @ self.B.negmatrix
@@ -197,6 +393,19 @@ def runtrace(self, xprior, yprior):
         return tuple(equilibrium)
 
     def trace_uniform_prior(self):
+        """Run the tracing procedure once, starting from a uniform prior.
+
+        Prints the resulting equilibrium.
+
+        Examples
+        --------
+        Run the tracing procedure from a uniform prior
+        and print the equilibrium found:
+
+        >>> from lemke.bimatrix import bimatrix
+        >>> game = bimatrix.from_file("game.txt")
+        >>> game.trace_uniform_prior()
+        """
         m = self.A.numrows
         n = self.A.numcolumns
 
@@ -210,6 +419,35 @@ def trace_uniform_prior(self):
         )
 
     def trace_random_priors(self, trace, seed=None, accuracy=1000):
+        """Run the tracing procedure from several random priors.
+
+        Parameters
+        ----------
+        trace : int
+            Number of random priors to try. Must be positive.
+        seed : int, optional
+            Random seed. If given, each prior `k` is seeded
+            deterministically from `seed` and `k`, so results are
+            reproducible. If None, priors are unseeded.
+            Default is None.
+        accuracy : int, optional
+            Denominator used to round each random prior's coordinates
+            to rational numbers. Default is 1000.
+
+        Raises
+        ------
+        ValueError
+            If `trace` is not a positive integer.
+
+        Examples
+        --------
+        Run the tracing procedure from one random prior
+        and print the equilibrium found:
+
+        >>> from lemke.bimatrix import bimatrix
+        >>> game = bimatrix.from_file("game.txt")
+        >>> game.trace_random_priors(1)
+        """
         if trace <= 0:
             raise ValueError("Number of priors must be a positive integer")
         m = self.A.numrows
@@ -234,6 +472,7 @@ def trace_random_priors(self, trace, seed=None, accuracy=1000):
         self._print_trace_results(trset, trace)
 
     def _print_trace_results(self, equilibria, total_priors):
+        """Print a summary of equilibria found by a tracing run."""
         m = self.A.numrows
         n = self.A.numcolumns
 
@@ -243,6 +482,24 @@ def _print_trace_results(self, equilibria, total_priors):
         print(total_priors, "total priors,", len(equilibria), "equilibria found")
 
     def eqindex(self, eq, m, n):
+        """Compute the index (+1/-1/0) of an equilibrium.
+
+        Parameters
+        ----------
+        eq : sequence of fractions.Fraction
+            Equilibrium strategy pair, flat tuple of length ``m + n``.
+        m : int
+            Number of rows (row player's strategies).
+        n : int
+            Number of columns (column player's strategies).
+
+        Returns
+        -------
+        int
+            ``+1`` or ``-1`` for a regular equilibrium,
+            ``0`` if the equilibrium is degenerate
+            (supports of unequal size or singular submatrices).
+        """
         rowset, colset = supports(eq, m, n)
         k, ell = len(rowset), len(colset)
         if k != ell:
@@ -260,10 +517,12 @@ def eqindex(self, eq, m, n):
 
 
 def uniform(n):
+    """Return the uniform distribution over `n` outcomes as fractions."""
     return np.array([fractions.Fraction(1, n) for _ in range(n)])
 
 
 def str_eq(eq, m, n):
+    """Format an equilibrium strategy pair as a readable string."""
     x = "(" + ",".join([str(x) for x in eq[0:m]]) + ")"
     y = "(" + ",".join([str(x) for x in eq[m: m + n]]) + ")"
     rowset, colset = supports(eq, m, n)
@@ -271,12 +530,14 @@ def str_eq(eq, m, n):
 
 
 def supports(eq, m, n):
+    """Return the index sets where each player's strategy is nonzero."""
     rowset = [i for i in range(m) if eq[i] != 0]
     colset = [j for j in range(n) if eq[m + j] != 0]
     return rowset, colset
 
 
 def submatrix(A, rowset, colset):
+    """Extract a submatrix at the given row and column indices."""
     k, ell = len(rowset), len(colset)
     B = np.zeros((k, ell))
     for i in range(k):
diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index 5735c3d..91f5d14 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -1,4 +1,4 @@
-# LCP solver
+"""Lemke's algorithm for solving linear complementarity problems (LCPs)."""
 
 import fractions
 import math  # gcd
@@ -9,8 +9,46 @@
 from . import columnprint, utils
 
 
-# LCP data M,q,d
 class lcp:
+    r"""A linear complementarity problem instance `(M, q, d)`.
+
+    Represents the LCP :math:`w = q + Mz,\ w, z \ge 0,\ z^T w = 0`,
+    together with a covering vector `d` used to start Lemke's algorithm.
+
+    Can be constructed directly from `M`, `q`, `d`,
+    or via `from_file` for problems stored in the LCP file format.
+
+    Parameters
+    ----------
+    M : list of list of fractions.Fraction
+        Square matrix of shape `(n, n)`.
+    q : list of fractions.Fraction
+        Vector of length `n`.
+    d : list of fractions.Fraction
+        Covering vector of length `n`.
+
+    Attributes
+    ----------
+    M : list of list of fractions.Fraction
+        As passed in.
+    q, d : list of fractions.Fraction
+        As passed in.
+    n : int
+        Dimension of the problem.
+
+    Examples
+    --------
+    >>> from fractions import Fraction as F
+    >>> from lemke.lemke import lcp
+    >>> problem = lcp(
+    ...     M=[[F(2), F(1)], [F(1), F(2)]],
+    ...     q=[F(-1), F(-1)],
+    ...     d=[F(1), F(1)],
+    ... )
+    >>> problem.n
+    2
+    """
+
     def __init__(self, M, q, d):
         self.M = M
         self.q = q
@@ -19,7 +57,46 @@ def __init__(self, M, q, d):
 
     @classmethod
     def from_file(cls, filename):
-        # create LCP from file
+        """Create an LCP instance from a text file.
+
+        Expects a file starting with ``n= <dim>``,
+        followed by keyword-labeled blocks ``M=``, ``q=``, ``d=``
+        giving the matrix and vectors as whitespace-separated numbers.
+        Numbers may be given as integers, fractions (e.g. ``1/3``),
+        or decimals (e.g. ``0.3``).
+        ``#`` starts a comment that runs to the end of the line.
+
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Path to the LCP file.
+
+        Returns
+        -------
+        lcp
+            The LCP described by the file.
+
+        Raises
+        ------
+        ValueError
+            If the file does not start with ``n=``, doesn't contain
+            the expected number of values, or has an unrecognized keyword.
+
+        Examples
+        --------
+        Given a file ``lcp.txt`` containing::
+
+            n= 2
+            M= 2 1
+               1 2
+            q= -1 -1
+            d= 1 1
+
+        >>> from lemke.lemke import lcp
+        >>> problem = lcp.from_file("lcp.txt")
+        >>> print(problem.n)
+        2
+        """
         lines = utils.stripcomments(filename)
         # flatten into words
         words = utils.towords(lines)
@@ -85,7 +162,42 @@ def __str__(self):
 
 
 class tableau:
-    # filling the tableau from the LCP instance Mqd
+    """The tableau used to run Lemke's pivoting algorithm.
+
+    Stores the LCP data as a scaled integer tableau
+    (to allow exact arithmetic during pivoting)
+    and keeps track of which variables are currently basic/cobasic.
+
+    Parameters
+    ----------
+    Mqd : lcp
+        The LCP instance to build the tableau from.
+
+    Attributes
+    ----------
+    n : int
+        Problem dimension.
+    A : list of list of int
+        The scaled tableau, shape `(n, n + 2)`
+        (columns correspond to `n + 1` cobasic variables, plus RHS).
+    determinant : int
+        Current tableau determinant (always positive after pivoting).
+    scalefactor : list of int
+        Per-column scale factors (LCM of denominators) used to keep
+        the tableau in exact integer arithmetic.
+    bascobas : list of int
+        Location for each variable index: its tableau row (if basic) or
+        ``n + column`` (if cobasic).
+    whichvar : list of int
+        Inverse of `bascobas`: which variable occupies each row or column.
+    solution : list of fractions.Fraction
+        Current solution vector, filled in by `createsol`.
+    pivotcount : int
+        Number of pivots performed so far.
+    lextested, lexcomparisons : list of int
+        Per-column statistics from the lexicographic minimum-ratio test.
+    """
+
     def __init__(self, Mqd):
         self.n = Mqd.n
         n = self.n
@@ -170,13 +282,15 @@ def __str__(self):
         out += "\n" + "-----------------end of tableau-----------------"
         return out
 
-    def vartoa(self, v):  # variable as as string w1..wn or z0..zn
+    def vartoa(self, v):
+        """Return the display name of variable index `v`, e.g. "z0" or "w3"."""
         if v > self.n:
             return "w" + str(v - self.n)
         else:
             return "z" + str(v)
 
-    def createsol(self):  # get solution from current tableau
+    def createsol(self):
+        """Get solution from current tableau."""
         n = self.n
         for i in range(2 * n + 1):
             row = self.bascobas[i]
@@ -191,19 +305,22 @@ def createsol(self):  # get solution from current tableau
             else:  # i is nonbasic
                 self.solution[i] = fractions.Fraction(0)
 
-    def assertbasic(self, v, info):  # assert that v is basic
+    def assertbasic(self, v, info):
+        """Assert that variable v is basic."""
         if self.bascobas[v] >= self.n:
             raise RuntimeError(
                 f"({info}) Cobasic variable {self.vartoa(v)} should be basic"
             )
 
-    def assertcobasic(self, v, info):  # assert that v is cobasic
+    def assertcobasic(self, v, info):
+        """Assert that variable v is cobasic."""
         if self.bascobas[v] < self.n:
             raise RuntimeError(
                 f"({info}) Basic variable {self.vartoa(v)} should be cobasic"
             )
 
     def testtablvars(self):
+        """Check that `bascobas` and `whichvar` are consistent inverses."""
         n = self.n
         for i in range(2 * n + 1):
             if self.bascobas[self.whichvar[i]] != i:
@@ -218,7 +335,24 @@ def testtablvars(self):
                     )
                 raise RuntimeError(f"testtablvars() failed:\n{message}")
 
-    def complement(self, v):  # Z(i),W(i) are complements
+    def complement(self, v):
+        """Return the complementary variable of `v` (Z(i) <-> W(i)).
+
+        Parameters
+        ----------
+        v : int
+            A variable index other than ``z0``.
+
+        Returns
+        -------
+        int
+            The complementary variable index.
+
+        Raises
+        ------
+        RuntimeError
+            If `v` is ``z0`` (0), which has no complement.
+        """
         n = self.n
         if v == 0:
             raise RuntimeError("Attempt to find complement of z0")
@@ -227,14 +361,28 @@ def complement(self, v):  # Z(i),W(i) are complements
         else:
             return v + n
 
-    # returns leave,z0leave
-    # leave = leaving variable in VARS, given by lexmin row,
-    # when enter in VARS is entering variable
-    # only positive entries of entering column tested.
-    # Boolean z0leave indicates that z0 can leave the
-    # basis, but the lex-minratio test is performed fully,
-    # so  leave  might not be the index of z0
     def lexminvar(self, enter):
+        """Find the leaving variable via the lexicographic minimum-ratio test.
+
+        Parameters
+        ----------
+        enter : int
+            The entering variable index (must currently be cobasic).
+
+        Returns
+        -------
+        leave : int
+            The leaving variable index.
+        z0leave : bool
+            True if ``z0`` is among the (possibly tied) leaving candidates,
+            meaning the algorithm may terminate after this pivot.
+
+        Raises
+        ------
+        RayTermination
+            If no row has a positive entry in the entering column,
+            meaning the algorithm has found an unbounded ray instead of a solution.
+        """
         n = self.n
         A = self.A
         self.assertcobasic(enter, "Lexminvar")
@@ -298,21 +446,31 @@ def lexminvar(self, enter):
 
     # end of lexminvar(enter)
 
-    # negate tableau column  col
     def negcol(self, col):
+        """Negate every entry in the given column of the tableau."""
         for i in range(self.n):
             self.A[i][col] = -self.A[i][col]
 
-    # negate tableau row.  Used in  pivot()
     def negrow(self, row):
+        """Negate every entry in the given row of the tableau."""
         for j in range(self.n + 2):
             self.A[row][j] = -self.A[row][j]
 
-    # leave, enter in  VARS  defining  row, col  of  A
-    # pivot tableau on the element  A[row][col] which must be nonzero
-    # afterwards tableau normalized with positive determinant
-    # and updated tableau variables
     def pivot(self, leave, enter):
+        """Pivot the tableau, exchanging `leave` (basic) and `enter` (cobasic).
+
+        Performs an exact-arithmetic pivot on ``A[row][col]``
+        where `row`/`col` correspond to `leave`/`enter`,
+        then renormalizes the tableau to have positive determinant
+        and updates the basic/cobasic variables.
+
+        Parameters
+        ----------
+        leave : int
+            Variable currently basic that will become cobasic.
+        enter : int
+            Variable currently cobasic that will become basic.
+        """
         n = self.n
         A = self.A
         row = self.bascobas[leave]
@@ -354,6 +512,24 @@ def pivot(self, leave, enter):
 
 
 class RayTermination(Exception):
+    """Raised when Lemke's algorithm can't find a solution
+    and terminates on a secondary ray.
+
+    Parameters
+    ----------
+    enter : int
+        The variable that could not enter the basis
+        (no positive pivot candidate).
+    tableau : tableau
+        The tableau at the point of termination.
+        Its `solution` is populated before the exception is raised.
+
+    Attributes
+    ----------
+    tableau : tableau
+        The incomplete tableau at termination.
+    """
+
     def __init__(self, enter, tableau):
         tableau.createsol()
         self.tableau = tableau
@@ -362,7 +538,8 @@ def __init__(self, enter, tableau):
         )
 
 
-def outsol(tableau):  # string giving solution, after createsol()
+def outsol(tableau):
+    """Format the solution vector [z0, z1, ..., zn, w1, ..., wn] into columns."""
     # printout in columns to check complementarity
     n = tableau.n
     sol = columnprint.columnprint(n + 2)
@@ -384,8 +561,8 @@ def outsol(tableau):  # string giving solution, after createsol()
     return str(sol)
 
 
-# output statistics of minimum ratio test
 def outstatistics(tableau):
+    """Helper to output statistics of minimum ratio test."""
     n = tableau.n
     lext = tableau.lextested
     stats = columnprint.columnprint(n + 2)
@@ -411,23 +588,49 @@ def outstatistics(tableau):
 
 
 class LemkeCallback:
-    def on_start(self, lcp, tableau): pass
-    def on_negcol(self, tableau): pass
-    def on_pivot_start(self, tableau, leave, enter): pass
-    def on_pivot_end(self, tableau): pass
-    def on_done(self, tableau): pass
-    def on_ray_termination(self, tableau, message): pass
+    """Callback interface for observing the progress of `runlemke`.
+
+    Subclass and override any of these methods to log, print, or
+    otherwise react to the algorithm's progress (e.g. `PrintingCallback`).
+    All methods are no-ops by default.
+    """
+
+    def on_start(self, lcp, tableau):
+        """Called once, after the initial tableau is built."""
+
+    def on_negcol(self, tableau):
+        """Called after the RHS column is negated to start pivoting."""
+
+    def on_pivot_start(self, tableau, leave, enter):
+        """Called before each pivot, with the chosen leave/enter variables."""
+
+    def on_pivot_end(self, tableau):
+        """Called after each pivot completes."""
+
+    def on_done(self, tableau):
+        """Called once a complementary solution is found."""
+
+    def on_ray_termination(self, tableau, message):
+        """Called if the algorithm terminates on a secondary ray."""
 
 
 class PrintingCallback(LemkeCallback):
-    # z0: printout value of z0
-    # flags.maxcount   = 0;
-    # flags.bdocupivot = 1;
-    # flags.binitabl   = 1;
-    # flags.bouttabl   = 0;  (= verbose)
-    # flags.boutsol    = 1;
-    # flags.binteract  = 0;
-    # flags.blexstats  = 0;
+    """A `LemkeCallback` that prints tableaus and progress to a stream.
+
+    Parameters
+    ----------
+    stream : file-like, optional
+        Where to print output. Default is `sys.stdout`.
+    verbose : bool, optional
+        If True, print the full tableau after every pivot,
+        not just at the start and end. Default is False.
+    z0 : bool, optional
+        If True, print the current value of ``z0`` before each pivot.
+        Default is False.
+    lexstats : bool, optional
+        If True, print lexicographic minimum-ratio test statistics
+        when the algorithm finishes successfully. Default is False.
+    """
 
     def __init__(
         self,
@@ -442,9 +645,11 @@ def __init__(
         self.lexstats = lexstats
 
     def printout(self, *args):
+        """Print `args` to this callback's stream."""
         print(*args, file=self.stream)
 
     def on_start(self, lcp, tableau):
+        """Print the LCP instance and the initial tableau."""
         self.printout(f"verbose={self.verbose} z0={self.z0} lexstats={self.lexstats}")
         self.printout(lcp)
         self.printout("==================================")
@@ -454,12 +659,16 @@ def on_start(self, lcp, tableau):
         self.printout(tableau)
 
     def on_negcol(self, tableau):
+        """Print the tableau after negating the RHS column, if `verbose`."""
         # if (flags.binitabl)
         if self.verbose:
             self.printout("After negcol:")
             self.printout(tableau)
 
     def on_pivot_start(self, tableau, leave, enter):
+        """Print the chosen leaving/entering variables,
+        and `z0`'s value if `z0` is set.
+        """
         if self.z0:  # printout progress of z0
             z0_value = 0.0
             if tableau.bascobas[0] < tableau.n:  # z0 is basic
@@ -470,10 +679,14 @@ def on_pivot_start(self, tableau, leave, enter):
         self.printout(f"leaving: {leave.ljust(5)} entering: {enter}")
 
     def on_pivot_end(self, tableau):
+        """Print the tableau after the pivot completes, if `verbose`."""
         if self.verbose:
             self.printout(tableau)
 
     def on_done(self, tableau):
+        """Print the final tableau, solution,
+        and lex-stats if `lexstats` is set.
+        """
         if self.z0:
             self.printout(f"pivot count = {tableau.pivotcount + 1}, z0 = 0.0")
 
@@ -489,6 +702,9 @@ def on_done(self, tableau):
             self.printout(outstatistics(tableau))
 
     def on_ray_termination(self, tableau, message):
+        """Print the ray-termination message, tableau,
+        and current (incomplete) solution.
+        """
         self.printout(message)
         self.printout(tableau)
         self.printout("Current basis not an LCP solution:")
@@ -496,6 +712,45 @@ def on_ray_termination(self, tableau, message):
 
 
 def runlemke(*, lcp, callback=None):
+    """Solve an LCP using Lemke's complementary pivoting algorithm.
+
+    Parameters
+    ----------
+    lcp : lcp
+        The LCP instance to solve.
+    callback : LemkeCallback, optional
+        Hook invoked at each step of the algorithm to observe or report progress.
+        Pass `PrintingCallback` for a ready-made implementation. Default is None.
+
+    Returns
+    -------
+    list of fractions.Fraction or None
+        The full solution vector (length ``2n + 1``,
+        indices ``z0, z1, ..., zn, w1, ..., wn``)
+        if a complementary solution was found,
+        or None if the algorithm couldn't find a solution
+        and terminated on a secondary ray.
+
+    Examples
+    --------
+    Basic usage with no callbacks:
+
+    >>> from lemke.lemke import lcp, runlemke
+    >>> problem = lcp.from_file("lcp.txt")
+    >>> solution = runlemke(lcp=problem)
+    >>> print(solution)
+    [Fraction(0, 1), Fraction(2, 1), Fraction(1, 1), Fraction(0, 1), Fraction(0, 1)]
+
+    Usage with a printing callback:
+
+    >>> from lemke.lemke import lcp, runlemke, PrintingCallback
+    >>> problem = lcp.from_file("lcp.txt")
+    >>> cb = PrintingCallback(verbose=True, z0=True)
+    >>> solution = runlemke(lcp=problem, callback=cb)
+    # prints the given lcp, tableau and z0 at each step, and the final solution
+    >>> print(solution)
+    [Fraction(0, 1), Fraction(2, 1), Fraction(1, 1), Fraction(0, 1), Fraction(0, 1)]
+    """
     callback = callback or LemkeCallback()
 
     try: