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].