title: Get all positive integral solutions for a linear equation
body: <p>A game I played has a riddle that involves solving the following equation:</p>
  
  <pre><code>x*411 + y*295 + z*161 = 3200
  </code></pre>
  
  <p>Not wanting to think I just slapped it into <code>sympy</code>, which I haven’t really used up to that point:</p>
  
  <pre><code>&gt;&gt;&gt; from sympy import *
  &gt;&gt;&gt; x, y, z = symbols('x y z', integer=True, positive=True)
  &gt;&gt;&gt; solve(x*411 + y*295 + z*161 - 3200, [x, y, z])
  [{x: -295*y/411 - 161*z/411 + 3200/411}]
  </code></pre>
  
  <p>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) <code>[{x: 4, y: 2, z:6}]</code> or <code>[(4, 2, 6)]</code></p>
  
  <p>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.</p>