'Subject: Basic library for large-integer arithmetic.
'         Static array version.
'Author : Sjoerd.J.Schaper - vspickelen@zonnet.nl
'Date   : 02-02-2004
'Code   : QuickBasic, PDS, VBdos
REM $INCLUDE: 'largeint.bi'
REM $STATIC
DIM SHARED n(-2 TO uj, -4 TO ui) AS INTEGER
DIM SHARED i(-4 TO ui) AS INTEGER
'n() holds the numbers, n(-2, ..) their lengths,
'n(-1, ..) the signs, the i()'s are column indices.
DIM SHARED ln(-2 TO uj) AS LONG
'Used in Decf%, CnvSt and Mult only.
DIM SHARED Lognr AS INTEGER, Prmnr AS INTEGER

SUB Add (p AS INTEGER, q AS INTEGER)
   n(-1, i(q)) = -n(-1, i(q))
   Subt p, q
   n(-1, i(q)) = -n(-1, i(q))
END SUB

SUB Andf (p AS INTEGER, q AS INTEGER)
   lx = n(-2, i(p)) - 1
   lm = n(-2, i(q)) - 1
   IF lx < lm THEN SWAP lx, lm
   FOR j = 0 TO lm
      n(j, i(p)) = n(j, i(p)) AND n(j, i(q))
   NEXT j
   Lftj p, lm
   IF n(-1, i(p)) = -1 THEN
      n(-1, i(p)) = n(-1, i(q))
   END IF
END SUB

SUB Chs (p AS INTEGER)
   r = CLNG(n(0, i(p))) + n(1, i(p)) = 0
   IF NOT (r AND (n(-2, i(p)) = 1)) THEN
      n(-1, i(p)) = -n(-1, i(p))
   END IF
END SUB

FUNCTION Cmp% (p AS INTEGER, q AS INTEGER)
   x = n(-1, i(p))
   IF x = n(-1, i(q)) THEN '           equal signs
      s = n(-2, i(p)) - n(-2, i(q))
      IF s = 0 THEN '                  equal lengths
	 FOR j = n(-2, i(p)) - 1 TO 0 STEP -1
	    s = n(j, i(p)) - n(j, i(q)) ' L=>R check
	    IF s THEN EXIT FOR
	 NEXT j
	 IF s = 0 THEN x = 0
      END IF
      IF s < 0 THEN x = -x
   END IF
Cmp = x
END FUNCTION

SUB CnvSt (g AS STRING)
   k = LEN(g): t = 0
   IF ln(-1) = -1 THEN
      MID$(g, 1) = "-"
      t = 1: k = k - 1
   END IF
   lp = k \ LD: m = k - lp * LD
   IF m = 0 THEN m = LD: lp = lp - 1
   k = k + t + 1: t = LD
   FOR j = 0 TO lp
      s$ = LTRIM$(STR$(ln(j)))
      IF j = lp THEN t = m
      MID$(g, k - LEN(s$), t) = s$ '   stuff
      k = k - LD
   NEXT j
END SUB

SUB Copyf (p AS INTEGER, q AS INTEGER)
   FOR j = -2 TO n(-2, i(p))
      n(j, i(q)) = n(j, i(p))
   NEXT j
END SUB

SUB Dcr (p AS INTEGER, a AS INTEGER)
   a = ABS(a)
   c = CLNG(n(0, i(p))) - n(-1, i(p)) * a
   IF c > 0 AND c < MB THEN
      n(0, i(p)) = c
   ELSE
      IF c = 0 THEN
	 n(0, i(p)) = 0
	 IF Isf(p, 0) THEN n(-1, i(p)) = 1
      ELSE
	 Letf t2, CLNG(a)
	 Subt p, t2
      END IF
   END IF
END SUB

FUNCTION Decf% (p AS INTEGER)
   ln(-1) = n(-1, i(p))
   Letf t3, MY
   Copyf p, t0: j = 0
   DO
      Divd t0, t3, t1 '                decimal base t3
      c = n(0, i(t0)) '                unrolled single digit division
      Divd t1, t3, t0
      ln(j) = n(0, i(t1)) * MY + c
      j = j + 1
   LOOP UNTIL Isf(t0, 0)
   ln(j) = 0: ln(-2) = j
   j = j - 1: c = ln(j): j = j * LD
   DO
      c = c \ 10: j = j + 1
   LOOP WHILE c
   IF ln(-1) = -1 THEN j = j + 1
Decf = j
END FUNCTION

FUNCTION DecLSB% (p AS INTEGER)
   Letf t3, MY
   Copyf p, t0
   Divd t0, t3, t1
DecLSB = n(-1, i(p)) * n(0, i(t0))
END FUNCTION

SUB Divd (p AS INTEGER, d AS INTEGER, q AS INTEGER)
   s = n(-1, i(p)) * n(-1, i(d))
   n(-1, i(q)) = s
   l0 = n(-2, i(p)): l1 = n(-2, i(d))
   m = l0 - l1
   IF m < 0 THEN
      Letf q, 0
      n(-1, i(q)) = s: EXIT SUB
   END IF
   IF l1 = 1 THEN
      IF Isf(d, 0) THEN
	 Errorh "div. by zero in Sub Divd"
	 Letf q, 0: ERROR 11: EXIT SUB
      END IF
      t = l0 - 1: c = 0
      cs = n(0, i(d)) '                short divisor
      FOR j = t TO 0 STEP -1
	 SftL c, LB
	 c0 = c OR n(j, i(p)) '        paste carry w/dividend
	 c = c0 \ cs: n(j, i(q)) = c ' save quotient
	 c = c0 - c * cs '             remainder
      NEXT j
      n(-2, i(p)) = 1
      n(0, i(p)) = c
      n(1, i(p)) = 0: Lftj q, t
   ELSE
      cs = n(l1 - 2, i(d)) + n(l1 - 1, i(d)) * MB
      FOR j = -1 TO m '                MSB(divisor)
	 n(j, i(t2)) = n(j, i(p))
      NEXT j
      '
      FOR k = m TO 0 STEP -1
	 t = l1 + k
	 c2 = EstQ(n(t, i(p)), cs)
	 'machine language replacement for:
	 'nr# = n(t, i(p)) '           MSB(intermediate dividend)
	 'FOR j = t - 1 TO t - 2 STEP -1
	 '   nr# = nr# * MB + n(j, i(p))
	 'NEXT j
	 'c2 = INT(nr# / cs) '         estimate partial quotient
	 'IF c2 > M1 THEN c2 = M1

	 n(k, i(q)) = c2
	 IF c2 THEN
	    DO
	       c = MB
	       FOR j = 0 TO l1 '       subtract scaled divisor,
		  c = c + M2 - c2 * n(j, i(d)) + n(j + k, i(p))
		  n(j + k, i(t2)) = c AND M1: SftR c, LB
	       NEXT j
	       IF c = MB THEN
		  SWAP i(p), i(t2) '   swap dividend & remainder
		  EXIT DO
	       ELSE
		  IF c2 > 1 THEN '     or decrease quotient
		     c2 = c2 - 1
		     n(k, i(q)) = c2
		  ELSE
		     n(k, i(q)) = 0: EXIT DO
		  END IF
	       END IF
	    LOOP
	 END IF
      NEXT k
      Lftj p, l1 - 1
      Lftj q, m
   END IF
END SUB

SUB Errorh (g AS STRING)
   Printf CHR$(13) + CHR$(10), "Error: ", g, 1
END SUB

SUB Euclid (p AS INTEGER, q AS INTEGER, d AS INTEGER)
   IF Isf(q, 0) THEN
      Copyf p, d: EXIT SUB
   END IF
   s = n(-1, i(q)): n(-1, i(q)) = 1
   IF n(-1, i(p)) = -1 THEN Moddiv p, q
   IF Isf(p, 0) THEN
      SWAP i(q), i(d)
      Letf q, 1: EXIT SUB
   END IF
   '
   Letf d, 1: Letf t0, 0
   DO
      SWAP i(p), i(q)
      SWAP i(d), i(t0)
      Divd p, q, t1 '                  Euclidean division & inversion
      Mult t1, t0, t2: Subt d, t1 '    r_t = r_t-2 - q_t * r_t-1
   LOOP UNTIL Isf(p, 0)
   SWAP i(p), i(t0)
   SWAP i(q), i(d)
   IF n(-1, i(p)) = -1 THEN Add p, q
   n(-1, i(q)) = s
   n(-1, i(d)) = 1
END SUB

SUB Fctl (p AS INTEGER, a AS INTEGER)
   IF a < 3 THEN
      IF a < 0 THEN
	 Errorh "illegal argument Sub Fctl"
	 Letf p, 1: ERROR 5: EXIT SUB
      END IF
      IF a < 2 THEN a = 1
      Letf p, CLNG(a): EXIT SUB
   END IF
   Letf p, 1: r = 1
   DO
      nr# = r
      DO
	 r = r + 1
	 c = CLNG(nr#): nr# = nr# * r
      LOOP UNTIL nr# > Fx OR r > a
      Letf t1, c '                     partial product
      IF n(-2, i(p)) + n(-2, i(t1)) >= uj THEN
	 Errorh "overflow in Sub Fctl"
	 Letf p, 1: ERROR 6: EXIT SUB
      END IF
      Mult p, t1, t2
   LOOP UNTIL r > a
END SUB

FUNCTION Gete% (p AS INTEGER, j AS INTEGER)
   Gete = n(j, i(p))
END FUNCTION

FUNCTION Getl% (p AS INTEGER)
   Getl = n(-2, i(p))
END FUNCTION

FUNCTION Gets% (p AS INTEGER)
   Gets = n(-1, i(p))
END FUNCTION

FUNCTION Hp2% (p AS INTEGER)
   a = n(n(-2, i(p)) - 1, i(p))
   a0 = MH
   DO UNTIL a AND a0
      a0 = a0 \ 2
   LOOP
   Hp2 = a0
END FUNCTION

SUB Inc (p AS INTEGER, a AS INTEGER)
   a = ABS(a)
   c = CLNG(n(0, i(p))) + n(-1, i(p)) * a
   IF c > 0 AND c < MB THEN
      n(0, i(p)) = c
   ELSE
      IF c = 0 THEN
	 n(0, i(p)) = 0
	 IF Isf(p, 0) THEN n(-1, i(p)) = 1
      ELSE
	 Letf t2, CLNG(a): Add p, t2
      END IF
   END IF
END SUB

FUNCTION Init% (ix AS INTEGER, g AS STRING)
   IF ix > ui THEN
      PRINT "fail: Ubound(n, 2) ui too small,"
      PRINT "recompile with a larger value."
      END
   END IF
   FOR t = -4 TO ui '                  initialize indices,
      i(t) = t: Letf t, 0 '             set numbers to zero
   NEXT t
   '
   s$ = LTRIM$(RTRIM$(ENVIRON$("LARGEINT")))
   IF s$ = "" THEN s$ = ".\" '         working directory
   IF RIGHT$(s$, 1) <> "\" THEN s$ = s$ + "\"
   '
   Lognr = 0: Prmnr = FREEFILE
   OPEN s$ + "PrimFlgs.bin" FOR BINARY SHARED AS Prmnr
   IF LEN(g) THEN
      Lognr = FREEFILE
      OPEN s$ + g + ".log" FOR OUTPUT AS Lognr
   END IF
   '
   CALL Ffix
   RANDOMIZE TIMER
Init = -1
END FUNCTION

FUNCTION Isf% (p AS INTEGER, a AS INTEGER)
   r = CLNG(n(1, i(p))) + n(0, i(p)) = ABS(a)
   s = SGN(a): IF a = 0 THEN s = n(-1, i(p))
   r = r AND (n(-1, i(p)) = s)
   Isf = r AND (n(-2, i(p)) = 1)
END FUNCTION

FUNCTION Ismin1% (p AS INTEGER, m AS INTEGER)
   IF n(-2, i(p)) > n(-2, i(m)) THEN
      IF n(-1, i(p)) = 1 THEN
	 Ismin1 = 0: EXIT FUNCTION
      END IF
   END IF
   Inc p, 1
   Ismin1 = Cmp(p, m) = 0
   Dcr p, 1
END FUNCTION

FUNCTION IsPPrm% (p AS INTEGER, d AS INTEGER, w AS INTEGER)
   x = -1
   IF n(-1, i(p)) = -1 OR Isf(p, 1) THEN
      x = 0: sw = -1: GOTO jmp
   END IF
   sw = n(-1, i(d)) = -128 AND n(-2, i(p)) > 1
   IF NOT sw THEN '                       try small divisors
      Letf d, 2: a = 2: df = 1
      k = n(-2, i(p)) * 51
      IF k > 1547 THEN k = 1547
      FOR j = 1 TO k
	 Copyf p, t0
	 Divd t0, d, t1
	 IF Isf(t0, 0) THEN '             d | N
	    x = Isf(t1, 1)
	    sw = 0: GOTO jmp
	 END IF
	 IF Cmp(d, t1) = 1 THEN GOTO jmp 'd > sqrt(N)
	 IF sw THEN
	    DO
	       df = 6 - df: a = a + df '  sieve multiples of 2, 3 and 5
	    LOOP UNTIL a MOD 5
	 ELSE
	    a = a + df: df = df * 2: sw = a = 5
	 END IF
	 n(0, i(d)) = a
      NEXT j
   END IF
   '
   Copyf p, d '                        Miller-Rabin test
   Dcr d, 1: k = 0 '                   determine N - 1 = (2^ k) * odd(m)
   DO
      k = k + 1: Rsft d, 1
   LOOP UNTIL n(0, i(d)) AND 1
   a = 5 + INT(RND * 313) '            random witness
   r = a MOD 6
   IF r < 3 THEN
      a = a + 1 - r: df = 2
   ELSE
      a = a + 5 - r: df = 4
   END IF
   FOR r = 1 TO w
      Letf t3, CLNG(a)
      Modpwr t3, d, p '                if N is prime then
      x = Isf(t3, 1) OR Ismin1(t3, p) 'a^ m Mod N = 1 or
      IF x = 0 THEN '         (a^ m)^ (2^ j) Mod N = -1
	 FOR j = 1 TO k - 1 '   for some j in [0, k - 1]
	    Modsqu t3, p
	    IF Ismin1(t3, p) THEN x = -1: EXIT FOR
	 NEXT j
	 IF x = 0 THEN GOTO jmp
      END IF
      DO
	 df = 6 - df: a = a + df
      LOOP UNTIL a MOD 5
   NEXT r
jmp:
   IF x OR sw THEN Letf d, 1 '         Else return small factor d
IsPPrm = x
END FUNCTION

SUB Isqrt (p AS INTEGER, q AS INTEGER)
   IF n(-1, i(p)) = -1 THEN
      Errorh "illegal argument Isqrt"
      Letf q, 0: ERROR 5: EXIT SUB
   END IF
   IF Isf(p, 0) OR Isf(p, 1) THEN
      Copyf p, q: EXIT SUB
   END IF
   k = n(-2, i(p)) - 1
   s = n(k, i(p)): r! = 32 '           argument a in p
   DO
      tr! = s / r!
      r! = (r! + tr!) / 2 '            compute MSB(sqrt)
   LOOP UNTIL ABS(tr! - r!) < 1
   l1 = 1 + k \ 2
   n(-2, i(q)) = l1 '                  compose seed x0 in q
   n(-1, i(q)) = 1
   FOR j = 0 TO l1: n(j, i(q)) = 0: NEXT
   s = (k AND 1) * 180 '               int(sqrt(MB))
   n(l1 - 1, i(q)) = CINT(r! * (1 + s))
   '
   DO '                                Heron's algorithm
      Copyf p, t0
      Divd t0, q, t1 '                 a / x0
      Add q, t1: Rsft q, 1 '           x1:= (x0 + a/x0) / 2
      Subt t1, q
      s = CLNG(n(1, i(t1))) + n(0, i(t1)) < 2
      s = s AND (n(-1, i(t1)) = n(-2, i(t1)))
   LOOP UNTIL s '                      If a/x0 - x1 = 0 Or 1 Then Exit
END SUB

FUNCTION IsSqr% (p AS INTEGER, q AS INTEGER) STATIC
'Square test from Henri Cohen's `Computational Number Theory', 1.7.2
IF NOT sw THEN
   DIM q11(10), q63(62), q64(63), q65(64)
   FOR k = 0 TO 32
      k2 = k * k
      q11(k2 MOD 11) = -1
      q63(k2 MOD 63) = -1
      q64(k2 AND 63) = -1
      q65(k2 MOD 65) = -1
   NEXT k: sw = -1
END IF
   IsSqr = 0
   IF q64(n(0, i(p)) AND 63) THEN
      Copyf p, t0
      Letf q, 45045
      Divd t0, q, t1
      c = n(0, i(t0)) + n(1, i(t0)) * MB
      IF q63(c MOD 63) THEN
	 IF q65(c MOD 65) THEN
	    IF q11(c MOD 11) THEN
	       Isqrt p, q: Squ q, t1 ' p is possibly square
	       IsSqr = Cmp(p, t1) = 0
	    END IF
	 END IF
      END IF
   END IF
END FUNCTION

FUNCTION Kronec% (p AS INTEGER, q AS INTEGER, d AS INTEGER)
   x = 1
   Copyf p, t0
   Copyf q, t1
   IF n(-1, i(t1)) = -1 THEN '         negative(N)
      n(-1, i(t1)) = 1
      x = n(-1, i(t0))
   END IF
   IF (n(0, i(t1)) AND 1) = 0 THEN '   even(N)
      IF Isf(t1, 0) THEN
	 SWAP i(d), i(t0)
	 x = -128: GOTO jump '         (a/0) undefined
      END IF
      IF (n(0, i(t0)) AND 1) = 0 THEN
	 Letf d, 2
	 IF Isf(t0, 0) THEN SWAP i(d), i(t1)
	 x = 0: GOTO jump '            least common divisor 2
      END IF
      IF Isf(t1, 2) THEN
	 Letf d, 1: GOTO jump
      END IF
      r = 0
      DO
	 r = r + 1: Rsft t1, 1
      LOOP UNTIL n(0, i(t1)) AND 1 '   make odd(N)
      IF r AND 1 THEN
	 r = n(0, i(t0)) AND 7
	 IF r = 3 OR r = 5 THEN x = -x
      END IF
   END IF
   IF n(-1, i(t0)) = -1 THEN '         negative(a)
      n(-1, i(t0)) = 1
      IF (n(0, i(t1)) AND 3) = 3 THEN x = -x
   END IF
   IF Isf(t1, 1) THEN
      SWAP i(d), i(t1)
      GOTO jump
   END IF
   '
   DO
      Divd t0, t1, d
      IF (n(0, i(t0)) AND 1) = 0 THEN
	 IF Isf(t0, 0) THEN
	    SWAP i(d), i(t1) '         N | a