'******************************************************************************
'Subject: Solution of an m x n linear Diophantine system
' A*x = b using LLL reduction.
'Ref. : G. Havas, B. Majewski, K. Matthews,
' 'Extended gcd and Hermite normal form
' algorithms via lattice basis reduction,'
' Experimental Mathematics 7 (1998), no.2, pp.125-136
'Author : Sjoerd.J.Schaper
'Date : 03-12-2009
'Code : QBasic, PDS, VBdos, FreeBasic (-lang qb)
'******************************************************************************
DEFDBL A-Z
DEFINT I-K, M-N, R-T
DECLARE SUB Inpsys ()
'input linear system
DECLARE SUB Outvcs ()
'output relations
DECLARE SUB Reduce (k, t)
'LLL reduce rows k
DECLARE SUB Swop (k)
'exchange rows k and k-1
DECLARE SUB Negate (t)
'negate rows t
DECLARE FUNCTION chcks% (r)
'check if A*x = b
DECLARE SUB PrntV (r)
'print row P(r,)
DECLARE SUB PrntM (g AS STRING)
'print matrix A(,)
DECLARE SUB PrntG (g AS STRING)
'print Gram coefficients
CONST Verb = 1 ' level 0, 1, 2
CONST na = 99, da = 100 ' alpha = na / da
CONST half = .4999999999999999#
DIM SHARED m1, mn, nx, m, n ' rows & columns A
DIM SHARED tl, tc, kc ' indices
tim# = TIMER
begin:
PRINT
DO: INPUT ; " rows ", g$
LOOP WHILE INSTR(g$, "'")
n = VAL(g$)
IF n < 1 THEN GOTO eind
INPUT " cols ", g$: m = VAL(g$)
IF m < 1 THEN
FOR r = 1 TO n: INPUT g$: NEXT
GOTO begin
END IF
CLS : PRINT
m1 = m + 1: mn = m1 + n: nx = mn + 1
REDIM SHARED la(m, m), d(-1 TO m) ' mu_rs = lambda_rs / d_s
REDIM SHARED a(m, mn), x(m, n) ' Br = d_r / d_r-1
CALL Inpsys
k = 1: tl = 1
WHILE k <= m ' main limiting sequence
t = k - 1
Reduce k, t ' partial size reduction
'
sw = (tc = nx AND kc = nx)
IF sw THEN ' zero rows k-1, k
lk = la(k, t) ' Lovasz condition
db = d(t - 1) * d(k) + lk * lk ' Bk >= (alpha - mu_kt^2) * Bt
sw = db * da < d(t) * d(t) * na ' not satisfied
END IF
IF Verb > 1 THEN
PRINT "k, tc, kc, sw:"; k + 1; tc - m; kc - m; sw
ELSE
tl = tl + 1
END IF
IF sw OR (tc <= kc AND tc < nx) THEN ' force echelon form
Swop k
IF k > 1 THEN k = t ' decrease k
ELSE
FOR i = t - 1 TO 0 STEP -1 ' complete size reduction
Reduce k, i
NEXT i
k = k + 1 ' increase k
END IF
WEND
CALL Outvcs
PRINT : GOTO begin
eind:
g$ = STR$(CSNG(TIMER - tim#))
PRINT : PRINT "Timer: " + g$ + " s"
SYSTEM
FUNCTION chcks% (k)
DIM g AS STRING
sw = -1
FOR r = 0 TO n - 1 ' check v * Ab~
p = 0: FOR s = 0 TO m
p = p + a(k, s) * x(s, r)
NEXT s
IF p THEN
g = "fail: db(" + LTRIM$(STR$(r + 1))
PRINT g + ") ="; p: sw = 0
END IF
NEXT r
chcks = sw
END FUNCTION
SUB Inpsys
DIM g AS STRING, h AS STRING
'
FOR r = 0 TO n - 1 ' 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 m
FOR t = -1 TO 0
j = i
WHILE (INSTR(h, MID$(g, i, 1)) > 0) = t
i = i + 1: WEND
NEXT t
'
a(s, m1 + r) = VAL(MID$(g, j, i - j))
IF s < m THEN
x(s, r) = a(s, m1 + r)
IF i > k THEN a(m, m1 + r) = 0: EXIT FOR
END IF
NEXT s: PRINT
x(m, r) = a(m, m1 + r)
NEXT r
a(m, mn) = 1 ' augment Ab~ with column e_m
'
FOR r = 0 TO m: a(r, r) = 1: NEXT ' prefix standard basis
FOR r = -1 TO m: d(r) = 1: NEXT ' Gram subdeterminants
'
g = "Ab~:": IF Verb THEN g = "I | G:"
PrntM g
END SUB
SUB Negate (t)
FOR s = 0 TO mn
a(t, s) = -a(t, s)
NEXT s
FOR r = 1 TO m
FOR s = 0 TO r - 1
IF r = t OR s = t THEN
la(r, s) = -la(r, s)
END IF
NEXT s
NEXT r
IF Verb > 1 THEN
PRINT "Row"; t + 1; "= -Row"; t + 1
tl = tl + 1: PrntM "loop" + STR$(tl)
END IF
END SUB
SUB Outvcs
DIM g AS STRING
sx = n: IF m1 < n THEN sx = m1
p = 1: FOR s = 0 TO sx - 1
q = a(m - s, m1 + s)
IF q = 0 THEN EXIT FOR
p = p * q
NEXT s
g = "|Det| =" + STR$(p) ' HNF determinant
IF Verb < 2 THEN
IF Verb THEN
g = "P | HNF(G): " + g
ELSE
g = "HNF: " + g
END IF
PrntM g
ELSE
PRINT g
END IF
'
sw = 0: k = 0
FOR r = 0 TO m
IF a(r, m) THEN
IF a(r, m) = 1 THEN
sw = -1
FOR s = m1 TO mn - 1
IF a(r, s) THEN sw = 0: EXIT FOR
NEXT s
END IF
k = r: EXIT FOR
END IF
NEXT r
'
IF sw THEN ' zero HNF row k
FOR s = 0 TO m - 1
IF a(k, s) THEN sw = 0
a(k, s) = -a(k, s)
NEXT s: a(k, s) = -a(k, s)
IF chcks(k) AND (sw = 0) THEN
PRINT "solution:" ' nontrivial only
PrntV k
END IF
ELSE
PRINT "inconsistent"
END IF
FOR r = 0 TO k - 1 ' kernel basis
IF r = 0 THEN PRINT "nullspace:"
PrntV r
NEXT r
PRINT "rank ="; m - k
IF Verb THEN
PrntG "mu:"
IF Verb = 1 THEN PRINT "loop"; tl
END IF
END SUB
SUB PrntG (g AS STRING)
DIM l(m1, m) AS INTEGER, p(m) AS INTEGER
'
PRINT g
FOR s = 0 TO m
p(s) = 1: FOR r = 0 TO m1
IF r < m1 THEN
g = STR$(ABS(la(r, s)))
ELSE
g = STR$(ABS(d(s)))
END IF
l(r, s) = LEN(g) ' store lengths and max. length in column
IF l(r, s) > p(s) THEN p(s) = l(r, s)
NEXT r
NEXT s
'
FOR r = 0 TO m1
IF r = m1 THEN PRINT " /"
FOR s = 0 TO m
PRINT SPACE$(p(s) - l(r, s) + 1);
IF r < m1 THEN
PRINT la(r, s);
ELSE
PRINT d(s);
END IF
NEXT s: PRINT
NEXT r
END SUB
SUB PrntM (g AS STRING)
DIM l(m, mn) AS INTEGER, p(mn) AS INTEGER
'
PRINT g
s0 = m1: sx = mn - 1
IF Verb THEN s0 = 0: sx = mn
'
FOR s = s0 TO sx
p(s) = 1: FOR r = 0 TO m
l(r, s) = LEN(STR$(ABS(a(r, s))))
IF l(r, s) > p(s) THEN p(s) = l(r, s)
NEXT r
NEXT s
'
FOR r = 0 TO m
FOR s = s0 TO sx
IF s = m1 AND s0 = 0 THEN PRINT " |";
PRINT SPACE$(p(s) - l(r, s) + 1); a(r, s);
NEXT s: PRINT
NEXT r
END SUB
SUB PrntV (r)
DIM l(m) AS INTEGER, p AS INTEGER
'
p = 1: FOR s = 0 TO m - 1
l(s) = LEN(STR$(ABS(a(r, s))))
IF l(s) > p THEN p = l(s)
NEXT s
'
PRINT "[";
q = 0: FOR s = 0 TO m - 1
PRINT SPACE$(p - l(s) + 1); a(r, s);
q = q + a(r, s) * a(r, s)
NEXT s: PRINT " ]";
IF Verb THEN PRINT " |v|^2 ="; q;
PRINT
END SUB
SUB Reduce (k, t)
tc = nx: kc = nx ' 1st nonzero row elements
FOR s = m1 TO mn
IF a(t, s) THEN tc = s: EXIT FOR
NEXT s
FOR s = m1 TO mn
IF a(k, s) THEN kc = s: EXIT FOR
NEXT s
'
q = 0
IF tc < nx THEN
IF a(t, tc) < 0 THEN Negate t
q = INT(a(k, tc) / a(t, tc)) ' floor
ELSE
lk = la(k, t)
IF 2 * ABS(lk) > d(t) THEN ' 2|lambda_kt| > d_t
q = INT(lk / d(t) + half) ' round
END IF
END IF
IF q THEN
sx = mn: IF tc = nx THEN sx = m
FOR s = 0 TO sx ' reduce row k
a(k, s) = a(k, s) - q * a(t, s)
NEXT s
la(k, t) = la(k, t) - q * d(t)
FOR s = 0 TO t - 1
la(k, s) = la(k, s) - q * la(t, s)
NEXT s
IF Verb > 1 THEN
PRINT "Row"; k + 1; "-="; q; "x Row"; t + 1
tl = tl + 1: PrntM "loop" + STR$(tl)
END IF
END IF
END SUB
SUB Swop (k)
t = k - 1
FOR s = 0 TO mn ' exchange rows k, k-1
SWAP a(k, s), a(t, s)
NEXT s
FOR s = 0 TO t - 1
SWAP la(k, s), la(t, s)
NEXT s
'
lk = la(k, t)
db = (d(t - 1) * d(k) + lk * lk) / d(t)
FOR r = k + 1 TO m ' update Gram coefficients
lr = la(r, k) ' columns k, k-1 for r > k
la(r, k) = (d(k) * la(r, t) - lk * lr) / d(t)
la(r, t) = (db * lr + lk * la(r, k)) / d(k)
NEXT r
d(t) = db
IF Verb > 1 THEN
PRINT "Swapping Rows"; t + 1; "and"; k + 1
tl = tl + 1: PrntM "loop" + STR$(tl)
END IF
END SUB