'Subject: Gauss-Bareiss fraction free elimination.
'         Bring square A into upper triangular form and solve A*X = B
'Author : Sjoerd.J.Schaper
'Date   : 05-08-2005
'Code   : QuickBasic, PDS, VBdos; FreeBasic extendable
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DECLARE SUB InptM (sw)
'input matrices A & B
DECLARE FUNCTION eliminat% (d)
'Gaussian elimination
DECLARE FUNCTION numstr% (g AS STRING, n)
'convert n to string, return length(g)
DECLARE FUNCTION ratstr% (g AS STRING, n, d)
'convert rational n / d to string
DECLARE SUB Backsubs (d)
'back substitution
DECLARE SUB Solve (n, sw)
'output solution A*X = B
DECLARE SUB PrntM (A())
'print matrices A & B

CONST Mmax = 172 '                      maximum size (-larger in FB)
CONST x = 0, Le = 28 '                  decimal fraction length
DIM SHARED p, q, rk '                   matrix dimensions, rank
DIM g AS STRING

r = LargeInit(Mmax, "gausseli")
CLS : tim# = TIMER: lin = 2

begin:
Prints "", 1: LOCATE lin, 1
DO: INPUT ; " order A ", g
LOOP WHILE INSTR(g, "'")
sw = INSTR(LCASE$(g), "r") > 0 '        suffix "r" for rational output
p = VAL(g): IF p < 1 THEN GOTO eind

INPUT " columns in B ", g
fI = INSTR(UCASE$(g), "I") > 0 '        input "I" for B:= identity matrix
q = VAL(g): IF fI THEN q = p
IF q < 1 THEN
   FOR r = 1 TO p: INPUT g: NEXT
   GOTO begin
END IF

t = 2 + p * (p + q)
IF t > Mmax THEN '                      check size
   FOR r = 1 TO p: INPUT g: NEXT
   Prints "out of memory " + STR$(t), 2
   GOTO begin
END IF

t = 2: p = p - 1
REDIM SHARED A(p, p + q), pt(p)
FOR r = 0 TO p
   FOR s = 0 TO p + q
      t = t + 1: A(r, s) = t '          initialize indices
   NEXT s: pt(r) = r '                  row pointers
NEXT r

InptM fI
PrntM A()
rk = eliminat(x) '                      triangularize
IF rk = p THEN
   Printn x, " Det = ", "", 1
ELSE
   Prints "Rank =" + STR$(rk + 1), 1 '  reduced rank
END IF
Solve x, sw '                           print solution

Prints "", 1
lin = CSRLIN: GOTO begin

eind:
LOCATE lin, 1
g = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + g + " s", 1
CALL Term
SYSTEM

SUB Backsubs (d)
m = 1
   '
   FOR r = rk TO 0 STEP -1 '            solve triangular system
      i = pt(r)
      FOR s = p + 1 TO p + q
	 Letf m, 1
	 Lmul A(i, s), d
	 FOR t = r + 1 TO p '           back substitution
	    Mult A(i, t), A(pt(t), s), t0
	    Mult m, t0, t1
	    Subt A(i, s), t1
	    Lmul m, A(pt(t), t) '       lcm A(t, t)
	 NEXT t
      NEXT s
      Lmul d, A(i, r) '                 return d:= lcm A(r, r)
   NEXT r
END SUB

FUNCTION eliminat% (d)
m = 1
   '
   Letf d, 1
   e = 1: k = p
   FOR s = 0 TO p
      FOR r = s TO k
	 sw = NOT Isf(A(pt(r), s), 0)
	 IF sw THEN EXIT FOR '          find non-zero pivot
      NEXT r

      IF sw = 0 THEN
	 SWAP pt(k), pt(s): k = k - 1
      ELSE
	 IF r > s THEN '                exchange rows,
	    SWAP pt(r), pt(s): e = -e ' change permutation sign
	 END IF
	 '
	 Dup A(pt(s), s), m '           pivot
	 FOR r = s + 1 TO k '           triangularize column s
	    i = pt(r)
	    FOR t = s + 1 TO p + q '    Gaussian elimination
	       Mult A(i, t), m, t0 '    using pseudo-division
	       Mult A(i, s), A(pt(s), t), t1
	       Subt t0, t1
	       Divd t0, d, A(i, t) '    divide by previous pivot
	    NEXT t
	    Letf A(i, s), 0
	 NEXT r: Swp d, m
      END IF
   NEXT s
   '
   Dup A(pt(k), k), d '                 Det(A)
   Sets d, SGN(Gets(d) OR 1) * e '      adjust sign
   '
eliminat = k
END FUNCTION

SUB InptM (sw)
DIM g AS STRING, h AS STRING
   '
   FOR r = 0 TO p '                     input rows A and B
      h = LTRIM$(STR$(r + 1))
      PRINT " row A"; h; " and B"; h;
      INPUT ; " ", g: h = "| " '        valid separators
      '
      i = 1: j = 1
      k = LEN(g): g = g + " # "
      FOR s = 0 TO p + q
	 IF i < k THEN
	    FOR t = -1 TO 0
	       j = i
	       WHILE (INSTR(h, MID$(g, i, 1)) > 0) = t
	       i = i + 1: WEND
	    NEXT t
	 END IF
	 '
	 n = A(r, s)
	 IF sw AND (s > p) THEN '       build matrix I
	    IF r = s - p - 1 THEN
	       Letf n, 1
	    ELSE
	       Letf n, 0
	    END IF
	 ELSE
	    Readst n, MID$(g, j, i - j)
	 END IF
      NEXT s: PRINT
   NEXT r
END SUB

FUNCTION numstr% (g AS STRING, n)
DIM h AS STRING
   lx = LEN(g): k = Decf(n)
   h = STRING$(k, "0"): CnvSt h '       LargeInt conversion
   IF h = "-0" THEN h = "0": k = 1
   g = g + h
numstr = lx + k
END FUNCTION

SUB PrntM (A())
DIM g(p, p + q) AS STRING
DIM l(p, p + q), w(p + q)
   '
   FOR s = 0 TO p + q
      w(s) = 1: FOR r = 0 TO p
	 k = numstr(g(r, s), A(r, s))
	 l(r, s) = k: IF k > w(s) THEN w(s) = k
      NEXT r
   NEXT s
   '
   Prints " A & B:", 1
   FOR r = 0 TO p
      i = pt(r)
      FOR s = 0 TO p + q
	 Prints SPACE$(w(s) - l(i, s) + 2) + g(i, s), 0
	 IF s = p THEN Prints "|", 0
      NEXT s: Prints "", 1
   NEXT r
END SUB

FUNCTION ratstr% (g AS STRING, n, d)
DIM h AS STRING: t = 2
   '
   Divd n, d, t: k = Decf(t)
   g = STRING$(k, "0"): CnvSt g
   IF NOT Isf(n, 0) THEN '              nonzero remainder
      IF g = "0" AND Gets(t) <> 0 THEN
	 g = "-0": k = k + 1
      END IF
      h = STRING$(Le, "0")
      Ratdec h, n, d
      FOR i = Le TO 1 STEP -1
	 IF MID$(h, i, 1) <> "0" THEN EXIT FOR
      NEXT i
      IF i THEN
	 g = g + "." + LEFT$(h, i) '    trim trailing zeroes
	 k = k + 1 + i
      END IF
   END IF
ratstr = k
END FUNCTION

SUB Solve (d, sw)
DIM g(p, q) AS STRING
DIM l(p, q), w(q)
m = 1: t = 2
   '
   Backsubs d
   FOR r = 0 TO rk
      i = pt(r)
      FOR s = 1 TO q
	 n = A(i, p + s)
	 IF sw THEN
	    Dup d, m
	    IF Absf(m) THEN Chs n
	    Gcd m, n, t
	    IF NOT Isf(t, 1) THEN '     common factor
	       Ldiv m, t: Ldiv n, t
	    END IF
	    '
	    k = numstr(g(i, s), n)
	    IF NOT Isf(m, 1) THEN '     rational
	       g(i, s) = g(i, s) + "/"
	       k = numstr(g(i, s), m)
	    END IF
	 ELSE
	    Dup n, m
	    k = ratstr(g(i, s), m, d) ' decimal fraction
	 END IF
	 l(i, s) = k: IF k > w(s) THEN w(s) = k
      NEXT s
      '
      Ldiv d, A(i, r) '                 lcm A(p-r, p-r)
   NEXT r
   '
   Prints " x:", 1
   FOR r = 0 TO rk
      i = pt(r)
      FOR s = 1 TO q
	 Prints SPACE$(w(s) - l(i, s) + 2) + g(i, s), 0
      NEXT s: Prints "", 1
   NEXT r
END SUB

