A game I played has a riddle that involves solving the following equation:
x*411 + y*295 + z*161 = 3200
Not wanting to think I just slapped it into sympy, which I haven’t really used up to that point:
>>> from sympy import *
>>> x, y, z = symbols('x y z', integer=True, positive=True)
>>> solve(x*411 + y*295 + z*161 - 3200, [x, y, z])
[{x: -295*y/411 - 161*z/411 + 3200/411}]
Hmm, this only gave me a dependent solution, but I want all possible solutions in the domain I constrained the variables to, e.g. (assuming there are no other solutions) [{x: 4, y: 2, z:6}] or [(4, 2, 6)]
Of course I could now manually substitute two variables in a nested loop, or solve it by hand (as I did to get the solution above), but I want to know how to get sympy (or another library) to do it for me.
SymPy can solve Diophantine equations but doesn't have a built-in way to generate positive solutions. With Sage one can do that easily: here is four-line code that generates all nonnegative integer solutions of your equation.
p = MixedIntegerLinearProgram()
w = p.new_variable(integer=True, nonnegative=True)
p.add_constraint(411*w[0] + 295*w[1] + 161*w[2] == 3200)
p.polyhedron().integral_points()
The output is ((4, 2, 6),)
Behind the scenes, integral_points will most likely just run a multiple loop; although when that doesn't seem to work it tries to use Smith normal form.
I know you wanted positive solutions, but (a) it's easy to exclude any zero-containing tuples from the answer; (b) it's also easy to replace x by x-1, etc, prior to solving; (c) sticking to "nonnegative" makes it easy to create a polyhedron using Mixed Integer Linear Programming module
as above.
According to documentation one can also build a Polyhedron object directly from a system of inequalities ("Hrep"). This would allow one to explicitly say x >= 1, etc, but I haven't succeeded at this route.
With SymPy
The output of SymPy's Diophantine module is a parametric solution, like
(t_0, 2627*t_0 + 161*t_1 - 19200, -4816*t_0 - 295*t_1 + 35200)
in your example. This can be used in a loop to generate solutions in a pretty efficient way. The sticky point is finding bounds for parameters t_0 and t_1. Since this is just an illustration, I looked at the last expression above and plugged the limits 35200/4816 and 35200/295 directly in the loops below.
from sympy import *
x, y, z = symbols('x y z')
[s] = diophantine(x*411 + y*295 + z*161 - 3200)
print(s)
t_0, t_1 = s[2].free_symbols
for t0 in range(int(35200/4816)+1):
for t1 in range(int(35200/295)+1):
sol = [expr.subs({t_0: t0, t_1: t1}) for expr in s]
if min(sol) > 0:
print(sol)
The output is [4, 2, 6].
