<< Computer Algorithms
Use  Ctrl +  or    to enlarge or reduce text size.

How to Solve a System of Linear Equations
in the Integers


The problem is to solve the system of linear equations
A x  =  b
where
A  is a square matrix,
x  is an unknown column array,
b  is a column array,
and every number is an integer.  In other words, a Diophantine system of linear equations.

The program below solves this problem by finding a diagonal form of A (call it D) and an associated pair of matrices whose determinants are 1 (call them U and V) such that  U A V  =  D.

Then there is an integral solution when (and only when) the elements of
t  =  D-1 U b
are integers, and in that case the solution is
x  =  V t.
(Multiplying  y = U b  by  D-1  just divides y by the diagonal elements of  D  because  D  is diagonal.  The difficult part is computing D.)

No fractions or real numbers are used in the computation, all variables are declared long or quad integers.

NOTE:  In general D is not the Smith normal form of A, which requires that each diagonal element divide the next.

The program is written for PBCC and can be modified for Visual Basic, FORTRAN and any other BASIC-like compiler.

' Copyright 2015.

' Problem: Solve the system of linear equations
'   Ax = b
' for x where
'   A is a known non-singular  n x n  matrix
'   x is an unknown column vector of n elements
'   b is a known column vector of n elements
' and every number is an integer,
' or determine that there is no solution using integers.
'
' Method:  Find a precursor to the Smith normal form of A and associated matrices.

%loquacious = 1            '1 means print intermediate results, 0 means do not

' Given A(), b(), n.
' Return x() by reference, -1 if no integral solution, -2 if A() is singular.
Function SolveDiophantineSystem(A() As Long, b() As Long, ByVal n As Long, x() As Long) As Long
 Local i, j As Long        'indices
 Local k As Quad           'a matrix diagonal value
 Local s As Long           'stage of computation, the row and/or column we are working on
 Local p, q As Quad        'pair of matrix values
 Local f As Quad           'subtraction factor
 Local minval, minvalabs As Quad   'value of smallest non-zero element in current submatrix
 Local minrow, mincol As Long      'its location
 Local sum As Quad         'used in matrix multiplication
 Dim D(n,n) As Quad        'at first a copy of the input matrix A()
 Dim t(n) As Quad          'pre-solution vector
 Dim U(n,n) As Quad        'left auxiliary matrix, at first the identity matrix
 Dim V(n,n) As Quad        'right auxiliary matrix, at first the identity matrix

 #If %loquacious
 Local ii, jj As Long
 Dim Q(n,n) As Long
 #EndIf

 'Start with  U = I,  D = A,  V = I  at stage 0.
 'Execute row and column operations on D, using only integers, to transform it to diagonal form.
 '(D does not end up as the Smith Normal Form, which would require additional row and column
 'operations so that each diagonal value divides the next.)
 'Simultaneously do each row operation on U and each column operations on V.
 'Let s = the current stage, s = 0 to n.

 For i = 1 To n
  U(i,i) = 1               'at first  U = I  and  V = I
  V(i,i) = 1
  For j = 1 To n
   D(i,j) = A(i,j)         'work with a copy of A
  Next
 Next

 s = 0                     'stage

 Do
  Incr s                   'next stage

Pivot:
  'Look at the submatrix of D from s,s to n,n.  Find the smallest
  '(in magnitude) element, call it minval at minrow,mincol.
  If D(s,s) = 0 Then              'then swap in a row or column so it's not
   For j = s + 1 To n             'see if can swap columns to fix it and do same to V
    If D(s,j) Then
     For i = s To n
      Swap D(i,s), D(i,j)
     Next
     For i = 1 To n
      Swap V(i,s), V(i,j)
     Next
     Exit For
    End If
   Next
   If j <=n Then Exit If         'fixed it
   For i = s + 1 To n            'see if can swap rows to fix it and do same to U
    If D(i,s) Then
     For j = s To n
      Swap D(s,j), D(i,j)
     Next
     For j = 1 To n
      Swap U(s,j), U(i,j)
     Next
     Exit For
    End If
   Next
   If i > n Then                 'failed, singular matrix
    Function = -2
    Exit Function
   End If
  End If

  minval = D(s,s) :minvalabs = Abs(minval)
  minrow = s :mincol = s
  For i = s To n      'rows
  For j = s To n      'columns
   If D(i,j) Then
    If Abs(D(i,j)) < minvalabs Then
     minval = D(i,j) :minvalabs = Abs(minval)
     minrow = i :mincol = j
    End If
   End If
  Next
  Next
  If minval = 0 Then               'singular matrix
   Function = -2
   Exit Function
  End If

  'Swap rows and columns of D so this element is moved to s,s.
  'Do the same row operation on U, the same column operation on V.
  If minrow <> s Then              'swap rows
   For j = s To n
    Swap D(s,j), D(minrow,j)
   Next
   For j = 1 To n
    Swap U(s,j), U(minrow,j)
   Next
  End If
  If mincol <> s Then              'swap columns
   For i = s To n
    Swap D(i,s), D(i,mincol)
   Next
   For i = 1 To n
    Swap V(i,s), V(i,mincol)
   Next
  End If

  #If %loquacious
  Print "pivot D ="
  For ii = 1 To n
  For jj = 1 To n
   Print D(ii, jj);
  Next
  Print
  Next
  Print
  #EndIf

  'Kill elements in the s-row to the right of s using column operations.
  'As this is done do the same thing to the columns of V.
  For j = s + 1 To n                       'columns
   p = D(s,s) :q = D(s,j)                  'p <> 0, want q --> 0
   Do
    If     q = 0 Then
     f = 0  :Exit Do
    ElseIf p = 1 Then
     f = q  :Exit Do
    ElseIf p = -1 Then
     f = -q :Exit Do
    ElseIf p = q Then
     f = 1  :Exit Do
    ElseIf p = -q Then
     f = -1 :Exit Do
    End If
    If Abs(p) < Abs(q) Then
     f = q \ p
     q = q - p*f
     For i = s To n                         'columns s and j
      D(i,j) = D(i,j) - D(i,s) * f
     Next
     For i = 1 To n
      V(i,j) = V(i,j) - V(i,s) * f
     Next
    Else
     f = p \ q
     If q*f = p Then Decr f
     p = p - q*f
     For i = s To n
      D(i,s) = D(i,s) - D(i,j) * f
     Next
     For i = 1 To n
      V(i,s) = V(i,s) - V(i,j) * f
     Next
    End If
   Loop
   If f Then                                'q = q - p*f
    For i = s To n
     D(i,j) = D(i,j) - D(i,s) * f
    Next
    For i = 1 To n
     V(i,j) = V(i,j) - V(i,s) * f
    Next
   End If

   #If %loquacious
   Print "column operation D ="
   For ii = 1 To n
   For jj = 1 To n
    Print D(ii, jj);
   Next
   Print
   Next
   Print
   #EndIf

  Next

  'Kill elements in the s-column below s using row operations.
  For i = s + 1 To n                       'rows
   p = D(s,s) :q = D(i,s)
   Do
    If     q = 0 Then
     f = 0  :Exit Do
    ElseIf p = 1 Then
     f = q  :Exit Do
    ElseIf p = -1 Then
     f = -q :Exit Do
    ElseIf p = q Then
     f = 1  :Exit Do
    ElseIf p = -q Then
     f = -1 :Exit Do
    End If
    If Abs(p) < Abs(q) Then
     f = q \ p
     q = q - p*f
     For j = s To n                         'rows s and i
      D(i,j) = D(i,j) - D(s,j) * f
     Next
     For j = 1 To n
      U(i,j) = U(i,j) - U(s,j) * f
     Next
    Else
     f = p \ q
     If q*f = p Then Decr f
     p = p - q*f
     For j = s To n
      D(s,j) = D(s,j) - D(i,j) * f
     Next
     For j = 1 To n
      U(s,j) = U(s,j) - U(i,j) * f
     Next
    End If
   Loop
   If f Then                                'q = q - p*f
    For j = s To n
     D(i,j) = D(i,j) - D(s,j) * f
    Next
    For j = 1 To n
     U(i,j) = U(i,j) - U(s,j) * f
    Next
   End If

   #If %loquacious
   Print "row operation D ="
   For ii = 1 To n
   For jj = 1 To n
    Print D(ii, jj);
   Next
   Print
   Next
   Print
   #EndIf

  Next

  'check if lines at s,s are 0 except at s,s
  For j = s + 1 To n             'columns
   If D(s,j) <> 0 Then GoTo Pivot
  Next
  For i = s + 1 To n             'rows
   If D(i,s) <> 0 Then GoTo Pivot
  Next

 Loop Until s = n

 #If %loquacious
 Locate CursorY - 1
 Print "____________________________________" :Print
 Print " U ="
 For i = 1 To n
 For j = 1 To n
  Print U(i, j);
 Next
 Print
 Next
 Print

 Print " V ="
 For i = 1 To n
 For j = 1 To n
  Print V(i, j);
 Next
 Print
 Next
 Print

 Print " UAV ="
 For i = 1 To n                                 'Q = AV
 For j = 1 To n
  sum = 0
  For k = 1 To n
   sum = sum + A(i,k) * V(k, j)
  Next
  Q(i,j) = sum
 Next
 Next
 For i = 1 To n                                 'UQ
 For j = 1 To n
  sum = 0
  For k = 1 To n
   sum = sum + U(i,k) * Q(k, j)
  Next
  Print sum;
 Next
 Print
 Next
 Print

 #EndIf

 For i = 1 To n
  sum = 0
  For j = 1 To n
   sum = sum + U(i,j) * b(j)
  Next
  k = D(i,i)            'n = n0 And
  If (sum \ k) * k <> sum Then    'if D(i,i) doesn't divide sum = c(i), no integral solution
   Function = -1
   Exit Function
  End If
  t(i) = sum \ k
 Next

 #If %loquacious
 Print " t ="
 For i = 1 To n
  Print t(i)
 Next
 Print "____________________________________" :Print
 #EndIf

 For i = 1 To n
  sum = 0
  For j = 1 To n
   sum = sum + V(i,j) * t(j)
  Next
  x(i) = sum
 Next
End Function
'---------------------------------------------------

' | a11  a12  a13 ... a1n |   | x1 |       | b1 |
' | a21  a22  a23 ... a2n |   | x2 |       | b2 |
' | ..   ..   ..      ..  |   | .. |   =   | .. |
' | an1  an2  an3 ... ann |   | xn |       | bn |
'
Function PBMain
 Local i, j As Long        'row and column indices
 Local n As Long           'number of rows, columns of A()
 Local solvable As Long    'error code
 Local sum As Long

 'For example
 n = 3

 Dim A(1 To n, 1 To n) As Long    'given matrix
 Dim b(1 To n) As Long            'given vector
 Dim x(1 To n) As Long            'solution vector

 Dim xsave(n) As Long

 'For example
' A(1,1) = 2 :A(1,2) = 1  :A(1,3) = 2      :b(1) = 30
' A(2,1) = 3 :A(2,2) = -1 :A(2,3) = 0      :b(2) = 9
' A(3,1) = 4 :A(3,2) = 1  :A(3,3) = 5      :b(3) = 61
' A(1,1) = 5 :A(1,2) = 1  :A(1,3) = 2      :b(1) = 29
' A(2,1) = 3 :A(2,2) = 3  :A(2,3) = 2      :b(2) = 31
' A(3,1) = 2 :A(3,2) = 1  :A(3,3) = 5      :b(3) = 35
 A(1,1) = 5 :A(1,2) = 1  :A(1,3) = 2      :b(1) = 13
 A(2,1) = 3 :A(2,2) =  1 :A(2,3) = 2      :b(2) = 7
 A(3,1) = 2 :A(3,2) = 1  :A(3,3) = 5      :b(3) = 19
'---------------------------------------------------

 Print " Ax  =  b"
 Print
 For i = 1 To n
 For j = 1 To n
  Print A(i, j);
 Next
 Print ,"x";Format$(i),b(i)
 Next
 Print "____________________________________" :Print

 solvable = SolveDiophantineSystem ( A(), b(), n, x() )

 If     solvable = -2 Then
  Print "The matrix is singular."
 ElseIf solvable = -1 Then
  Print "There is no solution in integers."
 Else
  Print " x ="
  For i = 1 To n
   Print "    ";x(i)
  Next
 End If
 Print
 WaitKey$
End Function