From f34598abfde04cb4fcc45d31f4435638448cc4d1 Mon Sep 17 00:00:00 2001 From: Stanislav Date: Fri, 16 Oct 2026 12:04:06 +0400 Subject: [PATCH 1/2] Add the Modulus option to RowReduce, NullSpace, MatrixRank, Inverse and LinearSolve With Modulus -> p for a prime p, these compute over the integers mod p with plain Gaussian elimination on Python ints. Entries may be integers or rationals whose denominator is invertible mod p. Modulus -> 0 keeps ordinary arithmetic; a modulus that is not a prime, or a symbolic entry, leaves the expression unevaluated. Co-Authored-By: Claude Opus 5.5 --- mathics/builtin/numbers/linalg.py | 108 +++++++++++++++-- mathics/eval/numbers/linalg_modular.py | 153 +++++++++++++++++++++++++ test/builtin/numbers/test_linalg.py | 110 ++++++++++++++++++ 3 files changed, 361 insertions(+), 10 deletions(-) create mode 100644 mathics/eval/numbers/linalg_modular.py 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..ac555d7471 --- /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. + + As with the exact NullSpace, there is one vector for each column + without a pivot, in increasing order of that column; it has 1 there + and 0 in the other columns without a pivot. + """ + reduced, pivots = rref_mod(rows, p) + ncols = len(rows[0]) + pivot_set = set(pivots) + basis = [] + for free in 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..817d4f4459 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, 1, 0}, {6, 0, 1}}"), + ("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[" + " 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[" + " m = RandomInteger[{0, 6}, {5, 5}];" + " inv = Quiet[Inverse[m, Modulus -> 7]];" + " Head[inv] === Inverse || Mod[m . inv, 7] === IdentityMatrix[5]," + " {k, 30}]", + "True", + ) From 4e65388f006bcde7c305dc9a5c90871987d712e2 Mon Sep 17 00:00:00 2001 From: Stanislav Date: Fri, 16 Oct 2026 12:22:26 +0400 Subject: [PATCH 2/2] List modular NullSpace vectors as WMA does; scope test variables Strickland's LambdaAlgebra.m, run in WMA, shows that NullSpace with Modulus gives the vectors in decreasing order of their free column, so do the same. The random tests assigned x, m and p in the shared test session, which broke later tests that use x as a symbol; use Module. Co-Authored-By: Claude Opus 5.5 --- mathics/eval/numbers/linalg_modular.py | 8 ++++---- test/builtin/numbers/test_linalg.py | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/mathics/eval/numbers/linalg_modular.py b/mathics/eval/numbers/linalg_modular.py index ac555d7471..e915bb87d0 100644 --- a/mathics/eval/numbers/linalg_modular.py +++ b/mathics/eval/numbers/linalg_modular.py @@ -107,15 +107,15 @@ def rref_mod(rows: ModMatrix, p: int) -> Tuple[ModMatrix, List[int]]: def nullspace_mod(rows: ModMatrix, p: int) -> ModMatrix: """Return a basis of the null space mod p. - As with the exact NullSpace, there is one vector for each column - without a pivot, in increasing order of that column; it has 1 there - and 0 in the other columns without a pivot. + 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 range(ncols): + for free in reversed(range(ncols)): if free in pivot_set: continue vector = [0] * ncols diff --git a/test/builtin/numbers/test_linalg.py b/test/builtin/numbers/test_linalg.py index 817d4f4459..cdba04fd08 100644 --- a/test/builtin/numbers/test_linalg.py +++ b/test/builtin/numbers/test_linalg.py @@ -297,7 +297,7 @@ def test_linalg(str_expr, msgs, str_expected, fail_msg): ("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, 1, 0}, {6, 0, 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]", @@ -369,7 +369,7 @@ def test_linalg_modulus_random(): """Check the modular results against their defining properties.""" check_evaluation( "SeedRandom[42]; " - "And @@ Table[" + "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];" @@ -379,16 +379,16 @@ def test_linalg_modulus_random(): " 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]," + " Count[rr, Except[{0 ..}]] == rk]]," " {k, 30}]", "True", ) check_evaluation( "SeedRandom[7]; " - "And @@ Table[" + "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]," + " Head[inv] === Inverse || Mod[m . inv, 7] === IdentityMatrix[5]]," " {k, 30}]", "True", )