Get all positive integral solutions for the linear equation

Question

The game I played has a riddle that involves solving the following equation:

x*411 + y*295 + z*161 = 3200 

I don't want to think that I just hit it in sympy , which I really hadn't used until this point:

 >>> from sympy import * >>> x, y, z = symbols('xy 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 to be limited by variables, for example, (if there are no other solutions) [{x: 4, y: 2, z:6}] or [(4, 2, 6)]

Of course, now I can manually replace the two variables in a nested loop or solve it manually (as I did to get the solution above), but I want to know how to get the simps (or another library) to do this for me.

Answer 1

SymPy can solve Diophantine equations , but has no built-in way to create positive solutions. With Sage, you can do this easily: here is a four-line code that generates all non-negative integer solutions to 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() 

Output signal ((4, 2, 6),)

Behind the scenes, integral_points will most likely just run a few loops; although when it does not work, he tries to use Smith's normal form.

I know that you need positive solutions, but (a) it is easy to exclude any null tuple-containing answers; (b) it is also easy to replace x with x-1, etc. before the decision; (c) stick to “non-negative”, simplifies the creation of a polyhedron using the Mixed Integer Linear Programming Module as described above.

According to the documentation, you can also build a polyhedron object directly from the system of inequalities ("Hrep"). This would allow us to say unambiguously x> = 1, etc., but I did not succeed.

With sympy

SymPy Diophantine module output is a parametric solution, for example

 (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 create solutions in a rather efficient way. The sticky point defines the boundaries of the parameters t_0 and t_1. Since this is just an illustration, I looked at the last expression above and connected the limitations of 35200/4816 and 35200/295 directly in the loops below.

 from sympy import * x, y, z = symbols('xy 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 signal [4, 2, 6] .