How would I solve a linear Diophantine congruence in Python?

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 1

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