Question: 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).

Answer: 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