diff --git a/pyproject.toml b/pyproject.toml
index 9c08bc7..95352ce 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,8 +6,9 @@ build-backend = "setuptools.build_meta"
 name = "lemke"
 version = "0.0.1"
 dependencies = [
-    "numpy>=2.2,<2.3",
-    "matplotlib>=3.10,<3.11"
+    "numpy>=2.2",
+    "matplotlib>=3.10",
+    "click>=8.1",
 ]
 
 requires-python = ">=3.10"
diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py
index 241e34c..b66f339 100644
--- a/src/lemke/bimatrix.py
+++ b/src/lemke/bimatrix.py
@@ -253,7 +253,7 @@ def runLH(self, droppedlabel):
         lcp.d[droppedlabel - 1] = 0  # subsidize this label
         tabl = lemke.tableau(lcp)
         # tabl.runlemke(verbose=True, lexstats=True, z0=gz0)
-        tabl.runlemke(silent=True)
+        tabl.runlemke()
         return tuple(getequil(tabl))
 
     def LH(self, LHstring):
@@ -281,7 +281,7 @@ def runtrace(self, xprior, yprior):
         xB = xprior @ self.B.negmatrix
         lcp.d = np.hstack((Ay, xB, [1, 1]))
         tabl = lemke.tableau(lcp)
-        tabl.runlemke(silent=True)
+        tabl.runlemke()
         return tuple(getequil(tabl))
 
     def tracing(self, trace):
diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index 4aad705..c65e042 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -4,49 +4,9 @@
 import math  # gcd
 import sys
 
-from . import columnprint, utils
-
-# global defaults
-lcpfilename = "lcp"
-outfile = lcpfilename + ".out"
-filehandle = sys.stdout
-verbose = False
-silent = False
-z0 = False
-
-
-# process command-line arguments
-def processArguments():
-    global lcpfilename, outfile, verbose, silent, z0
-    helpstring = """usage: lemke.py [options]
-options: -v, -verbose : printout intermediate tableaus
-         -s, -silent  : send output to <lcpfilename>.out
-         -z0 : show value of z0 at each step
-         -?, -help:    show this help
-         <lcpfilename> (default: "lcp", must not start with "-")
-         Example: [python] lemke.py -v 2lcp"""
-    arglist = sys.argv[1:]
-    showhelp = False
-    for s in arglist:
-        if s == "-v" or s == "-verbose":
-            verbose = True
-        elif s == "-s" or s == "-silent":
-            silent = True
-        elif s == "-z0":
-            z0 = True
-        elif s[0] == "-":
-            showhelp = True
-        else:
-            lcpfilename = s
-            outfile = s + ".out"
-    if showhelp:
-        printout(helpstring)
-        exit(0)
-    return
-
+import click
 
-def printout(*s):
-    print(*s, file=filehandle)
+from . import columnprint, utils
 
 
 # LCP data M,q,d
@@ -70,9 +30,9 @@ def __init__(self, arg):
             # flatten into words
             words = utils.towords(lines)
             if words[0] != "n=":
-                printout("lcp file", repr(filename),
-                         "must start with 'n=' lcpdim, e.g. 'n= 5', not",
-                         repr(words[0]))
+                print("lcp file", repr(filename),
+                      "must start with 'n=' lcpdim, e.g. 'n= 5', not",
+                      repr(words[0]))
                 exit(1)
             n = int(words[1])
             self.n = n
@@ -87,9 +47,9 @@ def __init__(self, arg):
             needfracs = n * n + 2 * n
             if len(words) != needfracs + 5:
                 # printout("in lcp file '",filename,"':")
-                printout("in lcp file " + repr(filename) + ":")
-                printout("n=", n, ", need keywords 'M=' 'q=' 'd=' and n*n + n + n =",
-                         needfracs, "fractions, got", len(words) - 5)
+                print("in lcp file " + repr(filename) + ":")
+                print("n=", n, ", need keywords 'M=' 'q=' 'd=' and n*n + n + n =",
+                      needfracs, "fractions, got", len(words) - 5)
                 exit(1)
             k = 2  # index in words
             while k < len(words):
@@ -106,8 +66,8 @@ def __init__(self, arg):
                     self.d = utils.tovector(n, words, k)
                     k += n
                 else:
-                    printout("in lcp file " + repr(filename) + ":")
-                    printout("expected one of 'M=' 'q=' 'd=', got", repr(words[k]))
+                    print("in lcp file " + repr(filename) + ":")
+                    print("expected one of 'M=' 'q=' 'd=', got", repr(words[k]))
                     exit(1)
             return
 
@@ -139,6 +99,7 @@ def __str__(self):
 class tableau:
     # filling the tableau from the LCP instance Mqd
     def __init__(self, Mqd):
+        self.output_stream = None
         self.n = Mqd.n
         n = self.n
         self.scalefactor = [0] * (n + 2)  # 0 for z0, n+1 for RHS
@@ -266,15 +227,15 @@ def outsol(self):  # string giving solution, after createsol()
 
     def assertbasic(self, v, info):  # assert that v is basic
         if self.bascobas[v] >= self.n:
-            printout(info, "Cobasic variable", self.vartoa(v),
-                     "should be basic")
+            self.printout(info, "Cobasic variable", self.vartoa(v),
+                          "should be basic")
             exit(1)
         return
 
     def assertcobasic(self, v, info):  # assert that v is cobasic
         if self.bascobas[v] < self.n:
-            printout(info, "Cobasic variable", self.vartoa(v),
-                     "should be cobasic")
+            self.printout(info, "Cobasic variable", self.vartoa(v),
+                          "should be cobasic")
             exit(1)
         return
 
@@ -283,15 +244,15 @@ def docupivot(self, leave, enter):  # leave, enter in VARS
         self.assertcobasic(enter, "docupivot")
         s = "leaving: " + self.vartoa(leave).ljust(5)
         s += "entering: " + self.vartoa(enter)
-        printout(s)
+        self.printout(s)
         return
 
     def raytermination(self, enter):
-        printout("Ray termination when trying to enter", self.vartoa(enter))
-        printout(self)
-        printout("Current basis not an LCP solution:")
+        self.printout("Ray termination when trying to enter", self.vartoa(enter))
+        self.printout(self)
+        self.printout("Current basis not an LCP solution:")
         self.createsol()
-        printout(self.outsol())
+        self.printout(self.outsol())
         exit(1)
 
     def testtablvars(self):  # msg only if error, continue
@@ -301,9 +262,9 @@ def testtablvars(self):  # msg only if error, continue
                 # injective suffices
                 for j in range(2 * n + 1):
                     if j == i:
-                        printout("First problem for j=", j, ":")
+                        self.printout("First problem for j=", j, ":")
                     # printout (f"{j=} {self.bascobas[j]=} {self.whichvar[j]=}")
-                    printout(
+                    self.printout(
                         f"j={j} self.bascobas[j]={self.bascobas[j]} "
                         f"self.whichvar[j]={self.whichvar[j]}"
                     )
@@ -312,7 +273,7 @@ def testtablvars(self):  # msg only if error, continue
     def complement(self, v):  # Z(i),W(i) are complements
         n = self.n
         if v == 0:
-            printout("Attempt to find complement of z0")
+            self.printout("Attempt to find complement of z0")
             exit(1)
         if v > n:
             return v - n
@@ -342,7 +303,7 @@ def outstatistics(self):
                     stats.sprint(str(x / 10.0))
                 else:
                     stats.sprint("-")
-        printout(stats)
+        self.printout(stats)
 
     # returns leave,z0leave
     # leave = leaving variable in VARS, given by lexmin row,
@@ -377,7 +338,7 @@ def lexminvar(self, enter):
         j = 0  # going through j = 0..n
         while len(leavecand) > 1:
             if j > n:  # impossible, perturbed RHS should have full rank
-                printout("lex-minratio test failed")
+                self.printout("lex-minratio test failed")
                 exit(1)
             self.lextested[j] += 1
             self.lexcomparisons[j] += len(leavecand)
@@ -468,8 +429,19 @@ def pivot(self, leave, enter):
 
     # end of  pivot (leave, enter)
 
-    def runlemke(self, *, verbose=False, lexstats=False, z0=False, silent=False):
-        global filehandle
+    def printout(self, *s):
+        if self.output_stream is None:
+            return
+        print(*s, file=self.output_stream)
+
+    def runlemke(
+        self,
+        *,
+        verbose=False,
+        lexstats=False,
+        z0=False,
+        output_stream=None,
+    ):
         # z0: printout value of z0
         # flags.maxcount   = 0;
         # flags.bdocupivot = 1;
@@ -479,14 +451,14 @@ def runlemke(self, *, verbose=False, lexstats=False, z0=False, silent=False):
         # flags.binteract  = 0;
         # flags.blexstats  = 0;
 
-        if silent:
-            filehandle = open(outfile, "w")  # noqa: SIM115
+        self.output_stream = output_stream
+
         n = self.n
         self.pivotcount = 1
         # check if d is ok - TBC
         # if (flags.binitabl)
-        printout("After filltableau:")
-        printout(self)
+        self.printout("After filltableau:")
+        self.printout(self)
 
         # z0 enters the basis to obtain lex-feasible solution
         enter = 0
@@ -495,52 +467,75 @@ def runlemke(self, *, verbose=False, lexstats=False, z0=False, silent=False):
         self.negcol(n + 1)
         # if (flags.binitabl)
         if verbose:
-            printout("After negcol:")
-            printout(self)
+            self.printout("After negcol:")
+            self.printout(self)
 
         while True:  # main loop of complementary pivoting
             self.testtablvars()
             if z0:  # printout progress of z0
                 if self.bascobas[0] < n:  # z0 is basic
-                    printout(
+                    self.printout(
                         "step,z0=",
                         self.pivotcount, self.A[self.bascobas[0]][n + 1] / self.determinant
                     )
                 else:
-                    printout("step,z0=", self.pivotcount, 0.0)
+                    self.printout("step,z0=", self.pivotcount, 0.0)
             # if (flags.bdocupivot)
             self.docupivot(leave, enter)
             self.pivot(leave, enter)
             if z0leave:
                 if z0:
-                    printout("step,z0=", self.pivotcount + 1, 0.0)
+                    self.printout("step,z0=", self.pivotcount + 1, 0.0)
                 break
             if verbose:
-                printout(self)
+                self.printout(self)
             enter = self.complement(leave)
             leave, z0leave = self.lexminvar(enter)
             self.pivotcount += 1
 
         # if (flags.binitabl)
-        printout("Final tableau:")
-        printout(self)
+        self.printout("Final tableau:")
+        self.printout(self)
         # if (flags.boutsol)
         self.createsol()
-        printout(self.outsol())
+        self.printout(self.outsol())
         if lexstats:
             self.outstatistics()
     #  end of class tableau
 
 
-def main():
-    processArguments()
-    printout(f"verbose={verbose} lcpfilename={lcpfilename} silent={silent} z0={z0}")
+@click.command(
+    context_settings={"help_option_names": ["-?", "-h", "--help"]},
+)
+@click.option(
+    "-v",
+    "--verbose",
+    is_flag=True,
+    help="Printout intermediate tableaus",
+)
+@click.option(
+    "--z0",
+    is_flag=True,
+    help="Show value of z0 at each step",
+)
+@click.argument(
+    "lcpfilename",
+    type=click.Path(exists=True, readable=True, file_okay=True, dir_okay=False),
+)
+def main(verbose, z0, lcpfilename):
+    """
+    Tool for solving linear complementarity problems using Lemke's algorithm.
+
+    LCPFILENAME is the path to the input file.
+    """
+
+    print(f"verbose={verbose} lcpfilename={lcpfilename} z0={z0}")
     # printout (f"{verbose}= {lcpfilename}= {silent}= {z0}=")
     m = lcp(lcpfilename)
-    printout(m)
-    printout("==================================")
+    print(m)
+    print("==================================")
     tabl = tableau(m)
-    tabl.runlemke(verbose=verbose, z0=z0, silent=silent)
+    tabl.runlemke(verbose=verbose, z0=z0, output_stream=sys.stdout)
 
 
 if __name__ == "__main__":
diff --git a/tests/test_lcp.py b/tests/test_lcp.py
index 4198b85..84baf85 100644
--- a/tests/test_lcp.py
+++ b/tests/test_lcp.py
@@ -26,7 +26,7 @@ def lemke_solver(lcp_instance: lcp) -> list[Fr]:
     """Runs Lemke's algorithm on the given LCP and returns the solution."""
 
     tabl = tableau(lcp_instance)
-    tabl.runlemke(verbose=False, z0=False, silent=False)
+    tabl.runlemke()
     return tabl.solution
 
 
diff --git a/tests/test_lemke.py b/tests/test_lemke.py
new file mode 100644
index 0000000..f52b73d
--- /dev/null
+++ b/tests/test_lemke.py
@@ -0,0 +1,182 @@
+from fractions import Fraction
+
+import pytest
+from click.testing import CliRunner
+
+from lemke.lemke import (
+    lcp,
+    main,
+    tableau,
+)
+
+
+# ---   LCP   -------------------------------------------------------
+@pytest.mark.parametrize("raw, expected", [
+    ("1",    Fraction(1)),
+    ("0.75", Fraction(3, 4)),
+    ("-1/3", Fraction(-1, 3)),
+    ("2/4",  Fraction(1, 2)),
+])
+def test_lcp_valid_file(tmp_path, raw, expected):
+    content = f"n= 1\nM= {raw}\nq= 1\nd= 1\n"
+    file_path = tmp_path / "lcp"
+    file_path.write_text(content)
+
+    m = lcp(str(file_path))
+    assert m.M[0][0] == expected
+
+
+@pytest.mark.parametrize("content", [
+    "M= 1 0 0 1 q= 1 1 d= 1 1\n",           # missing n=
+    "n= 2\nM= 1 0 0 1\nq= 1 1\nd= 1\n",     # wrong number of values
+    "n= 2\nM= 1 0 0 1\nq= 1 1\nx= 1 1\n",   # X= instead of d=
+    "invalid content",
+])
+def test_lcp_invalid_file(tmp_path, content):
+    file_path = tmp_path / "lcp"
+    file_path.write_text(content)
+
+    with pytest.raises(SystemExit) as exc_info:
+        lcp(str(file_path))
+
+    assert exc_info.value.code == 1
+
+
+# ---   TABLEAU INIT   ----------------------------------------------
+def test_tableau_dimensions():
+    n = 3
+    t = tableau(lcp(n))
+    assert t.n == n
+    assert len(t.A) == n
+    assert len(t.A[0]) == n + 2
+
+
+def test_tableau_initial_basis():
+    n = 3
+    t = tableau(lcp(n))
+    for i in range(2 * n + 1):
+        if i <= n:
+            assert t.bascobas[i] >= n, f"Z({i}) should be cobasic"
+        else:
+            assert t.bascobas[i] < n, f"W({i}) should be basic"
+
+
+def test_tableau_bascobas_whichvar_are_inverses():
+    n = 3
+    t = tableau(lcp(n))
+    for i in range(2 * n + 1):
+        assert t.bascobas[t.whichvar[i]] == i
+        assert t.whichvar[t.bascobas[i]] == i
+
+
+# ---   VARTOA   ----------------------------------------------------
+def test_vartoa():
+    n = 3
+    t = tableau(lcp(n))
+    assert t.vartoa(0) == "z0"
+    assert t.vartoa(n + 1) == "w1"
+    assert t.vartoa(n) == f"z{n}"
+    assert t.vartoa(n + n) == f"w{n}"
+
+
+# ---   COMPLEMENT   ------------------------------------------------
+@pytest.mark.parametrize("i", [1, 2, 3])
+def test_complement_pairs(i):
+    n = 3
+    t = tableau(lcp(n))
+    assert t.complement(i) == n + i
+    assert t.complement(n + i) == i
+
+
+def test_complement_z0_fails():
+    t = tableau(lcp(2))
+    with pytest.raises(SystemExit):
+        t.complement(0)
+
+
+# ---   PIVOT   -----------------------------------------------------
+def test_pivot_correctly_swaps_basis_variables():
+    m = lcp(2)
+    m.M = [
+        [Fraction(1), Fraction(0)],
+        [Fraction(0), Fraction(1)],
+    ]
+    m.q = [Fraction(-1), Fraction(-2)]
+    m.d = [Fraction(1), Fraction(1)]
+
+    t = tableau(m)
+
+    enter = 0              # z0 (cobasic in col 0)
+    leave = t.whichvar[0]  # w1 (basic in row 0)
+
+    t.pivot(leave, enter)
+
+    assert t.bascobas[enter] < t.n   # entering variable is basic
+    assert t.bascobas[leave] >= t.n  # leaving variable is cobasic
+
+
+# ---   LEXMINVAR   -------------------------------------------------
+def test_lexminvar_without_positive_entry():
+    t = tableau(lcp(2))
+
+    t.A[0][0] = -1
+    t.A[1][0] = -3
+
+    enter = 0  # z0, cobasic in col 0
+
+    with pytest.raises(SystemExit) as exc_info:
+        t.lexminvar(enter)
+
+    assert exc_info.value.code == 1
+
+
+@pytest.mark.parametrize(
+    "row0, row1",
+    [
+        pytest.param([-1, 0, 0, 0], [3, 0, 0, 0], id="single_candidate"),
+        pytest.param([2, 0, 0, 6], [2, 0, 0, 4], id="multiple_candidates"),
+    ],
+)
+def test_lexminvar_with_positive_entry(row0, row1):
+    t = tableau(lcp(2))
+    t.A[0] = row0
+    t.A[1] = row1
+
+    enter = 0  # z0, cobasic in col 0
+    leave, z0leave = t.lexminvar(enter)
+
+    assert leave == t.whichvar[1]  # w2 leaves
+    assert z0leave is False
+
+
+# ---   CLI   -------------------------------------------------------
+@pytest.mark.parametrize("extra_args", [
+    [],
+    ["--verbose", "--z0"],
+])
+def test_cli_runs_without_error(tmp_path, extra_args):
+    file_path = tmp_path / "test.lcp"
+    file_path.write_text("n= 2\nM= 1 0 0 1\nq= 1 1\nd= 1 1\n")
+
+    runner = CliRunner()
+    result = runner.invoke(main, [str(file_path)] + extra_args)
+
+    assert result.exit_code == 0
+
+
+def test_cli_rejects_missing_file(tmp_path):
+    missing_path = tmp_path / "missing"
+
+    runner = CliRunner()
+    result = runner.invoke(main, [str(missing_path)])
+
+    assert result.exit_code == 2
+    assert "does not exist" in result.output.lower()
+
+
+def test_cli_rejects_directory(tmp_path):
+    runner = CliRunner()
+    result = runner.invoke(main, [str(tmp_path)])
+
+    assert result.exit_code == 2
+    assert "is a directory" in result.output.lower()