I was given a challenge where the solution involves solving a series of linear modular equations in 14 variables. The following is a selection of these equations:
3a + 3b + 3c + 3d + 3e + 3f + 3g + h + i + j + k + l + m + n = 15
7a + 9b + 17c + 11d + 6e + 5f + g = 3
13a + 2b + 9c + 8d + 12f + 13g = 17
5a + 2b + 16c + 12d + 5e + 7f + g = 11
6a + 4b + 9c + 6d + 4e + 9f + 6g + h + 7i + 11j + k + 7l + 11m + n = 8
10a + 15b + 13c + 10d + 15e + 13f + 10g + 12h + 18i + 8j + 12k + 18l + 8m + 12n = 18
9a + 12b + 14c + 4d + 9e + 16f + 3g + 7h + 17i + 11j + 14k + 3l + 18m + n = 15
9a + 12b + 16c + 15d + e + 14f + 6g + 11h + 2i + 9j + 12k + 16l + 15m + n = 14
I ended up copying them into this modular equation solver to get a parameterized solution. However, I want to be able to do this automatically in a Python program, without depending on that website.
Preferably, I should be able to do this with an arbitrary group of (linear) equations, not just these specific ones. Part of my solution for the challenge required me to write a few different equations for different scenarios, and swap them out as needed.
I'm aware that SymPy has a Diophantine equation solver that almost does what I want it to do, but in the docs I didn't see a way to get it to enforce a certain modulus (mod 19, mod 23, etc).
The word diophantine is misleading here. What you want to do is linear algebra over Z_p meaning the integers modulo a prime p. SymPy can do this if you use DomainMatrix instead of Matrix. This is how to do it in SymPy 1.12:
In [1]: from sympy import *
In [2]: from sympy.polys.matrices import DomainMatrix
In [3]: a, b, c, d, e, f, g, h, i, j, k, l, m, n = syms = symbols('a:n')
In [4]: eqs = '''
...: 3a + 3b + 3c + 3d + 3e + 3f + 3g + h + i + j + k + l + m + n = 15
...: 7a + 9b + 17c + 11d + 6e + 5f + g = 3
...: 13a + 2b + 9c + 8d + 12f + 13g = 17
...: 5a + 2b + 16c + 12d + 5e + 7f + g = 11
...: 6a + 4b + 9c + 6d + 4e + 9f + 6g + h + 7i + 11j + k + 7l + 11m + n = 8
...: 10a + 15b + 13c + 10d + 15e + 13f + 10g + 12h + 18i + 8j + 12k + 18l + 8m + 12n = 18
...: 9a + 12b + 14c + 4d + 9e + 16f + 3g + 7h + 17i + 11j + 14k + 3l + 18m + n = 15
...: 9a + 12b + 16c + 15d + e + 14f + 6g + 11h + 2i + 9j + 12k + 16l + 15m + n = 14
...: '''
In [5]: eqs = [parse_expr(eq, transformations='all').subs(I, i) for eq in eqs.strip().splitlines()]
In [6]: M, b = linear_eq_to_matrix(eqs, syms)
In [7]: M
Out[7]:
⎡3 3 3 3 3 3 3 1 1 1 1 1 1 1 ⎤
⎢ ⎥
⎢7 9 17 11 6 5 1 0 0 0 0 0 0 0 ⎥
⎢ ⎥
⎢13 2 9 8 0 12 13 0 0 0 0 0 0 0 ⎥
⎢ ⎥
⎢5 2 16 12 5 7 1 0 0 0 0 0 0 0 ⎥
⎢ ⎥
⎢6 4 9 6 4 9 6 1 18 0 1 7 11 1 ⎥
⎢ ⎥
⎢10 15 13 10 15 13 10 12 26 0 12 18 8 12⎥
⎢ ⎥
⎢9 12 14 4 9 16 3 7 28 0 14 3 18 1 ⎥
⎢ ⎥
⎣9 12 16 15 1 14 6 11 11 0 12 16 15 1 ⎦
In [8]: b
Out[8]:
⎡15⎤
⎢ ⎥
⎢3 ⎥
⎢ ⎥
⎢17⎥
⎢ ⎥
⎢11⎥
⎢ ⎥
⎢8 ⎥
⎢ ⎥
⎢18⎥
⎢ ⎥
⎢15⎥
⎢ ⎥
⎣14⎦
In [9]: Z_17 = GF(17, symmetric=False)
In [12]: dM = DomainMatrix.from_Matrix(M).convert_to(Z_17).to_dense()
In [13]: db = DomainMatrix.from_Matrix(b).convert_to(Z_17).to_dense()
In [14]: dM
Out[14]: DomainMatrix([[3 mod 17, 3 mod 17, 3 mod 17, 3 mod 17, 3 mod 17, 3 mod 17, 3 mod 17, 1 mod 17, 1 mod 17, 1 mod 17, 1 mod 17, 1 mod 17, 1 mod 17, 1 mod 17], [7 mod 17, 9 mod 17, 0 mod 17, 11 mod 17, 6 mod 17, 5 mod 17, 1 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17], [13 mod 17, 2 mod 17, 9 mod 17, 8 mod 17, 0 mod 17, 12 mod 17, 13 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17], [5 mod 17, 2 mod 17, 16 mod 17, 12 mod 17, 5 mod 17, 7 mod 17, 1 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17, 0 mod 17], [6 mod 17, 4 mod 17, 9 mod 17, 6 mod 17, 4 mod 17, 9 mod 17, 6 mod 17, 1 mod 17, 1 mod 17, 0 mod 17, 1 mod 17, 7 mod 17, 11 mod 17, 1 mod 17], [10 mod 17, 15 mod 17, 13 mod 17, 10 mod 17, 15 mod 17, 13 mod 17, 10 mod 17, 12 mod 17, 9 mod 17, 0 mod 17, 12 mod 17, 1 mod 17, 8 mod 17, 12 mod 17], [9 mod 17, 12 mod 17, 14 mod 17, 4 mod 17, 9 mod 17, 16 mod 17, 3 mod 17, 7 mod 17, 11 mod 17, 0 mod 17, 14 mod 17, 3 mod 17, 1 mod 17, 1 mod 17], [9 mod 17, 12 mod 17, 16 mod 17, 15 mod 17, 1 mod 17, 14 mod 17, 6 mod 17, 11 mod 17, 11 mod 17, 0 mod 17, 12 mod 17, 16 mod 17, 15 mod 17, 1 mod 17]], (8, 14), GF(17))
In [15]: db
Out[15]: DomainMatrix([[15 mod 17], [3 mod 17], [0 mod 17], [11 mod 17], [8 mod 17], [1 mod 17], [15 mod 17], [14 mod 17]], (8, 1), GF(17))
Now you have dM and db as matrices over Z_17 and you can use whatever matrix operations you want and convert back to Matrix with to_Matrix:
In [16]: dM.hstack(db).rref()[0].to_Matrix() # RREF of augmented matrix
Out[16]:
⎡1 0 0 0 0 0 0 0 11 11 5 2 2 6 0 ⎤
⎢ ⎥
⎢0 1 0 0 0 0 0 0 14 2 11 12 14 13 6 ⎥
⎢ ⎥
⎢0 0 1 0 0 0 0 0 3 7 10 11 15 0 5 ⎥
⎢ ⎥
⎢0 0 0 1 0 0 0 0 2 2 1 16 16 13 5 ⎥
⎢ ⎥
⎢0 0 0 0 1 0 0 0 13 13 16 7 14 0 15⎥
⎢ ⎥
⎢0 0 0 0 0 1 0 0 16 10 5 11 15 11 7 ⎥
⎢ ⎥
⎢0 0 0 0 0 0 1 0 8 10 6 13 1 0 7 ⎥
⎢ ⎥
⎣0 0 0 0 0 0 0 1 4 6 9 6 8 8 16⎦
In [17]: dM.nullspace().to_Matrix() # nullspace of M
Out[17]:
⎡15 16 1 12 10 11 14 7 11 0 0 0 0 0 ⎤
⎢ ⎥
⎢15 12 8 12 10 9 9 2 0 11 0 0 0 0 ⎥
⎢ ⎥
⎢13 15 9 6 11 13 2 3 0 0 11 0 0 0 ⎥
⎢ ⎥
⎢12 4 15 11 8 15 10 2 0 0 0 11 0 0 ⎥
⎢ ⎥
⎢12 16 5 11 16 5 6 14 0 0 0 0 11 0 ⎥
⎢ ⎥
⎣2 10 0 10 0 15 0 14 0 0 0 0 0 11⎦
You can use these to construct your parametric solution:
In [30]: p = dM.hstack(db).rref()[0].to_Matrix()[:,-1]
In [31]: p = Matrix.vstack(p, p.zeros(6, 1))
In [32]: p.transpose()
Out[32]: [0 6 5 5 15 7 7 16 0 0 0 0 0 0]
In [33]: (M*p - b).applyfunc(lambda e: e % 17).is_zero_matrix
Out[33]: True
In [34]: N = dM.nullspace().to_Matrix()
In [35]: N = dM.nullspace().transpose().to_Matrix()
In [36]: (M*N).applyfunc(lambda e: e % 17).is_zero_matrix
Out[36]: True
The parametric solution is generated from the particular solution and nullspace:
In [38]: sol = p + N*Matrix(symbols('alpha:6'))
In [39]: sol
Out[39]:
⎡ 15⋅α₀ + 15⋅α₁ + 13⋅α₂ + 12⋅α₃ + 12⋅α₄ + 2⋅α₅ ⎤
⎢ ⎥
⎢16⋅α₀ + 12⋅α₁ + 15⋅α₂ + 4⋅α₃ + 16⋅α₄ + 10⋅α₅ + 6⎥
⎢ ⎥
⎢ α₀ + 8⋅α₁ + 9⋅α₂ + 15⋅α₃ + 5⋅α₄ + 5 ⎥
⎢ ⎥
⎢12⋅α₀ + 12⋅α₁ + 6⋅α₂ + 11⋅α₃ + 11⋅α₄ + 10⋅α₅ + 5⎥
⎢ ⎥
⎢ 10⋅α₀ + 10⋅α₁ + 11⋅α₂ + 8⋅α₃ + 16⋅α₄ + 15 ⎥
⎢ ⎥
⎢11⋅α₀ + 9⋅α₁ + 13⋅α₂ + 15⋅α₃ + 5⋅α₄ + 15⋅α₅ + 7 ⎥
⎢ ⎥
⎢ 14⋅α₀ + 9⋅α₁ + 2⋅α₂ + 10⋅α₃ + 6⋅α₄ + 7 ⎥
⎢ ⎥
⎢ 7⋅α₀ + 2⋅α₁ + 3⋅α₂ + 2⋅α₃ + 14⋅α₄ + 14⋅α₅ + 16 ⎥
⎢ ⎥
⎢ 11⋅α₀ ⎥
⎢ ⎥
⎢ 11⋅α₁ ⎥
⎢ ⎥
⎢ 11⋅α₂ ⎥
⎢ ⎥
⎢ 11⋅α₃ ⎥
⎢ ⎥
⎢ 11⋅α₄ ⎥
⎢ ⎥
⎣ 11⋅α₅ ⎦
Let's check that:
In [47]: (M*sol)._rep.convert_to(Z_17[symbols('alpha:6')]).to_Matrix() == b.applyfunc(lambda e: e % 17)
Out[47]: True