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.