Algorithm for solving systems of linear Diophantine inequalities
Mathematica implements an algorithm: see the manual here:
(Added in response to a comment query.)
The paper
Hochbaum, Dorit S., and Anu Pathria. "Can a System of Linear Diophantine Equations be Solved in Strongly Polynomial Time?." (1994). (PDF download link)
says that "no strongly polynomial algorithm exists for the problem of finding the set of solutions to a system of linear diophantine equations in the complexity model allowing the operations $\{ +, -, \times, /, \bmod, < \}$." But if one assumes the gcd can be found quickly, then the system can be solved in strongly polynomial time.
GAP provides a function NullspaceIntMat
which solves systems
of linear diophantine equations. The documentation says:
25.1-2 SolutionIntMat
* SolutionIntMat( mat, vec ) ───────────────────────────────────── operation
If mat is a matrix with integral entries and vec a vector with integral
entries, this function returns a vector x with integer entries that is a
solution of the equation x * mat = vec. It returns fail if no such vector
exists.
──────────────────────────────── Example ─────────────────────────────────
gap> mat:=[[1,2,7],[4,5,6],[7,8,9],[10,11,19],[5,7,12]];;
gap> SolutionMat(mat,[95,115,182]);
[ 47/4, -17/2, 67/4, 0, 0 ]
gap> SolutionIntMat(mat,[95,115,182]);
[ 2285, -5854, 4888, -1299, 0 ]
────────────────────────────────────────────────────────────────────────────
The source code can be found in the file lib/matint.gi included in the GAP distribution.
There is also a function SolutionNullspaceIntMat
which additionally
computes a basis of the integral nullspace of the given matrix.