How to Solve a System of Linear Equations
in the Integers
The problem is to solve the system of linear equationsA x = b
whereA 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 oft = D-1 U b
are integers, and in that case the solution isx = 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 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 Main
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