diff --git a/mathics/builtin/numbers/linalg.py b/mathics/builtin/numbers/linalg.py
index e6296bba1d..9d71828383 100644
--- a/mathics/builtin/numbers/linalg.py
+++ b/mathics/builtin/numbers/linalg.py
@@ -19,6 +19,16 @@
 from mathics.core.list import ListExpression
 from mathics.core.number import is_inexact
 from mathics.core.symbols import SymbolList
+from mathics.eval.numbers.linalg_modular import (
+    get_prime_modulus,
+    inverse_mod,
+    linear_solve_mod,
+    matrix_mod,
+    nullspace_mod,
+    rref_mod,
+    to_mathics_matrix,
+    vector_mod,
+)
 
 
 class DesignMatrix(Builtin):
@@ -308,16 +318,23 @@ class Inverse(Builtin):
      : The matrix {{1, 0}, {0, 0}} is singular.
      = Inverse[{{1, 0}, {0, 0}}]
 
+    With 'Modulus' -> $p$, the inverse is computed modulo the prime $p$:
+    >> Inverse[{{12, -2}, {1, 3}}, Modulus -> 5]
+     = {{1, 4}, {3, 4}}
     """
 
     messages = {
         "sing": "The matrix `1` is singular.",
         "matsq": "Argument `1` at position 1 is not a non-empty square matrix.",
     }
+    options = {"Modulus": "0"}
     summary_text = "inverse matrix"
 
-    def eval(self, m, evaluation: Evaluation):
-        "Inverse[m_List]"
+    def eval(self, m, evaluation: Evaluation, options: dict):
+        "Inverse[m_List, OptionsPattern[Inverse]]"
+        p = get_prime_modulus(self.get_option(options, "Modulus", evaluation))
+        if p is None:
+            return None
         rows = m.elements
         nrows = len(rows)
         for row in rows:
@@ -331,6 +348,16 @@ def eval(self, m, evaluation: Evaluation):
                 evaluation.message("Inverse", "matsq", m)
                 return None
 
+        if p:
+            mod_matrix = matrix_mod(m, p)
+            if mod_matrix is None:
+                return None
+            mod_inverse = inverse_mod(mod_matrix, p)
+            if mod_inverse is None:
+                evaluation.message("Inverse", "sing", m)
+                return None
+            return to_mathics_matrix(mod_inverse)
+
         matrix = to_sympy_matrix(m)
         det = matrix.det()
         if det == 0:
@@ -492,6 +519,10 @@ class LinearSolve(Builtin):
     >> LinearSolve[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {1, -2, 3}]
      : Linear equation encountered that has no solution.
      = LinearSolve[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}, {1, -2, 3}]
+
+    With 'Modulus' -> $p$, the system is solved modulo the prime $p$:
+    >> LinearSolve[{{1, 1, 1}, {1, 2, 3}, {1, 4, 9}}, {1, 2, 3}, Modulus -> 47]
+     = {23, 2, 23}
     """
 
     messages = {
@@ -502,11 +533,15 @@ class LinearSolve(Builtin):
         "nosol": "Linear equation encountered that has no solution.",
         "matrix": "Argument `1` at position `2` is not a non-empty rectangular matrix.",
     }
+    options = {"Modulus": "0"}
     summary_text = "solves linear systems in matrix form"
 
-    def eval(self, m, b, evaluation: Evaluation):
-        "LinearSolve[m_, b_]"
+    def eval(self, m, b, evaluation: Evaluation, options: dict):
+        "LinearSolve[m_, b_, OptionsPattern[LinearSolve]]"
 
+        p = get_prime_modulus(self.get_option(options, "Modulus", evaluation))
+        if p is None:
+            return None
         matrix = matrix_data(m)
         if matrix is None:
             evaluation.message("LinearSolve", "matrix", m, 1)
@@ -522,6 +557,17 @@ def eval(self, m, b, evaluation: Evaluation):
                 evaluation.message("LinearSolve", "matrix", b, 2)
                 return
 
+        if p:
+            mod_matrix = matrix_mod(m, p)
+            mod_b = vector_mod(b, p)
+            if mod_matrix is None or mod_b is None:
+                return None
+            mod_solution = linear_solve_mod(mod_matrix, mod_b, p)
+            if mod_solution is None:
+                evaluation.message("LinearSolve", "nosol")
+                return None
+            return ListExpression(*(Integer(x) for x in mod_solution))
+
         system = [mm + [v.to_sympy()] for mm, v in zip(matrix, b.elements)]
         system = to_sympy_matrix(system)
         if system is None:
@@ -642,15 +688,29 @@ class MatrixRank(Builtin):
      = 3
     >> MatrixRank[{{a, b}, {3 a, 3 b}}]
      = 1
+
+    The rank depends on the modulus used:
+    >> MatrixRank[{{1, 0, 4}, {2, 0, 3}, {2, 1, 2}}, Modulus -> 5]
+     = 2
     """
 
     messages = {
         "matrix": "Argument `1` at position `2` is not a non-empty rectangular matrix."
     }
+    options = {"Modulus": "0"}
     summary_text = "rank of a matrix"
 
-    def eval(self, m, evaluation: Evaluation):
-        "MatrixRank[m_]"
+    def eval(self, m, evaluation: Evaluation, options: dict):
+        "MatrixRank[m_, OptionsPattern[MatrixRank]]"
+
+        p = get_prime_modulus(self.get_option(options, "Modulus", evaluation))
+        if p is None:
+            return None
+        if p:
+            mod_matrix = matrix_mod(m, p)
+            if mod_matrix is None:
+                return None
+            return Integer(len(rref_mod(mod_matrix, p)[1]))
 
         matrix = to_sympy_matrix(m)
         if matrix is None:
@@ -678,15 +738,29 @@ class NullSpace(Builtin):
      = {}
     >> MatrixRank[A]
      = 3
+
+    With 'Modulus' -> $p$, the null space is computed modulo the prime $p$:
+    >> NullSpace[A, Modulus -> 2]
+     = {{1, 1, 1}}
     """
 
     messages = {
         "matrix": "Argument `1` at position `2` is not a non-empty rectangular matrix."
     }
+    options = {"Modulus": "0"}
     summary_text = "generators for the null space of a matrix"
 
-    def eval(self, m, evaluation: Evaluation):
-        "NullSpace[m_]"
+    def eval(self, m, evaluation: Evaluation, options: dict):
+        "NullSpace[m_, OptionsPattern[NullSpace]]"
+
+        p = get_prime_modulus(self.get_option(options, "Modulus", evaluation))
+        if p is None:
+            return None
+        if p:
+            mod_matrix = matrix_mod(m, p)
+            if mod_matrix is None:
+                return None
+            return to_mathics_matrix(nullspace_mod(mod_matrix, p))
 
         matrix = to_sympy_matrix(m)
         if matrix is None:
@@ -789,15 +863,29 @@ class RowReduce(Builtin):
      . 0   1   2
      .
      . 0   0   0
+
+    With 'Modulus' -> $p$, the reduction is done modulo the prime $p$:
+    >> RowReduce[{{1, 0, 4}, {2, 0, 3}, {2, 1, 2}}, Modulus -> 5]
+     = {{1, 0, 4}, {0, 1, 4}, {0, 0, 0}}
     """
 
     messages = {
         "matrix": "Argument `1` at position `2` is not a non-empty rectangular matrix."
     }
+    options = {"Modulus": "0"}
     summary_text = "matrix reduced row-echelon form"
 
-    def eval(self, m, evaluation: Evaluation):
-        "RowReduce[m_]"
+    def eval(self, m, evaluation: Evaluation, options: dict):
+        "RowReduce[m_, OptionsPattern[RowReduce]]"
+
+        p = get_prime_modulus(self.get_option(options, "Modulus", evaluation))
+        if p is None:
+            return None
+        if p:
+            mod_matrix = matrix_mod(m, p)
+            if mod_matrix is None:
+                return None
+            return to_mathics_matrix(rref_mod(mod_matrix, p)[0])
 
         matrix = to_sympy_matrix(m)
         if matrix is None:
diff --git a/mathics/eval/numbers/linalg_modular.py b/mathics/eval/numbers/linalg_modular.py
new file mode 100644
index 0000000000..e915bb87d0
--- /dev/null
+++ b/mathics/eval/numbers/linalg_modular.py
@@ -0,0 +1,153 @@
+"""
+Linear algebra over the integers modulo a prime, used for the ``Modulus``
+option of RowReduce, NullSpace, MatrixRank, Inverse and LinearSolve.
+
+Matrices are plain Python lists of lists of ``int`` in the range 0..p-1.
+"""
+
+from typing import List, Optional, Tuple
+
+import sympy
+
+from mathics.core.atoms import Integer, Rational
+from mathics.core.element import BaseElement
+from mathics.core.list import ListExpression
+from mathics.core.symbols import SymbolList
+
+ModMatrix = List[List[int]]
+
+
+def get_prime_modulus(modulus: Optional[BaseElement]) -> Optional[int]:
+    """Return the value of a ``Modulus`` option.
+
+    0 means ordinary arithmetic. None means the value is not a prime
+    or 0, and is not supported.
+    """
+    if not isinstance(modulus, Integer):
+        return None
+    p = modulus.value
+    if p == 0 or (p > 1 and sympy.isprime(p)):
+        return p
+    return None
+
+
+def element_mod(element: BaseElement, p: int) -> Optional[int]:
+    """Reduce an Integer or a Rational modulo p, or return None."""
+    if isinstance(element, Integer):
+        return element.value % p
+    if isinstance(element, Rational):
+        num, den = (int(x) for x in element.value.as_numer_denom())
+        if den % p == 0:
+            return None
+        return num * pow(den, -1, p) % p
+    return None
+
+
+def vector_mod(vector: BaseElement, p: int) -> Optional[List[int]]:
+    """Convert a List of rational numbers into a list of residues mod p."""
+    if not vector.has_form(SymbolList, None):
+        return None
+    result = []
+    for element in vector.elements:
+        value = element_mod(element, p)
+        if value is None:
+            return None
+        result.append(value)
+    return result
+
+
+def matrix_mod(matrix: BaseElement, p: int) -> Optional[ModMatrix]:
+    """Convert a non-empty rectangular matrix of rational numbers into
+    residues mod p, or return None."""
+    if not matrix.has_form(SymbolList, None) or len(matrix.elements) == 0:
+        return None
+    rows = []
+    for row in matrix.elements:
+        values = vector_mod(row, p)
+        if not values:
+            return None
+        rows.append(values)
+    if any(len(row) != len(rows[0]) for row in rows):
+        return None
+    return rows
+
+
+def to_mathics_matrix(rows: ModMatrix) -> ListExpression:
+    return ListExpression(*(ListExpression(*(Integer(x) for x in row)) for row in rows))
+
+
+def rref_mod(rows: ModMatrix, p: int) -> Tuple[ModMatrix, List[int]]:
+    """Return the reduced row echelon form of ``rows`` mod p and the
+    list of its pivot columns. ``rows`` is not modified."""
+    rows = [list(row) for row in rows]
+    nrows = len(rows)
+    ncols = len(rows[0]) if rows else 0
+    pivots: List[int] = []
+    r = 0
+    for c in range(ncols):
+        if r == nrows:
+            break
+        pivot = next((i for i in range(r, nrows) if rows[i][c]), None)
+        if pivot is None:
+            continue
+        rows[r], rows[pivot] = rows[pivot], rows[r]
+        inverse = pow(rows[r][c], -1, p)
+        if inverse != 1:
+            rows[r] = [x * inverse % p for x in rows[r]]
+        pivot_row = rows[r]
+        for i in range(nrows):
+            factor = rows[i][c]
+            if i != r and factor:
+                rows[i] = [(x - factor * y) % p for x, y in zip(rows[i], pivot_row)]
+        pivots.append(c)
+        r += 1
+    return rows, pivots
+
+
+def nullspace_mod(rows: ModMatrix, p: int) -> ModMatrix:
+    """Return a basis of the null space mod p.
+
+    There is one vector for each column without a pivot; it has 1 there
+    and 0 in the other columns without a pivot. As in WMA, the vectors
+    come in decreasing order of that column.
+    """
+    reduced, pivots = rref_mod(rows, p)
+    ncols = len(rows[0])
+    pivot_set = set(pivots)
+    basis = []
+    for free in reversed(range(ncols)):
+        if free in pivot_set:
+            continue
+        vector = [0] * ncols
+        vector[free] = 1
+        for i, c in enumerate(pivots):
+            vector[c] = -reduced[i][free] % p
+        basis.append(vector)
+    return basis
+
+
+def inverse_mod(rows: ModMatrix, p: int) -> Optional[ModMatrix]:
+    """Return the inverse of a square matrix mod p, or None if it is
+    singular mod p."""
+    n = len(rows)
+    augmented = [
+        row + [1 if i == j else 0 for j in range(n)] for i, row in enumerate(rows)
+    ]
+    reduced, pivots = rref_mod(augmented, p)
+    if pivots[:n] != list(range(n)):
+        return None
+    return [row[n:] for row in reduced]
+
+
+def linear_solve_mod(rows: ModMatrix, b: List[int], p: int) -> Optional[List[int]]:
+    """Return one solution x of ``rows . x == b`` mod p, with 0 for the
+    free variables, or None if there is no solution."""
+    ncols = len(rows[0])
+    augmented = [row + [v] for row, v in zip(rows, b)]
+    reduced, pivots = rref_mod(augmented, p)
+    if pivots and pivots[-1] == ncols:
+        return None
+    solution = [0] * ncols
+    for i, c in enumerate(pivots):
+        solution[c] = reduced[i][ncols]
+    return solution
diff --git a/test/builtin/numbers/test_linalg.py b/test/builtin/numbers/test_linalg.py
index 8f97ccdc61..cdba04fd08 100644
--- a/test/builtin/numbers/test_linalg.py
+++ b/test/builtin/numbers/test_linalg.py
@@ -2,6 +2,7 @@
 """
 Unit tests for mathics.builtins.numbers.linalg
 """
+
 from test.helper import check_arg_counts, check_evaluation
 
 import pytest
@@ -282,3 +283,112 @@ def test_linalg(str_expr, msgs, str_expected, fail_msg):
         failure_message=fail_msg,
         expected_messages=msgs,
     )
+
+
+@pytest.mark.parametrize(
+    ("str_expr", "msgs", "str_expected"),
+    [
+        (
+            "RowReduce[{{1, 0, 4}, {2, 0, 3}, {2, 1, 2}}, Modulus -> 5]",
+            None,
+            "{{1, 0, 4}, {0, 1, 4}, {0, 0, 0}}",
+        ),
+        ("RowReduce[{{1, 1}, {1, 1}}, Modulus -> 2]", None, "{{1, 1}, {0, 0}}"),
+        ("RowReduce[{{1/2, 1}, {1, 3}}, Modulus -> 3]", None, "{{1, 0}, {0, 1}}"),
+        ("RowReduce[{{-1, 2}, {3, -4}}, Modulus -> 0]", None, "{{1, 0}, {0, 1}}"),
+        ("NullSpace[{{1, 1}, {1, 1}}, Modulus -> 2]", None, "{{1, 1}}"),
+        ("NullSpace[{{1, 1, 1}}, Modulus -> 7]", None, "{{6, 0, 1}, {6, 1, 0}}"),
+        ("NullSpace[{{1, 0}, {0, 1}}, Modulus -> 3]", None, "{}"),
+        (
+            "NullSpace[{{1, 0, 4}, {2, 0, 3}, {2, 1, 2}}, Modulus -> 5]",
+            None,
+            "{{1, 1, 1}}",
+        ),
+        ("MatrixRank[{{1, 0, 4}, {2, 0, 3}, {2, 1, 2}}, Modulus -> 5]", None, "2"),
+        ("MatrixRank[{{2, 4}, {6, 8}}, Modulus -> 2]", None, "0"),
+        ("Inverse[{{1, 1}, {0, 1}}, Modulus -> 2]", None, "{{1, 1}, {0, 1}}"),
+        ("Inverse[{{12, -2}, {1, 3}}, Modulus -> 5]", None, "{{1, 4}, {3, 4}}"),
+        (
+            "Head[Inverse[{{1, 1}, {1, 1}}, Modulus -> 2]]",
+            ("The matrix {{1, 1}, {1, 1}} is singular.",),
+            "Inverse",
+        ),
+        (
+            "Head[Inverse[{{1, 2}, {3, 4}}, Modulus -> 2]]",
+            ("The matrix {{1, 2}, {3, 4}} is singular.",),
+            "Inverse",
+        ),
+        ("LinearSolve[{{1, 1}, {0, 1}}, {1, 0}, Modulus -> 2]", None, "{1, 0}"),
+        (
+            "LinearSolve[{{1, 1, 1}, {1, 2, 3}, {1, 4, 9}}, {1, 2, 3}, Modulus -> 47]",
+            None,
+            "{23, 2, 23}",
+        ),
+        ("LinearSolve[{{1, 1}, {1, 1}}, {1, 1}, Modulus -> 2]", None, "{1, 0}"),
+        (
+            "Head[LinearSolve[{{1, 1}, {1, 1}}, {1, 0}, Modulus -> 2]]",
+            ("Linear equation encountered that has no solution.",),
+            "LinearSolve",
+        ),
+        # Left unevaluated: a modulus that is not a prime, a symbolic entry,
+        # and a denominator that is not invertible.
+        (
+            "Head[NullSpace[{{1, 2}, {3, 2}}, Modulus -> 6]]",
+            None,
+            "NullSpace",
+        ),
+        (
+            "Head[MatrixRank[{{1, 2}, {3, 2}}, Modulus -> x]]",
+            None,
+            "MatrixRank",
+        ),
+        (
+            "Head[RowReduce[{{a, 1}, {1, 3}}, Modulus -> 3]]",
+            None,
+            "RowReduce",
+        ),
+        (
+            "Head[RowReduce[{{1/3, 1}, {1, 3}}, Modulus -> 3]]",
+            None,
+            "RowReduce",
+        ),
+    ],
+)
+def test_linalg_modulus(str_expr, msgs, str_expected):
+    check_evaluation(
+        str_expr,
+        str_expected,
+        to_string_expr=True,
+        to_string_expected=True,
+        hold_expected=True,
+        expected_messages=msgs,
+    )
+
+
+def test_linalg_modulus_random():
+    """Check the modular results against their defining properties."""
+    check_evaluation(
+        "SeedRandom[42]; "
+        "And @@ Table[Module[{m, p, ns, rr, rk, x},"
+        "  m = RandomInteger[{-5, 5}, {6, 8}]; p = {2, 3, 7}[[Mod[k, 3] + 1]];"
+        "  ns = NullSpace[m, Modulus -> p];"
+        "  rr = RowReduce[m, Modulus -> p];"
+        "  rk = MatrixRank[m, Modulus -> p];"
+        "  x = LinearSolve[m, m . Range[8], Modulus -> p];"
+        "  And["
+        "    Length[ns] + rk == 8,"
+        "    ns === {} || Union[Flatten[Mod[m . Transpose[ns], p]]] === {0},"
+        "    Mod[m . x - m . Range[8], p] === Table[0, 6],"
+        "    Count[rr, Except[{0 ..}]] == rk]],"
+        "  {k, 30}]",
+        "True",
+    )
+    check_evaluation(
+        "SeedRandom[7]; "
+        "And @@ Table[Module[{m, inv},"
+        "  m = RandomInteger[{0, 6}, {5, 5}];"
+        "  inv = Quiet[Inverse[m, Modulus -> 7]];"
+        "  Head[inv] === Inverse || Mod[m . inv, 7] === IdentityMatrix[5]],"
+        "  {k, 30}]",
+        "True",
+    )