Use Ctrl + or – to enlarge or reduce text size. |

in the Integers

The problem is to solve the system of linear equations

A **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 of

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