MATLAB functions for Diophantine equations

John Burkardt 21 May 2020

We suppose we are given a linear Diophantine equation in n variables:

∑_{i=1}^{n} a_{i}x_{i} = b

where the coefficients a_{i} and the right hand side b are assumed to be strictly positive integers. It is desired to determine all solution vectors x which are integers, and either nonnegative, or strictly positive.

Under the assumptions on a and b, there can no more than finitely many nonnegative or strictly positive integer solutions. We may use k to represent the number of solutions of a given problem. Of course, k may be zero.

This text describes a pair of MATLAB functions, with the signatures

  x = diophantine_nd_nonnegative(a,b)
  x = diophantine_nd_positive(a,b)

which accept the coefficient n-vector a and right hand side b of a Diophantine equation, and return the set of k solutions in an n × k array x.

We will regard the nonnegative solution procedure as the more fundamental of the two. If we have such a procedure, then we can easily construct a corresponding procedure for all strictly positive solutions. Of course, a brute force approach would generate all nonnegative solutions and then discard those that are not strictly positive. A more satisfactory approach defines an auxilliary problem that can be treated by the nonnegative procedure, whose result can then be transformed into a solution of the strictly positive problem, and this is the method we employ.

Thus, assume for a moment that we have a nonnegative solver. Say we are given the information for a Diophantine equation, but we wish to find only the strictly positive solutions. In that case, from the data a and b, we can define a new right hand side β, and use the nonnegative solution procedure to find solution vectors χ, from which the solution set x of the original problem can be recovered.

The following algorithm suggests how this is done:

Algorithm 1 Solver for positive solutions of Diophantine equation

procedure diophantine nd positive(a, b) n ← length of a β ← b − ∑_{i=1}^{n} a_{i} if β < 0 then x ← ∅

else χ ← diophantine nd nonnegative(a, β). for 1 ≤ j ≤ n do x_{j} ← χ_{j} + 1

end for end if return x

end procedure Now we can concentrate on the algorithm diophantine nd nonnegative() which seeks nonnegative solutions to a Diophantine problem. Given the data a and b, we plan to construct each possible solution vector x one component at a time. The step counter j starts at 1. We start by evaluating the “surplus”, r = b − ∑ i = 1^{j−1}a_{i}x_{j}. On step j, we assume we have already set components 1 through j. If j < n, then we can increment j and choose x_{j} to be any value between 0 and ⌊\frac{r}{a_{j}} ⌋, taking the greedy choice x_{j} = ⌊\frac{r}{a_{j}} ⌋, and we are ready for the next step.

However, if j = n, we have constructed a complete proposed solution. If the solution satisfies the equation, we add it to our list. In any case, we then search our current solution x for the largest previous index k for which x_{k} was nonzero. If we find such an index, we decrement x_{k} by 1, set j = k, and restart the algorithm from that index to consider the next possible solution. But if no such index k exists, then we exhausted our options and the search terminates.

An outline of this algorithm follows:

Algorithm 2 Solver for nonnegative solutions of Diophantine equation

procedure diophantine nd nonnegative(a, b) n ← length of a x ← ∅ j ← 0 r ← 0 while true do ▷ Construct a vector Y that is a possible solution.

r ← b − ∑_{i=1}^{j} a_{i} y_{i}

if j < n then ▷ We have a partial vector Y. Get next component. j ← j + 1 y(j) ← ⌊(r/a(j)⌋)

else ▷ We have a full vector Y. if r == 0 then ▷ Is Y a solution?

x ← x ∪ y end if while 0 < j do ▷ Find last nonzero Y entry, decrease by 1.

if 0 < y(j) then y(j) ← y(j) − 1 break

end if j ← j − 1

end while if j == 0 then return x

end if end if end while

end procedure

Should efficiency be vital, there are several ways to quicken the procedure. At an early step j < n, if the residual r is zero, then the remaining entries of y can be set to zero immediately and a solution accepted. If, at any step j, the residual is not divisible by the greatest common divisor of the coeffients indexed from j +1 through n, then the partial proposed solution is untenable; the last nonzero entry in y should be decremented and the search restarted from that index. For the moment, however, these possible efficiencies have not been implemented.

As a small demonstration of this procedure, consider the calculation of the nonnegative solutions of

12 x_{1} + 9 x_{1} + 7 x_{3} = 30

The following diagram suggests the progress of the algorithm. The first line suggests that algorithm begins by setting y_{1} to its maximum value of 2. It then produces the partial vector (2,0) and the full vector (2,0,0), which is rejected because it does not satisfy the equation. On the second line, we see that the algorithm steps back to the last nonzero component of y, namely y_{1}, reduces it by 1, and extends the vector to (1,2) and then (1,2,0). This vector satisfies the equation and so is accepted as the first solution. In the third line, we see that the algorithm steps back to the last nonzero component of y, namely y_{2}, decrements it to 1, and generates the vector (1,1,1), which is rejected. Subsequent lines show the orderly generation of tentative vectors, and the discovery of a second solution.

 2; 2 0; 2 0 0
 1; 1 2; 1 2 0  Solution #1
    1 1; 1 1 1
         1 1 0
    1 0; 1 0 2
         1 0 1
         1 0 0
 0; 0 3; 0 3 0
    0 2; 0 2 1
         0 2 0
    0 1; 0 1 3  Solution #2
         0 1 2
         0 1 1
         0 1 0
    0 0; 0 0 4
         0 0 3
         0 0 2
         0 0 1
         0 0 0

As a result of the way that the possibilities are examined, the accepted solutions are automatically accumulated in a sorted order. If you regard each solution as the digits of a number, they show up in descending numeric order. For example, the k = 10 strictly positive solutions of 2 x_{1} + 3 x_{2} + 5 x_{3} + 6 x_{4} + 7 x_{5} = 35 are returned in a 5 × 10 array x. To illustrate the ordering, we print the transpose of x, so that each solution appears as a row, and can be read as a 5 digit number:

  7,1,1,1,1
  4,3,1,1,1
  4,1,1,2,1
  3,2,2,1,1
  2,2,1,1,2
  2,1,3,1,1
  1,5,1,1,1
  1,3,1,2,1
  1,1,2,1,2
  1,1,1,3,1