I have a large, under-determined system (60 equations and 116 unknowns) of linear Diophantine equations. I am aware of the algorithms typically used to solve these systems, which is not my question.

I have solved the system in Mupad and Mathematica for a parametric solution set. I have confirmed that the system has a bounded, finite set of integer solutions, and I am only interested in non-negative solutions.

I am looking for a way to:

derive (or estimate) the size of the reduced (i.e., non-negative) solution set, it seems that this can possibly be done using generating functions, although I can't find a good reference for performing this calculation in this instance; and

fully enumerate all non-negative solutions.

Given the size of the system, Mathematica seems to give up no matter how much memory I allocate to enumerating all of the solutions using `Reduce[]`

.

I am aware of a paper by Papp & Vizvari (2006) J. Mathematical Chemistry that solved such a system using several different algorithms and then enumerated all solutions following the method laid out in Land and Doig, but their paper does not give a lot of detail in the approach.

I would appreciate any direction to texts or resources that describe the estimation of the solution set and an efficient enumeration approach if they exist.

Thank you for your help.