Consider the following example Diophantine equation

-118w + 989x + 918y  -512z == 0

You can solve this in sympy with:

import sympy
from sympy.solvers.diophantine import diophantine 
w, x, y, z = sympy.symbols('w x y z')
a, b, c, d = -118, 989, 918, -512
diof = list(diophantine(a*w + b*x +c*y +d*z ))[0]
print(diof)

which gives:

(t_0, 118*t_0 + 2*t_1, 1690468*t_0 + 28937*t_1 + 256*t_2, 3031184*t_0 + 51887*t_1 + 459*t_2)

Explanation: This solution allows you to pick three integers t_0, t_1, t_2 freely and from there you can compute the values of w, x, y, z.

This is fine except the coefficients are always much larger than they need to be. If you solve the same system in Mathematica you get:

x = 2 c_2, y = 9 c_2 + 256 c_3 + 81 w, z = 20 c_2 + 459 c_3 + 145 w, c_2 element Z and c_3 element Z

or:

Rewritten in the same style as the sympy solution this would be:

(t_0, 2*t_1, 9*t_1 + 256*t_2 + 81*t_0, 20*t_1 + 459*t_2 + 145*t_0)

Is it possible to get sympy to give solutions to linear Diophantine equations with similarly small coefficients?

Typos fixed in the sympy code.