Solve linear equations over the integers

Question

All variables in this question are integer valued.

Input

4 integers w, x, y, z. They can be positive or negative and will be less than 1048576 in absolute value.

Output

The general solution to the equation.

\$ aw+bx+cy+dz = 0 \$.

The variables \$a, b, c, d\$ must all be integer values.

Output format

Your output should consist of three tuples each with four parts, one for each of the values a, b, c, d. Let me explain by example:

Input: -118, 989, 918, -512

Mathematical answer:
        b = 2 n_0 
        c = 9 n_0 + 256 n_1 + 81 a
        d = 20 n_0 + 459 n_1 + 145 a

Explanation: n_0 and n_1 are integers that you can set to anything you like. The solution says: a can also be set to any integer value, b must be twice whatever you set n_0 to. This means that a can be set to any integer, c can now be calculated in terms of three variables we have already set and so can d.

The format of your output should be 3 tuples (#,#,#,#), (#,#,#,#), (#,#,#,#). We can assume three free integer variables n0, n1 and n2 and so (a,b,c,d) = (#,#,#,#)n0 + (#,#,#,#)n1 + (#,#,#,#)n2. In the example above the output would therefore be:

Output: (0, 2, 9, 20), (0, 0, 256, 459), (1, 0, 81, 145)

Examples

Example one:

 Input: -6, 3, 7, 8

 Mathematical answer: 
          c = 2a + 3b + 8n
          d = -a - 3b - 7n 
          n is any integer

 Output: (1, 0, 2, -1), (0, 1, 3, -3), (0, 0, 8, -7)

Example two:

Input: -116, 60, 897, 578

Mathematical answer: 
        c = 578 n + 158 a + 576 b
        d = -897 n - 245 a - 894 b 
        n is any integer

Output: (1, 0, 158, -245), (0, 1, 576, -894), (0, 0, 578, -897)

Example three:

Input: 159, -736, -845, -96

Output: (1, 0, 27, -236), (0, 1, 64, -571), (0, 0, 96, -845)

Discussion

To understand this challenge further it is worth looking at this possible general solution which does not work [(z, 0, 0, -w), (0, z, 0, -x), (0, 0, z, -y)]. The problem with this is that there are solutions to the problem instances above which are not the sum of any integer multiples of those tuples. For example: take input -6, 3, 7, 8 from Example 1. The proposed solution would give the tuples:

(8, 0, 0, 6), (0, 8, 0, -3), (0, 0, 8, -7)

Why doesn't this work?

There is a solution for this instance with a = 1, b = 1, c = 13, d = -11 because -6+3+7*13-11*8 = 0. However there are no integers n_0, n_1, n_2 to make n_0 * (8, 0, 0, 6) + n_1 * (0, 8, 0, -3) + n_2 * (0, 0, 8, -7) = (1, 1, 13, -11) .

Answer 1

Python 2, 182 bytes

def E(a,b):
 if b:x,y,d=E(b,a%b);return y,x-a/b*y,d
 return 1,0,a
def f(a,b,c,d):x,y,g=E(a,b);z,w,h=E(c,d);j=E(g,h)[2];return(b/g,-a/g,0,0),(0,0,d/h,-c/h),(-x*h/j,-y*h/j,z*g/j,w*g/j)

Try it online!

Yes, it is possible to solve it without fancy built-ins. And this gives relatively good results for a non-LLL one.

How it works

I derived a closed-form solution in the following steps: (I flipped the roles of \$a,b,c,d\$ and \$w,x,y,z\$ because it looked much more natural to me. So \$a,b,c,d\$ are coefficients and \$w,x,y,z\$ are unknowns here.)

$$ ax+by+cz+dw = 0 $$

If we define \$g=\gcd(a,b)\$ and \$h=\gcd(c,d)\$, the following substitutions can be done without loss of generality: (\$u,v\$ are temporary integer variables)

$$ ax+by = gu\\ cz+dw = hv $$

which results in a system of two heterogeneous and one homogeneous 2-variable linear Diophantine equations:

$$ ax+by = gu \\ cz+dw = hv \\ gu + hv = 0 $$

Then we can apply the textbook procedure:

$$ u = \frac{h}{j}n_3, \; v = \frac{-g}{j}n_3 \quad \text{where} \; j = \gcd(g,h) $$

If we let \$(x_0,y_0)\$ be a "special" solution to \$ax+by=g\$ and similarly \$(z_0,w_0)\$, then the solutions are

$$ \begin{align} x &= \frac{b}{g}n_1 + \frac{-x_0h}{j}n_3 \\ y &= \frac{-a}{g}n_1 + \frac{-y_0h}{j}n_3 \\ z &= \frac{d}{h}n_2 + \frac{z_0g}{j}n_3 \\ w &= \frac{-c}{h}n_2 + \frac{w_0g}{j}n_3 \end{align} $$

The special solutions can be found via Extended Euclidean algorithm. Special mention goes to this deceptively short and yet correct implementation of EGCD, once added to sympy but removed later:

def extended_euclid(a, b):
    if b == 0:
        return (1, 0, a)

    x0, y0, d = extended_euclid(b, a%b)
    x, y = y0, x0 - (a//b) * y0

    return x, y, d

Answer 2

Wolfram Language (Mathematica), 39 bytes

Still a built-in.

HermiteDecomposition gives the Hermite normal form decomposition of an integer matrix.

HermiteDecomposition[List/@#][[1,2;;]]&

Try it online!

---

Wolfram Language (Mathematica), without build-in, 118 bytes

Based on Algorithm 2.4.10 in A Course in Computational Algebraic Number Theory by Henri Cohen.

The output is sometimes extremely large. You can reduce them with the build-in LatticeReduce.

(a=#[[1]];U=IdentityMatrix@4;Do[{d,u}=ExtendedGCD[a,b=#[[j]]];U[[{1,j}]]={u,{-b,a}/d}.U[[{1,j}]];a=d,{j,2,4}];Rest@U)&

Try it online!

Answer 3

SageMath, 92 bytes

def f(x,y,z,w):
 for v in'abcd':var(v,domain=ZZ)
 return solve([a*x+b*y+c*z+d*w==0],a,b,c,d)

Try it online!

Answer 4

Pari/GP, 21 bytes

Another boring built-in answer.

matkerint(x,{flag=0}): LLL-reduced Z-basis of the kernel of the matrix x with integral entries. flag is deprecated, and may be set to 0 or 1 for backward compatibility.

a->matkerint(Mat(a))~

Try it online!