'Subject: Solution of a system of linear congruences.
'         Apply the Chinese remainder theorem to p
'         simultaneous congruences a_i*x = b_i Mod M_i.
'Author : Sjoerd.J.Schaper
'Date   : 21-12-2003
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
CONST X = 0, d = 1
r = LargeInit(92, "chineser")
DIM i(91, 3) ' pointer array
CLS : tim# = TIMER: lin = 9

begin:
LOCATE 1, 1
PRINT "   +------------------------------------+"
PRINT "   |  To solve a(i)*x = b(i) Mod M(i),  |"
PRINT "   |  input number of equations p,      |"
PRINT "   |  and p times a(i), b(i), M(i),     |"
PRINT "   |  p = 0 to exit.                    |"
PRINT "   |                                    |"
PRINT "   +------------------------------------+"
LOCATE 6, 7: INPUT ; "p"; ip$: pmx = 0
IF INSTR(ip$, "'") THEN GOTO begin
IF VAL(ip$) < 31 THEN pmx = VAL(ip$)
IF pmx < 1 THEN GOTO eind

p = 0
FOR r = 1 TO 1 + 3 * (pmx - 1) STEP 3
   LOCATE 6, 6: PRINT SPACE$(35); "|"
   LOCATE 6, 6: p = p + 1
   FOR s = 1 TO 3
      PRINT " (" + LTRIM$(STR$(p)) + "," + STR$(s) + ")";
      INPUT ; g$: i(p, s) = r + s
      Readst i(p, s), g$
   NEXT s
   Sets i(p, 3), 1
NEXT r

LOCATE lin, 1
p = 1: ps = -1
Euclid i(p, 1), i(p, 3), d '   i:= a(1)^ -1, M:= M(1)/gcd, d:= gcd(a(1), M(1))
Divd i(p, 2), d, X '                   Xt:= b(1)/gcd
IF NOT Isf(i(p, 2), 0) THEN
   ps = 0: GOTO jump '                 remainder > 0
END IF
Modmlt X, i(p, 1), i(p, 3) '           X:= Xt * i mod M
FOR p = 2 TO pmx
   Mult X, i(p, 1), t0
   Subt i(p, 2), t0 '                  bt:= b(p) - X * a(p)
   Lmul i(p, 1), i(p - 1, 3) '         at:= a(p) * M
   Euclid i(p, 1), i(p, 3), d ' i:= a(p)^ -1, M':= M(p)/gcd, d:= gcd(at, M(p))
   Divd i(p, 2), d, t0 '               t0:= bt / gcd
   IF NOT Isf(i(p, 2), 0) THEN
      ps = 0: EXIT FOR '               remainder > 0
   END IF
   Mult t0, i(p, 1), t1
   Mult t1, i(p - 1, 3), t0 '          t0 * i * M
   Add X, t0 '                         X:= X + (b(p) - X * a(p)) / gcd * i * M
   Lmul i(p, 3), i(p - 1, 3) '         M:= M'* M
   Moddiv X, i(p, 3) '                 X:= X mod M
NEXT p

jump:
IF NOT ps THEN
   Prints "Modulus = 0", 1 '           system unsolvable
   IF p > 1 THEN
      Prints "(" + LTRIM$(STR$(p - 1)) + ") ", 1
   END IF
END IF
IF ps OR p > 1 THEN
   Printn X, "X ð ", "", 2 '           output X
   Printn i(p - 1, 3), " mod ", "", 2 'output M
END IF
Prints "", 1
lin = CSRLIN
GOTO begin

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

