Solving linear system over integers with numpy

Question

I'm trying to solve an overdetermined linear system of equations with numpy. Currently, I'm doing something like this (as a simple example):

a = np.array([[1,0], [0,1], [-1,1]])
b = np.array([1,1,0])

print np.linalg.lstsq(a,b)[0]

[ 1.  1.]

This works, but uses floats. Is there any way to solve the system over integers only? I've tried something along the lines of

print map(int, np.linalg.lstsq(a,b)[0])

[0, 1]

in order to convert the solution to an array of ints, expecting [1, 1], but clearly I'm missing something. Could anyone point me in the right direction?

Answer 1

I needed to do this and ended up porting a PHP program written by Keith Matthews, which you can find on http://www.numbertheory.org/php/php.html, into Python. I initially used Numpy arrays but ran into problems with integer overflows so switched to Sympy matrices, which use arbitrary precision numerical representations.

The code is released on GitHub at https://github.com/tclose/Diophantine under the MIT licence, so feel free to use it and let me know if you run into any problems (sorry it isn't better documented). The master branch uses Sympy but you can access the original Numpy implementation in the 'numpy' branch, which seems to work okay for reasonably sparse systems.

If you do end up using it for a scientific publication please cite Keith's papers (and maybe add a link to the GitHub repo).