diff --git a/doc/conf.py b/doc/conf.py
index 4e0c586..2426932 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -3,6 +3,7 @@
 author = "Bernhard von Stengel"
 
 extensions = [
+    "matplotlib.sphinxext.plot_directive",
     "sphinx.ext.autodoc",
     "sphinx.ext.autosummary",
     "sphinx.ext.doctest",
@@ -34,3 +35,6 @@
     ],
     "shorten_urls": False,
 }
+
+plot_html_show_source_link = False
+plot_html_show_formats = False
diff --git a/doc/index.rst b/doc/index.rst
index 597118f..2a64a8e 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -141,3 +141,4 @@ Useful links
 
    api
    file_formats
+   tutorials/index
diff --git a/doc/tutorials/bimatrix.rst b/doc/tutorials/bimatrix.rst
new file mode 100644
index 0000000..a41428d
--- /dev/null
+++ b/doc/tutorials/bimatrix.rst
@@ -0,0 +1,188 @@
+Finding equilibria of a bimatrix game
+=====================================
+
+This tutorial shows how to define a two-player game in strategic form and
+compute its Nash equilibria with the Lemke-Howson algorithm and with the
+tracing procedure.
+
+In a *bimatrix game*, the row player has :math:`m` pure strategies, the
+column player has :math:`n`, and the payoffs are given by two
+:math:`m \times n` matrices :math:`A` (row player) and :math:`B` (column
+player). All payoffs are stored as exact fractions.
+
+Throughout we use this coordination game:
+
+.. math::
+
+   A = \begin{pmatrix} 3 & 0 \\ 0 & 2 \end{pmatrix}, \qquad
+   B = \begin{pmatrix} 2 & 0 \\ 0 & 7 \end{pmatrix}.
+
+
+Creating a game
+---------------
+
+A game is represented by the :class:`~lemke.bimatrix.bimatrix` class,
+which contains two :class:`~lemke.bimatrix.payoffmatrix` objects, one per player.
+A game instance can be created in two ways.
+
+1) From Python
+~~~~~~~~~~~~~~
+
+First wrap each payoff matrix in a :class:`~lemke.bimatrix.payoffmatrix`,
+then combine them:
+
+.. doctest::
+
+   >>> from fractions import Fraction as F
+   >>> from lemke.bimatrix import bimatrix, payoffmatrix
+   >>> game = bimatrix(
+   ...    A=payoffmatrix([
+   ...       [F(3), F(0)],
+   ...       [F(0), F(2)],
+   ...    ]),
+   ...    B=payoffmatrix([
+   ...       [F(2), F(0)],
+   ...       [F(0), F(7)],
+   ...    ]),
+   ... )
+   >>> print(game.A.numrows)
+   2
+
+A :class:`~lemke.bimatrix.payoffmatrix` also records useful facts about its entries,
+for example its largest entry:
+
+.. doctest::
+
+   >>> print(game.A.max)
+   3
+
+.. note::
+
+   Every entry must be an exact rational (:class:`fractions.Fraction`).
+   Other types cannot be passed directly.
+
+2) From a file
+~~~~~~~~~~~~~~
+
+Alternatively, describe the game in a text file and load it with
+:meth:`~lemke.bimatrix.bimatrix.from_file`.
+
+First, create a file such as ``examples/game``::
+
+   # m,n=
+   2 2
+
+   # A=
+   3 0
+   0 2
+
+   # B=
+   2 0
+   0 7
+
+then load it:
+
+.. doctest::
+
+   >>> game = bimatrix.from_file("examples/game")
+   >>> print(game.A.numrows)
+   2
+
+See details about how the input file should be formatted on
+the :doc:`file formats<../file_formats>` page.
+
+
+Lemke-Howson
+------------
+
+The Lemke-Howson algorithm finds one equilibrium for each starting *label*.
+There are :math:`m + n` labels, one for each pure strategy of either player.
+:meth:`~lemke.bimatrix.bimatrix.LH` takes a string saying which labels to
+try, for example ``"1-3,7"``. You can also pass ``"1-"`` for all labels:
+
+.. doctest::
+
+   >>> result = game.LH("1-")
+   label 1 found eq (1,0),(1,0)
+       supports: [0][0]
+   label 2 found eq (0,1),(0,1)
+       supports: [1][1]
+   -------- equilibria found: --------
+   (1,0),(1,0)
+       supports: [0][0] found by labels [1, 3]
+   (0,1),(0,1)
+       supports: [1][1] found by labels [2, 4]
+
+This prints a line each time a *new* equilibrium is found, followed by a
+summary of all distinct equilibria and the labels that found them.
+
+The same information is returned as a dictionary, which maps each
+equilibrium (a tuple) to the list of labels that produced it:
+
+.. doctest::
+
+   >>> for eq, labels in result.items():
+   ...    print(labels)
+   ...    print(eq)
+   [1, 3]
+   (Fraction(1, 1), Fraction(0, 1), Fraction(1, 1), Fraction(0, 1))
+   [2, 4]
+   (Fraction(0, 1), Fraction(1, 1), Fraction(0, 1), Fraction(1, 1))
+
+.. note::
+
+   Different labels can lead to the same equilibrium,
+   and Lemke-Howson may not find *all* equilibria of a game.
+   For this game, it found two pure strategy equilibria, but
+   not the mixed strategy one: :math:`(\frac{7}{9}, \frac{2}{9}), (\frac{2}{5}, \frac{3}{5})`.
+
+
+.. _tracing:
+
+The tracing procedure
+---------------------
+
+The tracing procedure is another way to find an equilibrium. It starts from
+a *prior:* a pair of mixed strategies, one per player, that represents
+initial beliefs. The equilibrium found depends on the prior.
+
+Uniform prior
+~~~~~~~~~~~~~
+
+:meth:`~lemke.bimatrix.bimatrix.trace_uniform_prior` runs the procedure once
+from the uniform prior, where each player mixes equally over all their
+pure strategies, and prints the equilibrium found:
+
+.. doctest::
+
+   >>> game.trace_uniform_prior()
+   -------- statistics of equilibria found: --------
+   1 times found (0,1),(0,1)
+       supports: [1][1]
+   1 total priors, 1 equilibria found
+
+Random priors
+~~~~~~~~~~~~~
+
+:meth:`~lemke.bimatrix.bimatrix.trace_random_priors` runs the procedure from
+several random priors, to see which equilibria they lead to and how often.
+To make a run reproducible, pass a ``seed``:
+
+.. doctest::
+
+   >>> game.trace_random_priors(100, seed=1, accuracy=100)
+   found eq (1,0),(1,0)
+       supports: [0][0] index 1
+   found eq (0,1),(0,1)
+       supports: [1][1] index 1
+   -------- statistics of equilibria found: --------
+   24 times found (1,0),(1,0)
+       supports: [0][0]
+   76 times found (0,1),(0,1)
+       supports: [1][1]
+   100 total priors, 2 equilibria found
+
+The priors are sampled uniformly from the simplex and
+rounded to fractions with denominator ``accuracy`` (default 1000).
+
+Check :doc:`this tutorial<randomstart>` to see how the sampling and rounding work.
diff --git a/doc/tutorials/index.rst b/doc/tutorials/index.rst
new file mode 100644
index 0000000..2cf92e2
--- /dev/null
+++ b/doc/tutorials/index.rst
@@ -0,0 +1,9 @@
+Tutorials
+=========
+
+.. toctree::
+   :maxdepth: 2
+
+   lcp
+   bimatrix
+   randomstart
diff --git a/doc/tutorials/lcp.rst b/doc/tutorials/lcp.rst
new file mode 100644
index 0000000..761ce3d
--- /dev/null
+++ b/doc/tutorials/lcp.rst
@@ -0,0 +1,300 @@
+Solving a linear complementarity problem
+========================================
+
+This tutorial walks through solving a *linear complementarity problem* (LCP)
+with Lemke's algorithm, from creating the problem to interpreting the result
+and watching the algorithm run.
+
+Given a vector :math:`q` and a matrix :math:`M`,
+an LCP asks to find a vector :math:`z \ge 0` such that
+
+.. math::
+
+   w = q + Mz \ge 0, \qquad z^T w = 0.
+
+To get started, Lemke's algorithm also needs a *covering vector* :math:`d`,
+which multiplies the extra variable :math:`z_0`.
+So, a problem is described by the triple :math:`(M, q, d)`.
+All arithmetic is exact, using rational numbers.
+
+
+Creating an LCP
+---------------
+
+An LCP is represented by the :class:`~lemke.lemke.lcp` class.
+It can be created in two ways.
+
+Let's consider an LCP, where
+
+.. math::
+
+   M = \begin{pmatrix} \frac{1}{2} & 1 \\ -1 & \frac{5}{2} \end{pmatrix},
+   \qquad
+   q = \begin{pmatrix} -1 \\ -1 \end{pmatrix},
+   \qquad
+   d = \begin{pmatrix} 1 \\ 1 \end{pmatrix}.
+
+
+1) From Python
+~~~~~~~~~~~~~~
+
+Pass the matrix ``M`` (as a list of its rows) and the vectors ``q`` and ``d``:
+
+.. doctest::
+
+   >>> from fractions import Fraction as F
+   >>> from lemke.lemke import lcp
+   >>> problem = lcp(
+   ...    M=[
+   ...       [F(1, 2), F(1)],
+   ...       [F(-1), F(5, 2)],
+   ...    ],
+   ...    q=[F(-1), F(-1)],
+   ...    d=[F(1), F(1)],
+   ... )
+   >>> print(problem.n)
+   2
+
+.. note::
+
+   Entries must be exact rationals, such as Python integers or
+   :class:`fractions.Fraction` instances. Floats are not supported.
+
+2) From a file
+~~~~~~~~~~~~~~
+
+Alternatively, describe the problem in a text file and load it with
+:meth:`~lemke.lemke.lcp.from_file`.
+
+First, create a file such as ``examples/lcp_2x2``::
+
+   n= 2
+   M= 1/2 1
+     -1 2.5
+   q= -1 -1
+   d= 1 1
+
+then load it:
+
+.. doctest::
+
+   >>> problem = lcp.from_file("examples/lcp_2x2")
+   >>> print(problem.n)
+   2
+
+See details about how the input file should be formatted on
+the :doc:`file formats<../file_formats>` page.
+
+
+Solving the problem
+-------------------
+
+Pass the problem to :func:`~lemke.lemke.runlemke`:
+
+.. doctest::
+
+   >>> from lemke.lemke import runlemke
+   >>> result = runlemke(lcp=problem)
+   >>> print(result.success)
+   True
+
+:func:`~lemke.lemke.runlemke` always returns an :class:`~lemke.lemke.LcpResult`,
+whether or not a solution was found. Its fields are the same in both cases,
+so you must check :attr:`~lemke.lemke.LcpResult.success` field first.
+
+.. note::
+
+   If every entry of ``q`` is non-negative, the trivial solution
+   :math:`z = 0,\ w = q` is returned immediately without any pivoting.
+
+
+Reading the result
+------------------
+
+The solution values are available as tuples of :class:`fractions.Fraction`:
+
+.. doctest::
+
+   >>> print(result.z)
+   (Fraction(2, 3), Fraction(2, 3))
+   >>> print(result.w)
+   (Fraction(0, 1), Fraction(0, 1))
+
+Printing the result shows the basis and the values in aligned columns,
+so that the complementarity condition (:math:`z_i w_i = 0`) can be checked by eye:
+
+.. doctest::
+
+   >>> print(result)
+   Process finished successfully after 3 pivots.
+   Solution found:
+   basis=     z1  z2
+       z=  0 2/3 2/3
+       w=      0   0
+
+The first column of ``z=`` is the artificial variable :math:`z_0`, which is
+always 0 for a successful run. See :class:`~lemke.lemke.LcpResult` for all
+fields, such as :attr:`~lemke.lemke.LcpResult.basis` and :attr:`~lemke.lemke.LcpResult.num_pivots`.
+
+Degenerate solutions
+~~~~~~~~~~~~~~~~~~~~
+
+The solution above is *non-degenerate*: for every :math:`i`, exactly one of
+:math:`z_i` and :math:`w_i` is positive.
+
+A solution is *degenerate* if some basic variable has value zero, meaning
+both :math:`z_i` and :math:`w_i` are zero for some :math:`i`.
+It is still a valid solution. For example:
+
+.. doctest::
+
+   >>> problem = lcp(
+   ...    M=[
+   ...       [F(1), F(0)],
+   ...       [F(0), F(1)],
+   ...    ],
+   ...    q=[F(-1), F(0)],
+   ...    d=[F(1), F(1)],
+   ... )
+   >>> result = runlemke(lcp=problem)
+   >>> print(result.z)
+   (Fraction(1, 1), Fraction(0, 1))
+   >>> print(result.w)
+   (Fraction(0, 1), Fraction(0, 1))
+
+Here :math:`z_2 = w_2 = 0`, so the second complementary pair is degenerate.
+Degeneracy is handled by the lexicographic minimum ratio test, which
+prevents the algorithm from cycling. You can inspect how often that test
+was needed with the ``lexstats`` option of the printing callback (see below).
+
+Since in this case it is not clear which variables are basic,
+you can explicitly check the :attr:`~lemke.lemke.LcpResult.basis` field:
+
+.. doctest::
+
+   >>> print(sorted(result.basis))
+   ['w2', 'z1']
+
+
+Ray termination
+---------------
+
+Lemke's algorithm can also end on a *secondary ray*, when the variable
+that should enter the basis has no positive pivot candidate.
+In that case no valid solution is found.
+
+Consider the following example, where :math:`w = -1 + 0 * z` can
+never be non-negative:
+
+.. doctest::
+
+   >>> problem = lcp(
+   ...   M=[[F(0)]],
+   ...   q=[F(-1)],
+   ...   d=[F(1)],
+   ... )
+   >>> result = runlemke(lcp=problem)
+   >>> print(result.success)
+   False
+   >>> print(result.ray_entering_variable)
+   z1
+
+When ``success`` is ``False``:
+
+* ``ray_entering_variable`` names the entering variable that led to the secondary ray,
+* ``z0``, ``z`` and ``w`` describe the final values reached,
+  which is **not** a solution of the LCP.
+
+.. doctest::
+
+   >>> print(result)
+   Terminated on a secondary ray after 1 pivot, when trying to enter z1.
+   Current basis not an LCP solution:
+   basis= z0
+       z=  1  0
+       w=     0
+
+.. note::
+
+   Ray termination does not, by itself, prove that the LCP has no solution.
+   It only means the algorithm could not find one from this covering vector.
+
+
+Watching the algorithm
+----------------------
+
+To see what happens at each pivot, pass a callback to :func:`~lemke.lemke.runlemke`.
+The ready-made :class:`~lemke.lemke.PrintingCallback` prints the LCP, the
+tableaus and the chosen leaving and entering variables:
+
+.. doctest::
+
+   >>> from lemke.lemke import PrintingCallback
+   >>> problem = lcp.from_file("examples/lcp_2x2")
+   >>> result = runlemke(lcp=problem, callback=PrintingCallback())
+   verbose=False z0=False lexstats=False
+   n= 2
+   M=
+   1/2   1
+   -1  5/2
+   q=
+   -1   -1
+   d=
+   1     1
+   ==================================
+   After filltableau:
+   Determinant: -1
+   var  z0 z1 z2 RHS
+   scfa  1  2  2   1
+   <BLANKLINE>
+   w1    1  1  2  -1
+   w2    1 -2  5  -1
+   -----------------end of tableau-----------------
+   leaving: w2    entering: z0
+   leaving: w1    entering: z2
+   leaving: z0    entering: z1
+   Final tableau:
+   Determinant: 9
+   var  w2 z0 w1 RHS
+   scfa  1  1  1   1
+   <BLANKLINE>
+   z2   -1  3 -2   3
+   z1    2  3 -5   3
+   -----------------end of tableau-----------------
+   Process finished successfully after 3 pivots.
+   Solution found:
+   basis=     z1  z2
+       z=  0 2/3 2/3
+       w=      0   0
+
+The options are:
+
+* ``verbose``: print the full tableau after every pivot,
+* ``z0``: print the current value of :math:`z_0` before each pivot,
+* ``lexstats``: print minimum-ratio test statistics at the end,
+* ``stream``: write to a file-like object instead of ``sys.stdout``.
+
+Writing your own callback
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+To react to the algorithm in your own way, subclass
+:class:`~lemke.lemke.LemkeCallback` and override only the methods you need.
+All of them do nothing by default.
+
+For example, to log entering variables during each pivot,
+you can override :func:`~lemke.lemke.LemkeCallback.on_pivot_start`:
+
+.. doctest::
+
+   >>> from lemke.lemke import LemkeCallback
+   >>> class LogEnteringVariables(LemkeCallback):
+   ...    def on_pivot_start(self, tableau, leave, enter):
+   ...        print(enter)
+   >>> result = runlemke(lcp=problem, callback=LogEnteringVariables())
+   z0
+   z2
+   z1
+
+.. note::
+
+   No callbacks are invoked in the trivial case :math:`q \ge 0`.
diff --git a/doc/tutorials/randomstart.rst b/doc/tutorials/randomstart.rst
new file mode 100644
index 0000000..c0d11b3
--- /dev/null
+++ b/doc/tutorials/randomstart.rst
@@ -0,0 +1,117 @@
+Visualizing random priors on the simplex
+========================================
+
+This tutorial shows how to use :func:`~lemke.randomstart.plot_simplex`
+function to visualize how random probability vectors are sampled and
+rounded to exact fractions. This is the same sampling used for
+random priors in the :ref:`tracing procedure <tracing>`.
+
+The set of all probability vectors with three entries (non-negative numbers
+that sum to 1) forms a triangle, called the *unit simplex*. Each probability
+vector is a point in this triangle, and its corners are the vectors
+:math:`(1, 0, 0)`, :math:`(0, 1, 0)` and :math:`(0, 0, 1)`.
+:func:`~lemke.randomstart.plot_simplex` samples random points
+from the simplex and draws two things:
+
+* the raw sampled points as **green dots**,
+* their rounded approximations as **red circles**.
+
+
+The rounding turns each point into fractions with a common denominator
+(``accuracy``) that still sum to exactly 1, which is needed for
+the exact arithmetic in the algorithms.
+
+
+A first plot
+------------
+
+.. plot::
+   :include-source:
+
+   import random
+   from lemke.randomstart import plot_simplex
+
+   random.seed(1)
+   plot_simplex(numpoints=300, accuracy=30)
+
+``random.seed`` is called first just to make the plot reproducible.
+
+Every red circle sits on a grid point of the triangle, and each green dot
+is rounded to a nearby grid point. With a small ``accuracy``, the grid is
+coarse, so many dots share one circle.
+
+The function also prints the parameters it was called with:
+
+.. code-block:: python
+
+   numpoints=300 accuracy=30 higherdim=3 naiveplot=False
+
+
+Changing the accuracy
+---------------------
+
+The ``accuracy`` is the denominator of the rounding: coordinates become
+multiples of ``1 / accuracy``. Increasing it gives a finer grid, so the
+red circles (which are also drawn smaller) lie closer to the dots:
+
+.. plot::
+   :include-source:
+
+   import random
+   from lemke.randomstart import plot_simplex
+
+   random.seed(1)
+   plot_simplex(numpoints=300, accuracy=100)
+
+
+.. note::
+
+   ``accuracy`` must be an integer between 1 and 10,000,000.
+
+
+Uniform versus naive sampling
+-----------------------------
+
+By default the points are sampled *uniformly* from the simplex.
+With ``naiveplot=True``, points are instead sampled by normalizing
+independent uniform random numbers. This is a common shortcut,
+but it is **not** uniform. The points are denser near the centre
+of the triangle and sparser near the corners:
+
+.. plot::
+   :include-source:
+
+   import random
+   from lemke.randomstart import plot_simplex
+
+   random.seed(1)
+   plot_simplex(numpoints=300, accuracy=30, naiveplot=True)
+
+
+Higher dimensions
+-----------------
+
+A triangle can only show three components. With ``higherdim`` argument greater
+than 3, points are sampled from the higher-dimensional simplex, and only the
+*middle three* components of each point are kept and renormalized to sum to 1.
+
+.. plot::
+   :include-source:
+
+   import random
+   from lemke.randomstart import plot_simplex
+
+   random.seed(1)
+   plot_simplex(numpoints=300, accuracy=30, higherdim=6)
+
+
+The function also prints which positions are shown:
+
+.. code-block:: python
+
+   numpoints=300 accuracy=30 higherdim=6 naiveplot=False
+   show positions 2 .. 4 of 0 .. 5
+
+.. note::
+
+   ``higherdim`` must be an integer between 3 and 10.
diff --git a/src/lemke/columnprint.py b/src/lemke/columnprint.py
index 11fc876..47fe5b7 100644
--- a/src/lemke/columnprint.py
+++ b/src/lemke/columnprint.py
@@ -37,7 +37,7 @@ def prline(self, line):  # prline(char *s)
             out += word
             if i < self.ncols - 1:
                 out += " "  # single spacing here
-        return out
+        return out.rstrip()
 
     # print integer i into the current column
     def iprint(self, i):  # colipr(i)
diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index e710b0b..6335991 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -619,12 +619,15 @@ class PrintingCallback(LemkeCallback):
 
     def __init__(
         self,
-        stream=sys.stdout,
+        stream=None,
         verbose=False,
         z0=False,
         lexstats=False,
     ):
-        self.stream = stream
+        if stream is not None:
+            self.stream = stream
+        else:
+            self.stream = sys.stdout
         self.verbose = verbose
         self.z0 = z0
         self.lexstats = lexstats