Problem solving symultaneous linear equations with integer variables

52 views
Skip to first unread message

Steve Collins

unread,
Aug 4, 2025, 10:06:13 AMAug 4
to sympy
Apologies if this is obvious, but I'm having trouble using solve to find out if solutions exist or not (I don't really care what they are, just whether or not there are any), for simultaneous equations with integer variables.

For the following example, which has no solutions, solve returns [] (no solution), or a solution, depending on the order of the symbols given. That is, it sometimes works but sometimes fails to realize that there is no solution.

For comparison, Mathematica seems to be correct, i.e. Reduce[l == 2 n && 2 h + l - 4 m - 3 == 0, {h, l, n, m}, Integers] correctly returns 'False'. Sympy linsolve, on the other hand, never gives a null result, so it seems to do worse than solve.


Example Python code below. Any advice would be gratefully appreciated!

from sympy import *

(h, l, m, n) = symbols(['h', 'l', 'm', 'n'], integer = True)

# these equations should have no solution as l must be both odd and even
eqnList = [Eq(l, 2*n), Eq(2*h + l - 4*m - 3, 0)]

# however, some orderings of the symbols to solve over fail to give no solution (empty list)
print(solve(eqnList, [l, n, h, m] )) # returns [] - correct
print(solve(eqnList, [h, m, n, l] )) # returns {h: -l/2 + 2*m + 3/2, n: l/2} - should give [] ?

Oscar Benjamin

unread,
Aug 4, 2025, 3:43:35 PMAug 4
to sy...@googlegroups.com
The linsolve function does not check assumptions like integer=True.
The solve function does but only as a final check. Neither of these
really solves integer systems in the way that you would want here.

The diophantine function is for integer equations but does not handle
systems of equations so is also not what you want.

I put some code that illustrates improvements to diophantine in this issue:
https://github.com/sympy/sympy/issues/20682

Using that code you have:

In [9]: diop_new(eqnList, [h, l, m, n])
Out[9]: ∅

In [11]: diop_new([eqnList[0]], [h, l, m, n])
Out[11]: {((t₁, 2⋅t₀, t₂, t₀), (t₀, t₁, t₂), (ℤ, ℤ, ℤ))}

In [13]: diop_new([eqnList[1]], [h, l, m, n])
Out[13]: {((t₀, 2⋅t₀ + 4⋅t₁ + 3, t₀ + t₁, t₂), (t₀, t₁, t₂), (ℤ, ℤ, ℤ))}

That code is not well tested but you might find it useful.

--
Oscar
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/sympy/0c23f1e5-15e8-4345-a961-f8be9728e127n%40googlegroups.com.

Oscar Benjamin

unread,
Aug 4, 2025, 3:51:57 PMAug 4
to sy...@googlegroups.com
Also you can compute the Smith normal form or Hermite normal form
which can be used to check for the existence of solutions to a linear
system of diophantine equations:
https://en.wikipedia.org/wiki/Diophantine_equation#System_of_linear_Diophantine_equations

In [14]: A, b = linear_eq_to_matrix(eqnList, [h, l, m, n])

In [15]: dM = A.to_DM()

In [16]: from sympy.matrices.normalforms import smith_normal_decomp

In [17]: A
Out[17]:
⎡0 1 0 -2⎤
⎢ ⎥
⎣2 1 -4 0 ⎦

In [18]: smith_normal_decomp(A)
Out[18]:
⎛ ⎡0 1 2 -1⎤⎞
⎜ ⎢ ⎥⎟
⎜⎡1 0 0 0⎤ ⎡0 1⎤ ⎢1 -2 0 2 ⎥⎟
⎜⎢ ⎥, ⎢ ⎥, ⎢ ⎥⎟
⎜⎣0 2 0 0⎦ ⎣-1 1⎦ ⎢0 0 1 0 ⎥⎟
⎜ ⎢ ⎥⎟
⎝ ⎣0 0 0 1 ⎦⎠

In [19]: from sympy.matrices.normalforms import hermite_normal_form

In [20]: hermite_normal_form(A)
Out[20]:
⎡2 1⎤
⎢ ⎥
⎣0 1⎦

--
Oscar

Steve Collins

unread,
Aug 5, 2025, 8:51:26 AMAug 5
to sy...@googlegroups.com
Thank you both for taking the time to give me some very helpful and interesting advice!

Reply all
Reply to author
Forward
0 new messages