A linear Diophantine equation solver

The Maude system includes a built-in linear Diophantine equation solver. The interface to the solver is defined in the file linear.maude which contains the functional module DIOPHANTINE. The current solver finds non-negative solutions of a system $S$ of $n$ simultaneous linear equations in $m$ variables having the form $M v = c$, where $M$ is an $n \times m$ integer coefficient matrix, $v$ is a column vector of $m$ variables and $c$ is a column vector of $n$ integer constants.

Both matrices and vectors are represented as sparse arrays with their dimensions implicit and their indices starting from 0. For this we make heavy use of the parameterized module ARRAY, described in Section 7.13.2.

First, a data type of pairs of natural numbers to be used as indices for matrices is created.

  fmod INDEX-PAIR is
    pr NAT .
    sort IndexPair .
    op _,_ : Nat Nat -> IndexPair [ctor] .
  endfm

Then, we instantiate (and rename as desired) the parameterized module ARRAY to obtain matrices of integers. Notice that Int0 is the view from DEFAULT to INT given in Section 7.11.2

  view IndexPair from TRIV to INDEX-PAIR is
    sort Elt to IndexPair .
  endv

  fmod MATRIX{X :: DEFAULT} is
    pr (ARRAY * (sort Entry{X,Y} to Entry{Y},
                 sort Array{X,Y} to Matrix{Y}))
       {IndexPair, X} .
  endfm

  fmod INT-MATRIX is
    pr MATRIX{Int0} * (sort Entry{Int0} to IntMatrixEntry,
                       sort Matrix{Int0} to IntMatrix,
                       op empty to zeroMatrix) .
  endfm

For example, the matrices

$$\left(\begin{array}{rr} 1 & 2 \\ 0 & -1 \end{array}\right) ... ...rrr} 1 & 2 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 0 \end{array}\right)$$

are both represented by the same term

  (0,0) |-> 1 ; (0,1) |-> 2 ; (1,1) |-> -1

Vectors are represented in a similar way as sparse arrays with natural numbers as indices. We use here the view Int0 already mentioned above and also the view Nat from TRIV to NAT given in Section 7.11.1. The view IntVector defined below will be used to construct sets of vectors later on.

  fmod VECTOR{X :: DEFAULT} is
    pr (ARRAY * (sort Entry{X,Y} to Entry{Y},
                 sort Array{X,Y} to Vector{Y}))
       {Nat, X} .
  endfm

  fmod INT-VECTOR is
    pr VECTOR{Int0} * (sort Entry{Int0} to IntVectorEntry,
                       sort Vector{Int0} to IntVector,
                       op empty to zeroVector) .
  endfm

  view IntVector from TRIV to INT-VECTOR is
    sort Elt to IntVector .
  endv

No distinction is made between row and column vectors, so, for example, both the row vector $\left( \begin{array}{cccc} -2 & 0 & 0 & 3 \end{array}\right)$ and its transpose $\left( \begin{array}{cccc} -2 & 0 & 0 & 3 \end{array}\right)^t$ are represented by the same term

  0 |-> -2 ; 3 |-> 3

The constants zeroMatrix and zeroVector denote the all zero matrix and vector, respectively.

The main module DIOPHANTINE begins defining pairs of sets of integer vectors, as follows:

  fmod DIOPHANTINE is
    pr STRING .
    pr INT-MATRIX .
    pr SET{IntVector} 
         * (sort NeSet{IntVector} to NeIntVectorSet,
            sort Set{IntVector} to IntVectorSet,
            op _,_ : Set{IntVector} Set{IntVector} -> Set{IntVector} 
                to (_,_) [prec 121 format (d d ni d)]) .

    sort IntVectorSetPair .
    op [_|_] : IntVectorSet IntVectorSet -> IntVectorSetPair
               [format (d n++i n ni n-- d)] .

Then, the solver is invoked with the built-in operator

    op natSystemSolve : IntMatrix IntVector String -> IntVectorSetPair
       [special (...)] .

which takes as arguments the coefficient matrix, the constant vector, and a string naming the algorithm to be used (see below), and returns the complete set of solutions encoded as a pair of sets of vectors [ $A$ | $B$ ]. The non-negative solutions of the linear Diophantine system correspond exactly to those vectors that can be formed as the sum of a vector from $A$ and a non-negative linear combination of vectors from $B$.

In particular, if the system $S$ is homogeneous (i.e., $c =$ zeroVector) then $A$ contains just the constant zeroVector and $B$ is the Diophantine basis of $S$ (which will be empty if $S$ only admits the trivial solution). A homogeneous system either has just the trivial solution or infinitely many solutions.

If $S$ is inhomogeneous (i.e., $c \neq$ zeroVector) then, if $S$ has no solution, both $A$ and $B$ will be empty; otherwise, $B$ will consist of the Diophantine basis of $S'$, the system formed by setting $c =$ zeroVector, while $A$ contains all solutions of $S$ that are not strictly larger than any element of $B$. An inhomogeneous system may have no solution (in this case $A$ and $B$ are both empty), a finite number of solutions (in this case $A$ is non-empty and $B$ is empty), or infinitely many solutions (in this case $A$ and $B$ are both non-empty).

In either case, the solution encoding [ $A$ | $B$ ] is unique.

Deciding whether a linear Diophantine system admits a non-negative, nontrivial solution is NP-complete (stated as known in [70]). Furthermore the size of the Diophantine basis of a homogeneous system can be very large. For example the equation: $x + y - k z = 0$, for constant $k > 0$, has a Diophantine basis (i.e., set of minimal, nontrivial solutions) of size $k + 1$.

There are currently two algorithms implemented.

The string "cd" specifies a version of the classical Contejean-Devie algorithm [22] with various improvements. The algorithm is based on incrementing a vector of counters, one for each variable, and so it can only solve systems where the answers involve fairly small numbers. It is fairly insensitive to the number of degrees of freedom in the problem. The improvements in this implementation take effect when an equation has zero or one unfrozen variables with nonzero coefficients and result in either forced assignments or early pruning of a branch of the search. It performs well on the following homogeneous system from [26],

$$\left(\begin{array}{rrrrrr} 1 & 2 & -1 & 0 & -2 & -1 \\ 0 &... ...ght) = \left(\begin{array}{c} 0 \\ 0 \\ 0 \end{array}\right)$$

which has a basis of size 13.

  Maude> red in DIOPHANTINE : 
           natSystemSolve(
             (0,0) |-> 1 ; (0,1) |-> 2 ; (0,2) |-> -1 ;
             (0,3) |-> 0 ; (0,4) |-> -2 ; (0,5) |-> -1 ;
             (1,0) |-> 0 ; (1,1) |-> -1 ; (1,2) |-> -2 ;
             (1,3) |-> 2 ; (1,4) |-> 0 ; (1,5) |-> 1 ;
             (2,0) |-> 2 ; (2,1) |-> 0 ; (2,2) |-> 1 ;
             (2,3) |-> -1 ; (2,4) |-> -2 ; (2,5) |-> 0,
             zeroVector,
             "cd") .
  rewrites: 1 in 10ms cpu (46ms real) (100 rews/sec)
  result IntVectorSetPair: 
    [
      zeroVector
    |
      0 |-> 1 ; 1 |-> 1 ; 4 |-> 1 ; 5 |-> 1,
      0 |-> 1 ; 1 |-> 4 ; 2 |-> 9 ; 3 |-> 11,
      0 |-> 10 ; 1 |-> 4 ; 3 |-> 2 ; 4 |-> 9,
      1 |-> 1 ; 2 |-> 1 ; 3 |-> 1 ; 5 |-> 1,
      1 |-> 8 ; 2 |-> 2 ; 4 |-> 1 ; 5 |-> 12,
      0 |-> 2 ; 1 |-> 4 ; 2 |-> 8 ; 3 |-> 10 ; 4 |-> 1,
      0 |-> 3 ; 1 |-> 4 ; 2 |-> 7 ; 3 |-> 9 ; 4 |-> 2,
      0 |-> 4 ; 1 |-> 4 ; 2 |-> 6 ; 3 |-> 8 ; 4 |-> 3,
      0 |-> 5 ; 1 |-> 4 ; 2 |-> 5 ; 3 |-> 7 ; 4 |-> 4,
      0 |-> 6 ; 1 |-> 4 ; 2 |-> 4 ; 3 |-> 6 ; 4 |-> 5,
      0 |-> 7 ; 1 |-> 4 ; 2 |-> 3 ; 3 |-> 5 ; 4 |-> 6,
      0 |-> 8 ; 1 |-> 4 ; 2 |-> 2 ; 3 |-> 4 ; 4 |-> 7,
      0 |-> 9 ; 1 |-> 4 ; 2 |-> 1 ; 3 |-> 3 ; 4 |-> 8
    ]

The string "gcd" specifies an original algorithm based on integer Gaussian elimination followed by a sequence of extended greatest common divisor (gcd) computations. It can ``home in'' quickly on solutions involving large numbers but it is very sensitive to the number of degrees of freedom and can easily degenerate into a brute force search. Furthermore, termination depends on the bound on the sum of minimal solutions established in [66], which can cause a huge amount of fruitless search after the last minimal solution has been found. It performs well on the ``sailors and monkey'' problem from [22]:

  red in DIOPHANTINE : 
    natSystemSolve(
      (0,0) |-> 1 ; (0,1) |-> -5 ; (1,1) |-> 4 ; (1,2) |-> -5 ;
      (2,2) |-> 4 ; (2,3) |-> -5 ; (3,3) |-> 4 ; (3,4) |-> -5 ;
      (4,4) |-> 4 ; (4,5) |-> -5 ; (5,5) |-> 4 ; (5,6) |-> -5,
      0 |-> 1 ; 1 |-> 1 ; 2 |-> 1 ; 3 |-> 1 ; 4 |-> 1 ; 5 |-> 1,
      "gcd") .
  result IntVectorSetPair: 
    [
      0 |-> 15621 ; 1 |-> 3124 ; 2 |-> 2499 ; 3 |-> 1999 ;
      4 |-> 1599 ; 5 |-> 1279 ; 6 |-> 1023
    |
      0 |-> 15625 ; 1 |-> 3125 ; 2 |-> 2500 ; 3 |-> 2000 ;
      4 |-> 1600 ; 5 |-> 1280 ; 6 |-> 1024
    ]

Finally, the string "" can be passed as third argument of natSystemSolve, thus allowing the system to choose which algorithm to use. For convenience, the operator

    op natSystemSolve : IntMatrix IntVector -> IntVectorSetPair .

is equationally defined to invoke the built-in operator with ""

    eq natSystemSolve(M:IntMatrix, V:IntVector) 
      = natSystemSolve(M:IntMatrix, V:IntVector, "") .
  endfm