Multimodular determinant in Python
Computes an integer determinant by evaluating modulo many primes and reconstructing with the Chinese remainder theorem.
- File
multimodular_determinant.py, 102 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/multimodular-determinant/. The authors are listed in prov/BETBZR.jsonl.
"""Determinant of an integer matrix by elimination modulo many primes and the Chinese remainder theorem."""
from math import isqrt
_BASES = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
def is_prime(n: int) -> bool:
"""Miller-Rabin with the first 13 prime bases, deterministic for n < 3.3 * 10^24."""
if n < 2:
return False
for q in _BASES:
if n % q == 0:
return n == q
d, s = n - 1, 0
while d % 2 == 0:
d //= 2
s += 1
for a in _BASES:
x = pow(a, d, n)
if x in (1, n - 1):
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False
return True
def next_prime(n: int) -> int:
"""Return the least prime that is at least n."""
n = max(n, 2)
while not is_prime(n):
n += 1
return n
def hadamard_int(a: list[list[int]]) -> int:
"""Return floor(sqrt(prod of squared row norms)), an integer bound on |det a|."""
prod = 1
for row in a:
prod *= sum(v * v for v in row)
return isqrt(prod)
def det_mod_p(a: list[list[int]], p: int) -> int:
"""Determinant of a square integer matrix modulo the prime p, by Gaussian elimination with row swaps."""
n = len(a)
m = [[v % p for v in row] for row in a]
det = 1
for k in range(n):
piv = next((i for i in range(k, n) if m[i][k] != 0), None)
if piv is None:
return 0
if piv != k:
m[k], m[piv] = m[piv], m[k]
det = -det
det = det * m[k][k] % p
inv = pow(m[k][k], -1, p)
for i in range(k + 1, n):
f = m[i][k] * inv % p
if f:
ri, rk = m[i], m[k]
for j in range(k + 1, n):
ri[j] = (ri[j] - f * rk[j]) % p
return det % p
def multimodular_det_trace(a: list[list[int]], start: int = 1 << 31) -> tuple[int, list[int]]:
"""Return (det a, the primes used), using the primes at least start in increasing order."""
n = len(a)
if any(len(row) != n for row in a):
raise ValueError("the matrix must be square")
if n == 0: # @line L1
return 1, []
h = hadamard_int(a) # @line L2
big_m, r, p = 1, 0, next_prime(start) # @line L3
used: list[int] = []
while big_m <= 2 * h: # @line L4
d = det_mod_p(a, p) # @line L5
r += big_m * ((d - r) * pow(big_m, -1, p) % p) # @line L6
big_m *= p
used.append(p)
p = next_prime(p + 1) # @line L7
if r > big_m // 2: # @line L8
r -= big_m
return r, used # @line L9
def multimodular_det(a: list[list[int]], start: int = 1 << 31) -> int:
"""Return the determinant of the square integer matrix a.
The determinant is computed modulo successive primes p >= start until their
product M exceeds twice the Hadamard bound H, then reconstructed by the
Chinese remainder theorem as the residue of absolute value at most M/2.
"""
return multimodular_det_trace(a, start)[0]
test_multimodular_determinant.py 114 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/multimodular-determinant/. The authors are listed in prov/BETBZR.jsonl.
import itertools
import random
from fractions import Fraction
from math import prod
import pytest
from multimodular_determinant import det_mod_p, hadamard_int, is_prime, multimodular_det, multimodular_det_trace
def leibniz(a):
n = len(a)
total = 0
for s in itertools.permutations(range(n)):
inv = sum(1 for i in range(n) for j in range(i + 1, n) if s[i] > s[j])
total += (-1) ** inv * prod(a[i][s[i]] for i in range(n))
return total
def det_fraction(a):
n = len(a)
m = [[Fraction(v) for v in row] for row in a]
det = Fraction(1)
for k in range(n):
piv = next((i for i in range(k, n) if m[i][k] != 0), None)
if piv is None:
return 0
if piv != k:
m[k], m[piv] = m[piv], m[k]
det = -det
det *= m[k][k]
for i in range(k + 1, n):
f = m[i][k] / m[k][k]
for j in range(k, n):
m[i][j] -= f * m[k][j]
return int(det)
def test_worked_example():
a = [[3, 1, 4], [1, 5, 9], [2, 6, 5]]
assert hadamard_int(a) == 425
assert [det_mod_p(a, p) for p in (7, 11, 13)] == [1, 9, 1]
d, used = multimodular_det_trace(a, start=7)
assert (d, used) == (-90, [7, 11, 13])
assert multimodular_det(a) == -90
def test_edges():
assert multimodular_det([]) == 1
assert multimodular_det([[7]]) == 7
assert multimodular_det([[-5]], start=2) == -5
assert multimodular_det([[0]]) == 0
assert multimodular_det([[0, 1], [1, 0]]) == -1 # zero pivot, needs a swap
assert multimodular_det([[0, 0], [3, 4]]) == 0 # zero row
assert multimodular_det([[0, 2], [0, 4]]) == 0 # zero column
assert multimodular_det([[1, 2], [2, 4]], start=2) == 0 # singular
assert multimodular_det([[2, 0, 0], [0, 3, 0], [0, 0, 35]], start=2) == 210 # det divisible by the first primes
with pytest.raises(ValueError):
multimodular_det([[1, 2]])
def test_is_prime():
small = [n for n in range(2000) if n > 1 and all(n % q for q in range(2, n))]
assert [n for n in range(2000) if is_prime(n)] == small
assert is_prime(2**61 - 1) and not is_prime(2**61 + 1)
def test_random_small_vs_leibniz():
rng = random.Random(7)
for _ in range(600):
n = rng.randrange(0, 6)
a = [[rng.randrange(-6, 7) for _ in range(n)] for _ in range(n)]
if n and rng.random() < 0.2: # make it singular
a[rng.randrange(n)] = [2 * v for v in a[rng.randrange(n)]]
start = rng.choice([2, 3, 50, 1 << 31])
assert multimodular_det(a, start) == leibniz(a)
def test_rank_deficient():
rng = random.Random(8)
for _ in range(100):
n = rng.randrange(2, 8)
r = rng.randrange(1, n)
u = [[rng.randrange(-4, 5) for _ in range(r)] for _ in range(n)]
v = [[rng.randrange(-4, 5) for _ in range(n)] for _ in range(r)]
a = [[sum(u[i][k] * v[k][j] for k in range(r)) for j in range(n)] for i in range(n)]
assert multimodular_det(a, rng.choice([2, 1 << 31])) == 0
def test_large_entries():
rng = random.Random(9)
for _ in range(40):
n = rng.randrange(1, 11)
bits = rng.choice([8, 64, 200])
a = [[rng.randrange(-2**bits, 2**bits) for _ in range(n)] for _ in range(n)]
d, used = multimodular_det_trace(a)
assert d == det_fraction(a)
assert abs(d) <= hadamard_int(a)
h = hadamard_int(a)
# the number of primes is the least s with p_1 ... p_s > 2H, and each prime exceeds 2^31
assert prod(used) > 2 * h and prod(used[:-1]) <= 2 * h
assert len(used) <= (2 * h).bit_length() // 31 + 1
def test_hadamard_tight():
# a Sylvester-Hadamard matrix of order 8 attains the bound: |det| = 8^4
h = [[1]]
for _ in range(3):
h = [r + r for r in h] + [r + [-v for v in r] for r in h]
assert hadamard_int(h) == 8**4
assert abs(multimodular_det(h, start=2)) == 8**4
What the next level needs
- the tests did not run in this build