'******************************************************************************
'Subject: Solution of an m x n linear Diophantine system
' A*x = b using Hermite and Smith normal forms.
'Author : Sjoerd.J.Schaper
'Date : 05-03-2007
'Code : QBasic, PDS, VBdos, FreeBasic (-lang qb)
'******************************************************************************
'This program is copyright (c) 2007 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
DEFDBL A-Z
DEFINT I-K, R-T
DECLARE SUB Inpts (r, s)
'input linear system
DECLARE FUNCTION bezout (x, y)
'Bezout's identity gcd(x,y) = ux + vy
DECLARE SUB Rmult (x(), rx, s, t)
'right multiply mtx X*mtx D
DECLARE SUB Colel (r, s, t)
'column elimination step
DECLARE SUB Lmult (x(), sx, r, t)
'left multiply mtx D*mtx X
DECLARE SUB Rowel (r, s)
'row elimination step
DECLARE SUB Negcol (x(), rx, s)
'negate column x(,s)
DECLARE SUB Redcol (r, t)
'reduce mod A(r, t)
DECLARE FUNCTION dchck% (r, k)
'check divisibility A(r, r)
DECLARE FUNCTION HNF% (r, s)
'put A in Hermite normal form,
'return rank(A)
DECLARE FUNCTION SNF# (r, s)
'put square A in Smith normal form,
'return |det(A)|
DECLARE SUB MultC (x(), rx, sx, y())
'let y:= mtx X*col y
DECLARE SUB Redsol (x(), k)
'find short solution x
DECLARE FUNCTION chcks% (x(), k)
'check if A*x = b
DECLARE SUB PrntC (g$, x(), rx)
'print column x()
DECLARE SUB PrntM (g$, x(), rx, sx)
'print matrix X(,)
CONST Verb = 1 ' print transition matrices V, U
CONST Mxd = 4503599627370496# ' full mantisse 2^ 52
DIM SHARED errsw AS INTEGER ' overflow in Rmult()
DIM SHARED m AS INTEGER, n AS INTEGER ' rows & columns A
DIM SHARED d(1, 1)
tim# = TIMER
begin:
PRINT : errsw = 0
DO: INPUT ; " rows ", g$
LOOP WHILE INSTR(g$, "'")
m = VAL(g$) - 1
IF m < 0 THEN GOTO eind
INPUT " cols ", g$: n = VAL(g$) - 1
IF n < 0 THEN
FOR r = 0 TO m: INPUT g$: NEXT
GOTO begin
END IF
CLS : PRINT
k = n: IF m > k THEN k = m
REDIM SHARED ab(m, n + 1), b(k)
REDIM SHARED v(m, m), a(m, k), u(n, n)
Inpts m, n
PrntM " A:", a(), m, n
PrntC " b:", b(), m
k = HNF(m, n) ' Hermite normal form
IF k < 0 THEN GOTO begin
IF errsw THEN PRINT "overflow": GOTO jump
PrntM " HNF:", a(), m, n
q = SNF(k, m) ' Smith normal form
IF errsw THEN PRINT "overflow": GOTO jump
g$ = " SNF: d(L) =" + STR$(q)
PrntM g$, a(), m, n
IF Verb THEN
PrntM " V:", v(), m, m ' unimodular manipulations
PrntM " U:", u(), n, n
END IF
FOR r = 0 TO m
sw = b(r) <> 0
IF sw THEN EXIT FOR
NEXT r
MultC v(), k, m, b() ' c = V*input vector b
'PrntC " c:", b(), k
FOR r = 0 TO k ' solve system SNF*y = c
q = INT(b(r) / a(r, r)) ' divide by invariant factors
IF b(r) - q * a(r, r) THEN ' nonzero remainder
PRINT "inconsistent": GOTO jump
END IF
b(r) = q
NEXT r
MultC u(), n, k, b() ' solution x = U*y
IF k < n THEN Redsol b(), k + 1 ' rank < n
IF sw THEN PrntC " x:", b(), n
IF k < n THEN
PrntM " nullspace:", u(), n, n - k - 1
END IF
IF chcks(b(), k + 1) THEN ' check solution
PRINT " rank ="; k + 1
END IF
jump:
PRINT : GOTO begin
eind:
g$ = STR$(CSNG(TIMER - tim#))
PRINT : PRINT "Timer: " + g$ + " s"
SYSTEM
FUNCTION bezout (x, y)
DIM c(2, 1)
'
i = 0: j = 1
c(0, i) = x: c(0, j) = y
c(1, i) = 1: c(1, j) = 0
c(2, i) = 0: c(2, j) = 1
'
WHILE c(0, i)
SWAP i, j
q = INT(c(0, i) / c(0, j)) ' Euclidean steps
FOR r = 0 TO 2
c(r, i) = c(r, i) - q * c(r, j)
NEXT r
WEND
'
IF c(0, j) < 0 THEN
FOR r = 0 TO 2
c(r, j) = -c(r, j)
NEXT r
END IF
'
x = c(1, j): y = c(2, j)
bezout = c(0, j)
END FUNCTION
FUNCTION chcks% (x(), k)
FOR r = 0 TO n
FOR s = 0 TO n - k
x(r) = x(r) + u(r, s) ' x:= kernel + x
NEXT s
NEXT r
'
MultC ab(), m, n, x(): s = 0
FOR r = 0 TO m
IF ab(r, n + 1) - x(r) THEN
g$ = "fail: b(" + LTRIM$(STR$(r + 1))
PRINT g$ + ") ="; x(r): s = -1
END IF
NEXT r
chcks = s
END FUNCTION
SUB Colel (r, s, t)
d(0, 0) = a(r, t): d(0, 1) = a(r, s)
d(1, 0) = a(r, s): d(1, 1) = a(r, t)
'
g = bezout(d(0, 0), d(1, 0))
IF g > 1 THEN
d(0, 1) = d(0, 1) / g
d(1, 1) = d(1, 1) / g
END IF
'
Rmult a(), m, t, s
Rmult u(), n, t, s
END SUB
FUNCTION dchck% (i, k)
c = a(i, i)
FOR r = i + 1 TO m ' scan lower right block
FOR s = i + 1 TO k
q = INT(a(r, s) / c)
'
IF a(r, s) - q * c THEN ' nonzero remainder
FOR t = 0 TO k ' add rows i and r
a(i, t) = a(i, t) + a(r, t)
NEXT t
FOR t = 0 TO m
v(i, t) = v(i, t) + v(r, t)
NEXT t
dchck = -1: EXIT FUNCTION
END IF
NEXT s
NEXT r
dchck = 0
END FUNCTION
FUNCTION HNF% (i, j)
t = 0 ' Hermite normal form A
FOR r = 0 TO i
FOR s = t + 1 TO j ' row reduction
IF a(r, s) THEN Colel r, s, t
NEXT s
IF errsw THEN EXIT FUNCTION
'
IF a(r, t) < 0 THEN ' column At = -At
Negcol a(), m, t
Negcol u(), n, t
END IF
'
IF a(r, t) THEN
Redcol r, t: t = t + 1 ' final reductions
END IF
NEXT r
HNF = t - 1
END FUNCTION
SUB Inpts (rx, sx)
DIM g AS STRING, h AS STRING
'
FOR r = 0 TO m ' input A*x = b
h = LTRIM$(STR$(r + 1))
PRINT " row A"; h; " and b"; h;
INPUT ; " ", g: h = "| " ' valid separators
'
i = 1: k = LEN(g): g = g + " # "
FOR s = 0 TO n + 1
FOR t = -1 TO 0
j = i
WHILE (INSTR(h, MID$(g, i, 1)) > 0) = t
i = i + 1: WEND
NEXT t
'
b(r) = VAL(MID$(g, j, i - j))
IF s < n + 1 THEN
ab(r, s) = b(r): a(r, s) = b(r)
IF i > k THEN b(r) = 0: EXIT FOR
END IF
NEXT s: PRINT
ab(r, n + 1) = b(r)
NEXT r
'
FOR r = 0 TO m: v(r, r) = 1: NEXT ' identity matrices
FOR r = 0 TO n: u(r, r) = 1: NEXT
END SUB
SUB Lmult (x(), sx, r, t)
FOR s = 0 TO sx ' rows r and t
c = d(0, 1) * x(t, s) + d(0, 0) * x(r, s)
x(t, s) = d(1, 1) * x(t, s) - d(1, 0) * x(r, s)
x(r, s) = c
NEXT s
END SUB
SUB MultC (x(), rx, sx, y())
DIM c(rx)
'
FOR r = 0 TO rx
FOR s = 0 TO sx ' multiply
c(r) = c(r) + x(r, s) * y(s)
NEXT s
NEXT r
'
FOR r = 0 TO rx: y(r) = c(r): NEXT ' copy back
END SUB
SUB Negcol (x(), rx, s)
FOR r = 0 TO rx
x(r, s) = -x(r, s)
NEXT r
END SUB
SUB PrntC (g$, x(), rx)
DIM l(rx) AS INTEGER, p AS INTEGER
'
p = 1: FOR r = 0 TO rx
l(r) = LEN(STR$(ABS(x(r)))) ' store lengths
IF l(r) > p THEN p = l(r) ' largest entry in column
NEXT r
'
PRINT g$
FOR r = 0 TO rx
PRINT SPACE$(p - l(r) + 1); x(r)
NEXT r
END SUB
SUB PrntM (g$, x(), rx, sx)
DIM l(rx, sx) AS INTEGER, p(sx) AS INTEGER
'
FOR s = 0 TO sx
p(s) = 1: FOR r = 0 TO rx
l(r, s) = LEN(STR$(ABS(x(r, s))))
IF l(r, s) > p(s) THEN p(s) = l(r, s)
NEXT r
NEXT s
'
PRINT g$
FOR r = 0 TO rx
FOR s = 0 TO sx
PRINT SPACE$(p(s) - l(r, s) + 1); x(r, s);
NEXT s: PRINT
NEXT r
END SUB
SUB Redcol (r, t)
FOR s = 0 TO t - 1 ' reduce left block
q = INT(a(r, s) / a(r, t))
'
FOR i = 0 TO m
a(i, s) = a(i, s) - q * a(i, t)
NEXT i
FOR i = 0 TO n
u(i, s) = u(i, s) - q * u(i, t)
NEXT i
NEXT s
END SUB
SUB Redsol (x(), k)
xs = 0
FOR s = 0 TO n - k
FOR r = 0 TO n
SWAP u(r, s), u(r, s + k) ' left-align kernel
'
IF s = 0 THEN
xs = xs + x(r) * x(r) ' sizeof(solution x)
END IF
NEXT r
NEXT s
'
DO
i = -1
FOR s = 0 TO n - k
FOR r = 0 TO n
IF u(r, s) THEN
q = INT(x(r) / u(r, s) + .5)
'
IF q THEN ' trial reductions
z = 0: FOR t = 0 TO n
p = x(t) - q * u(t, s)
z = z + p * p
NEXT t
IF z < xs THEN xs = z: i = r: j = s
END IF
END IF
NEXT r
NEXT s
'
IF i > -1 THEN
q = INT(x(i) / u(i, j) + .5)
FOR r = 0 TO n ' reduce solution
x(r) = x(r) - q * u(r, j)
NEXT r
END IF
LOOP WHILE i > -1
END SUB
SUB Rmult (x(), rx, s, t)
DIM c(3)
FOR r = 0 TO rx ' columns s and t
c(0) = x(r, s) * d(0, 0): c(2) = x(r, t) * d(1, 0)
c(1) = x(r, s) * d(0, 1): c(3) = x(r, t) * d(1, 1)
FOR i = 0 TO 3
IF ABS(c(i)) > Mxd THEN errsw = -1
NEXT i
IF errsw THEN EXIT FOR
x(r, s) = c(2) + c(0)
x(r, t) = c(3) - c(1)
NEXT r
END SUB
SUB Rowel (r, s)
d(0, 0) = a(r, r): d(0, 1) = a(s, r)
d(1, 0) = a(s, r): d(1, 1) = a(r, r)
'
g = bezout(d(0, 0), d(0, 1))
IF g > 1 THEN
d(1, 0) = d(1, 0) / g
d(1, 1) = d(1, 1) / g
END IF
'
Lmult a(), n, r, s
Lmult v(), m, r, s
END SUB
FUNCTION SNF# (i, j)
q = 1 ' Smith normal form A
FOR r = 0 TO i
DO
FOR s = r + 1 TO j ' row reduction
IF a(r, s) THEN Colel r, s, r
NEXT s
IF errsw THEN EXIT FUNCTION
'
sw = 0: FOR s = r + 1 TO j ' column reduction
IF a(s, r) THEN
Rowel r, s: sw = -1
END IF
NEXT s
' check the rest of A
IF sw = 0 THEN sw = dchck(r, i)
LOOP WHILE sw
q = q * a(r, r)
NEXT r
SNF = q
END FUNCTION