The Linear Diophantine Equation in Nonnegative Variables

Mathematica code is used to solve the linear diophantine equation for nonnegative integer values. Skiena's Mathematica implementation of  Dijkstra's shortest path algorithm is used for the problem solutions.

by  Harold Greenberg

Diophantus was a Greek mathematician who worked in Alexandria during the middle of the third century AD. The term "diophantine equation" is used when investigating solubility of algebraic equations in a particular number system, most often the set of integers.

We solve the diophantine equation 
					(∑_(j = 1)^n)(Subscript[a, j]) (Subscript[x, j])=L								(1)
for nonnegative integer values (Subscript[x, j]); the (Subscript[a, j]) and L are given positive integers. We consider the case where 1 < (Subscript[a, 1])<(Subscript[a, 2]) ⋯<(Subscript[a, n]) and the (Subscript[a, j]) are relatively prime. 
	If the (Subscript[x, j]) variables are unrestricted in sign, the solution to (1) is fairly easy (See Schrijver [1986] and the references cited therein.) The nonnegative integer condition on the (Subscript[x, j]) variables, however, prevents easy solution. Moreover, not all values of L lead to a feasible solution to (1). The integer (Subscript[x, j]) ≥ 0 condition also makes it difficult to determine whether or not a given value of L  allows (1) to have a solution. Schur has a theorem which states that there is an integer \(TraditionalForm\`\(L\^*\)\) such that for every \(TraditionalForm\`L > \(L\^*\)\) a solution exists. A problem of Frobenius is to find \(TraditionalForm\`\(L\^*\)\) as the largest value of L for which no solution exists. The history of the solution of problems like (1) for positive integers goes back to ancient times. See the references for an historical description and a list of works on diophantine problems and the Frobenius problem.
	We convert (1) to a dynamic programming recursive format.  We then obtain a network that is equivalent to the recursion and find the shortest path from a source node to all nodes in the network. The solution to (1) is obtained from the shortest path from the source node to a node that is characterized by the value L. The Frobenius value \(TraditionalForm\`\(L\^*\)\) is obtained from the longest of all the shortest paths to the nodes in the network. We adapt Skiena's Mathematica   implementation of Dijkstra's shortest path algorithm for our problem solution.

A Recursive Form for the Solution of (1)

We assert the following propositions:
1. There is an  \(TraditionalForm\`\(L\^*\)\) such that for all \(TraditionalForm\`L > \(L\^*\)\) a solution to (1) exists, while for \(TraditionalForm\`L = \(L\^*\)\) no solution exists.
2. An integer (Subscript[x, j]) ≥ 0 solution exists to the congruence
				(∑_(j = 2)^n)(Subscript[a, j]) (Subscript[x, j])≡y (mod (Subscript[a, 1]))								(2)
for any integer value y.
3. Consider y to be an integer variable and define the function F(y) by 
	F(y) = 
        Min [(∑_(j = 1)^n) (Subscript[a, j]) 
                (Subscript[x, j])  |  (∑_(j = 2)^n) (Subscript[a, j]) (Subscript[x, j]) ≡ 
              y  (mod (Subscript[a, 1])),  integer (Subscript[x, j]) ≥  0, 
          j = 1, 2, …, n].
The minimum exists for each y  from result 2 above; thus,
	F(y) = 
        Min [(∑_(j = 2)^n) (Subscript[a, j]) 
                (Subscript[x, j])  |  (∑_(j = 2)^n) (Subscript[a, j]) (Subscript[x, j]) ≡ 
              y  (mod (Subscript[a, 1])),  integer (Subscript[x, j]) ≥  0, 
          j = 2, …, n]				(3)
since optimal (Subscript[x, 1]) = 0.  Moreover, F(0) = 0. 
4. F(y)= F(y') when y  and y' are in the same residue class modulo (Subscript[a, 1]), i.e., y ≡ y'  (mod (Subscript[a, 1])). Hence, the calculation  for F(y) need be made only for the values y = 1,  2, …, (Subscript[a, 1]) -  1.
5. L and F(L) are in the same residue class modulo (Subscript[a, 1]); F(L)≡L (mod (Subscript[a, 1])). Thus,  given any integer L ≥ 0, a necessary and sufficient condition for a solution to (1) to exist is that L  ≥  F(L).
6. The Frobenius value is \(TraditionalForm\`\(L\^*\) = \(\(Max\)\(\ \)\([\)\)\)F(y)  |  y = 1,  
          2, …, (Subscript[a, 1]) -  1 ] - (Subscript[a, 1]).
7. The function F(y) may be written in the recursive form
				F\(TraditionalForm\`\((y)\) = \(\(Min\)\(\ \)\([\)\(\(a\_k\)\(+\)\)\)\) F(y - (Subscript[a, k]))  |  k = 2,  
      3, …, n]					(4)
where each y - (Subscript[a, k]) argument of F(y - (Subscript[a, k])) is taken modulo (Subscript[a, 1]) and  y = 1,  2, …,  
      (Subscript[a, 1]) -  1.
	The solution to the diophantine equation (1) is found below by first systematically finding F(y) for the values  y = 1,  2, …,  
      (Subscript[a, 1]) -  1 for (4). For any y, from the solution to (4), we will also find the (Subscript[x, j])values for j ≥ 2 for (3) . Hence, given L from (1), any integer solution to (1) must satisfy (2) with y ≡ L (mod (Subscript[a, 1])). With F(L), the corresponding values (Subscript[x, j]) produce the minimum possible value of (∑_(j = 2)^n) (Subscript[a, j]) (Subscript[x, j]). We are then left with (Subscript[a, 1]) (Subscript[x, 1])=L-F(L). If L ≥ F(L), then (Subscript[x, 1])=(L-F(L)/(Subscript[a, 1])to complete the solution. If L < F(L), then no feasible solution exists to (1).

The Network Analogy

The recursion for F(y) from (4) defines a network with nodes numbered 0, 1,  2, …, (Subscript[a, 1]) - 1. Each node x is incident to a directed arc of length (Subscript[a, k]), for k ≥ 2, that leads to node y = x + (Subscript[a, k]) (mod (Subscript[a, 1])), there being n - 1 arcs from node x. The required value F(y) is the shortest path distance to node y starting from node 0. The path follows a sequence of nodes and arcs from node 0 to node y.
Example. Consider  7  (Subscript[x, 1]) + 11  (Subscript[x, 2]) + 13  (Subscript[x, 3]) = 387. F(y) for (4) becomes F(y) = Min [11+  F(y - 11), 13+ F(y - 13)]. Hence, for y between 1 and 6, we obtain successively, 
F(1) = Min [11+ F(4), 13+ F(2)], F(2) = Min [11+ F(5), 13+ F(3)], F(3) = Min [11+ F(6), 13+ F(4)], F(4) = Min [11+ F(0), 13+ F(5)], F(5) = Min [11+ F(1), 13+ F(6)], F(6) = Min [11+ F(2), 13+ F(0)]. 
Node 0 is the start node and has a directed arc leading to node 4 of length 11 and a directed arc leading to node 6 of length 13. Similarly, node 4 has two entering arcs. One is the arc from node 0 with distance 11 and the second is the arc from node 5 with distance 13. In the same way, the entire network is built up from the recursion for F(y). Taking F(0) = 0, all other F(y) values are found using a shortest path algorithm. To solve the equation, we need to know F(y) for the node given by y ≡ 387  (mod 7), which results in y = 2 and represents node 2. Once the value of F(2) is determined, the equation has a solution if 387 ≥ F(2). When the equation has a solution, (Subscript[x, 2]) and (Subscript[x, 3]) are determined from the shortest path from node 0 to node 2. In addition,  (Subscript[x, 1]) = (387- F(2))/7.

The Algorithm

We adapt a Mathematica implementation [Skiena 1990] of Dijkstra's shortest path algorithm. Dijkstra's algorithm is also part of the Mathematica Standard Add-on Package DiscreteMath`Combinatorica`. For our example problem, we input the constants of \(TraditionalForm\`∑\+\(j = 1\)\%\(\(n\)\(\ \)\)	\(a\_j\) 
            x\_j = \ L\) as

a = {7, 11, 13}; L = 387;

We set up the problem for the format given by (4).

\(SetUp[a_] := \((b = a[[1]]; 
        c = Rest[a]; 
													e = 
          Table[Mod[c + x, b], {x, 0, b - 1}] + 
            1; 
														p = 
          Table[Infinity, {b}, {b}]; 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 		Do[
          p[[i, i]] = 0; 
          p[[i, 
              e[[
                i]]]] = c, {i, 
            b}])\)\)

The function Mod[c + x , b] gives a list of nodes that are adjacent to node x and are on the directed path from node x; hence, e produces the edge adjacency matrix for the network analogy. Whereas, in the example, node 0 is the start node and has directed arcs to nodes 4 and  6,  for ease of computation in Dijkstra's algorithm, we have made node 1 the start node which now leads to nodes 5 and 7;  this change is seen in e[[1]]={5, 7}. Hence, e[[k]] is a list of arcs leading from node k. The p matrix gives the directed distances between nodes; p[[i, j]] is the directed distance from node i to node j. The results for e and p in the sample problem are

SetUp[a];
e
p // MatrixForm

The setup for Dijkstra's algorithm is now complete. The algorithm follows as by Skiena where b equals the number of nodes and node 1 is the start node.

\(Skiena[e_, p_, b_] := 
      Module[{x}, parent = \(untraversed = Range[b]\); 

        dist = p[[1]]; 
        dist[[1]] = 0; 

        Scan[\((parent[[#]] = 1)\) &, 
          e[[1]]]; 

        While[untraversed != {}, 
	x = First[untraversed]; 
	
          Scan[\((If[
                  dist[[#]] < 
                    dist[[x]], 
                  x = #])\) &, untraversed]; 
	
          untraversed = Complement[untraversed, {x}]; 
	
          Scan[\((If[
                  dist[[#]] > 
                    dist[[x]] + 
                      p[[x, #]], 
                  dist[[#]] = 
                    dist[[x]] + 
                      p[[x, #]]; 
                  parent[[#]] = x])\) &, 
            e[[x]]]
]]\)

The results are lists of the parent nodes and the shortest distances along the shortest path to each node. Hence,

Skiena[e, p, b];
{parent, dist}

The parent list gives the predecessor node to a given node in the shortest path to the node. Hence, for example, node 3 has node 4 as predecessor node in the shortest path to node 3. Continuing to obtain the complete path to node 3: node 4 has node 5 as predecessor: node 5 has node 1 as predecessor node. Thus, the shortest path is given by 1->5->4->3. 
	The dist list gives the shortest distance to a given node. Hence, 37 is the shortest path distance to node 3. In terms of the original problem, the F(y) solutions to (4) are given by the dist list. Thus, F(0)=0, F(1)=22,  F(2)=37,F(3)=24,F(4)=11,F(5)=26,F(6)=13. Moreover, in terms of problem (1) and the definition of F(y) in (3) the dist list enables us to obtain all values of L for which no solution exists; e.g.,  for F(1)=22 there is no solution for L=1, 8, and 15, numbers smaller than 22 and congruent to 1 modulo 7.
	To determine the solution to a given diophantine equation after using SetUp and Skiena:

\(DiophSolution[L_, parent_, dist_, b_] := 
      Module[{y = Mod[L, b] + 1, x, path, r, s}, 
		If[
          L ≥ dist[[y]], 
          node[t_] := parent[[t]]; 
          path = Drop[FixedPointList[node, y], \(-1\)]; 
			\ \ r = 
            dist[[path]]; 
          s = Drop[r, \(-1\)] - 
              Drop[r, 1]; 
\ \ \ \ \ \ \ \ \ \ x = {\((L - 
                      dist[[y]])\)/b, 
                Map[Count[s, #] &, c]} // Flatten, {}]]\)

Note that we have changed the path calculation over that used by Skiena. Instead, we have used the function FixedPointList. The path obtained, however, is the reverse of Skiena's. The rest of the calculation finds the c values along the path and counts up how many times each individual c shows up to give the (Subscript[x, j])values. With L =387, b =7, and c ={11, 13}, we obtain

\(\(\(DiophSolution[L, parent, dist, b]\)\(			\)\)\)

Hence, (Subscript[x, 1]) = 50,  (Subscript[x, 2]) = 1, 
      (Subscript[x, 3]) = 2  solves  7  (Subscript[x, 1]) + 11  (Subscript[x, 2]) + 13  (Subscript[x, 4]) = 387. 
	The Frobenius value is then given by

Frobenius[dist_, b_] := Max[dist] - b
Frobenius[dist, b]

All values of L so that (1) has no solution are obtained from

f[x_] := Select[x - b, # > 0  &]
Sort[Flatten[FixedPointList[f,  Rest[dist] - b]]]

REFERENCES

GREENBERG, H. An Algorithm for a Linear Diophantine Equation and a Problem of Frobenius. Numer. Math, Vol. 34 (1980), pp. 349-352.
GREENBERG, H. Solution to a Linear Diophantine Equation for Nonnegative Integers. J. of Algorithms, Vol. 9 (1988), pp. 343-353.
JOHNSON, S. M. A Linear Diophantine Problem. Canad. J. Math., Vol. 12 (1960), pp. 390-8.
KNUTH, D. E. The Art of Computer Programming, Vol. 2, Seminumerical Algorithms. Addison-Wesley, Reading, MA (1981).
SCHRIJVER, A. The Theory of Linear and Integer Programming. John Wiley & Sons, New York (1986).
SKIENA, S. Implementing Discrete Mathematics: Combinatorics and Graph Theory with Mathematica, Addison-Wesley, Reading, MA  (1990), pp. 225-6.

ABOUT THE AUTHOR

Harold Greenberg uses Mathematica in his work on networks, cryptography, and mathematical programming. He is the author of the book  Integer Programming. In his spare time he writes humorous short stories. 


Harold Greenberg
Dept. of Statistics and Computer Information Systems
Baruch College/CUNY
17 Lexington Avenue
New York, NY 10010
dlobb@cunyvm.cuny.edu