Using SymPy's Diophantine module, it is possible to obtain a parameterised version for any linear diophantine equation. Solutions can then be found for the equation by iterating over each parameterised term, as demonstrated by the solution to this particular question.
However, the solution to that question was done by manually finding the ranges for each parameter to iterate over in a two-variable problem to obtain positive integral solutions. My program is supposed to find natural number solutions for a linear Diophantine equation of a variable number of terms (expected 3 to about 12 at most), so it must be able to determine the bounds to iterate over for each parameter by itself. All terms are additive barring the final integer target, for example:
12*A + B + 14*C + 16*D + 31*E - 507
(t_0, t_0 + t_1, t_0 + t_1 + t_2, -23*t_0 + t_1 + 3*t_2 + 31*t_3 + 1014, 11*t_0 - t_1 - 2*t_2 - 16*t_3 - 507)
Defining the number of nested loops is fairly simple (iterate over the solution set and use SymPy's free_symbols function to extract them all), and itertools has the product function for setting up that kind of nested loop, as shown here. What I want to find is the range for each t parameter, such that I find all natural number solutions and dont get stuck in an infinite loop. How might I go about this?