'Subject: Basic library for large-integer arithmetic.
'URL    : http://largeint.sourceforge.net/
'Author : Sjoerd.J.Schaper
'Update : 06-06-2009
'Code   : FreeBasic 0.20.0
#INCLUDE "largeint.bi"

CONST MH = MB \ 2, M1 = &H7FFF% '      bitmask MB - 1
CONST M2 = M1 * MB, SB = NOT M1 '      borrow, sign bit

CONST Asize = 97279 '                  number array size
DIM SHARED AS SHORT n(Asize), uj
'n() is the largeint number array,
'signs and lengths are at positions -1,
'followed by (length + 1) largeint-words
DIM SHARED i() AS INTEGER
'the i()'s are number indices
DIM SHARED AS INTEGER Bezsw, Lognr, Prmnr

FUNCTION Absf (p AS SHORT) AS SHORT export
DIM AS INTEGER s = n(i(p) - 1)
   n(i(p) - 1) = s AND M1
Absf = s AND SB
END FUNCTION

SUB Add (p AS SHORT, q AS SHORT) export
   n(i(q) - 1) XOR= SB
   Subt p, q
   n(i(q) - 1) XOR= SB
END SUB

SUB Bezout (p AS SHORT, q AS SHORT, d AS SHORT) export
   Dup p, t3: Bezsw = -1
   Euclid p, q, d: Bezsw = 0
   IF NOT Isf(d, 1) THEN Ldiv t3, d
   Mult t3, p, t0 '                                p * p^-1)
   Inc t0, -1: Chs t0 '                       (1 -
   Divd t0, q, t3: SWAP i(q), i(t3) '  q^-1 =                / q
END SUB

FUNCTION Bitl (p AS SHORT) AS INTEGER export
DIM AS INTEGER a, t = Getl(p) - 1
   a = n(i(p) + t): t *= LB
   WHILE a
      a SHR= 1: t += 1
   WEND
Bitl = t
END FUNCTION

SUB Boolf (p AS SHORT, q AS SHORT, k AS SHORT) export
DIM AS INTEGER lx = Getl(p) - 1, t = k AND 3
DIM AS INTEGER lm, im, j, o
   IF t THEN
      lm = Getl(q) - 1: im = i(q)
      IF lx < lm THEN SWAP lx, lm: im = i(p)
      IF t > 1 THEN
         FOR j = im + lm + 2 TO im + lx
           n(j) = 0: NEXT j
      ELSE
         lx = lm
      END IF
   END IF
   '
   o = i(q) - i(p)
   SELECT CASE t
   CASE 1
      FOR j = i(p) TO i(p) + lx
         n(j) AND= n(o + j)
      NEXT j
   CASE 2
      FOR j = i(p) TO i(p) + lx
         n(j) XOR= n(o + j)
      NEXT j
   CASE 3
      FOR j = i(p) TO i(p) + lx
         n(j) OR= n(o + j)
      NEXT j
   CASE ELSE
      FOR j = i(p) TO i(p) + lx '      one's complement
         n(j) = M1 AND NOT n(j) '      excepting sign bits
      NEXT j
   END SELECT
   Lftj p, lx
END SUB

SUB Chs (p AS SHORT) export
   n(i(p) - 1) XOR= SB '               signed magnitude
END SUB

FUNCTION Cmp (p AS SHORT, q AS SHORT) AS SHORT export
DIM AS INTEGER lp = n(i(p) - 1), lq = n(i(q) - 1)
DIM AS INTEGER j, o, s, x = SGN((lp AND SB) OR 1)
   '
   IF ((lp XOR lq) AND SB) = 0 THEN '  equal signs
      lp AND= M1
      s = lp - (lq AND M1)
      IF s = 0 THEN '                  equal lengths
         o = i(q) - i(p)
         FOR j = i(p) + lp - 1 TO i(p) STEP -1
            s = n(j) - n(o + j) '      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 ZSTRING PTR) export
DIM s AS STRING: DIM AS INTEGER lp, m, j
DIM AS INTEGER it = i(t2), k = LEN(*g), t = 0
   IF n(it - 1) THEN
      MID$(*g, 1) = "-"
      k -= 1: t = 1
   END IF
   lp = k \ LD: m = k - lp * LD
   IF m = 0 THEN m = LD: lp -= 1
   k += t + 1: t = LD
   FOR j = 0 TO lp
      s = LTRIM$(STR$(n(it + j)))
      IF j = lp THEN t = m
      MID$(*g, k - LEN(s), t) = s '    stuff
      k -= LD
   NEXT j
END SUB

FUNCTION Decf (p AS SHORT) AS INTEGER export
DIM AS INTEGER a, it = i(t2), j = 0
   Dup p, t1: n(it - 1) = Gets(p)
   DO
      n(it + j) = Divint(t1, MD) '     base MD digit
      j += 1
      IF j > uj THEN
         Errorh "conversion overflow Decf"
         j -= 1: EXIT DO
      END IF
   LOOP UNTIL Isf(t1, 0)
   n(it + j) = 0: j -= 1
   a = n(it + j): j *= LD
   DO: a \= 10: j += 1
   LOOP WHILE a
   '
   IF n(it - 1) THEN
      IF n(it) = 0 AND j = 1 THEN '   -0
         n(it - 1) = 0
      ELSE
         j += 1
      END IF
   END IF
Decf = j
END FUNCTION

SUB Divd (p AS SHORT, d AS SHORT, q AS SHORT) export
DIM AS INTEGER iq = i(q) - 1, ip = i(p) - 1
DIM AS INTEGER it = i(d) - 1, lt = n(it) AND M1
DIM AS INTEGER l2, j, k, t, c, c0, c2
DIM nr AS LONGINT
   '
   IF lt = 1 THEN '                    single-word divisor
      IF Isf(d, 0) THEN
         Errorh "div. by zero Sub Divd"
         ERROR 11: EXIT SUB
      END IF
      '
      n(iq) = (n(ip) AND SB) OR 1
      n(iq + 1) = Divint(p, n(it + 1))
      n(iq + 2) = 0: SWAP i(q), i(p)
      n(ip) XOR= n(it) AND SB
   ELSE
      '
      n(iq) = n(ip) XOR n(it)
      l2 = (n(ip) AND M1) - lt
      IF l2 < 0 THEN
         n(iq) = (n(iq) AND SB) OR 1 ' signed zero q
         n(iq + 1) = 0: n(iq + 2) = 0
         EXIT SUB
      END IF
      '
      c2 = CINT(n(it + lt)) SHL LB + n(it + lt - 1)
      iq += 1: it += 1 '               divisor prefix{2}
      '
      ip += 1 - i(t2)
      FOR j = i(t2) - 1 TO i(t2) + l2
        n(j) = n(ip + j): NEXT j '     initial remainder
      '
      FOR k = l2 TO 0 STEP -1 '        long division
         t = i(p) + lt + k
         nr = CLNGINT(n(t)) SHL LB '   current dividend prefix{3}
         nr = (nr + n(t - 1)) SHL LB + n(t - 2)
         c0 = nr \ c2 '                estimate quotient digit k
         IF c0 > M1 THEN c0 = M1
         '
         ip = i(p) - it + k
         t = i(t2) - it + k
         DO WHILE c0
            c = MB
            FOR j = it TO it + lt '    subtract scaled divisor,
               c += M2 - c0 * n(j) + n(ip + j)
               n(t + j) = c AND M1: c SHR= LB
            NEXT j
            IF c = MB THEN '           swap dividend and remainder
               SWAP i(p), i(t2): EXIT DO
            END IF
            c0 -= 1 '                       or decrease quotient
         LOOP
         n(iq + k) = c0
      NEXT k
      Lftj p, lt - 1
      Lftj q, l2
   END IF
END SUB

FUNCTION Divint (p AS SHORT, a AS SHORT) AS SHORT export
DIM AS INTEGER j, c, c0, c2, lp = Getl(p) - 1
   c = 0: c2 = ABS(a)
   FOR j = i(p) + lp TO i(p) STEP -1
      c0 = c SHL LB OR n(j) '          paste carry w/dividend
      n(j) = c0 \ c2 '                 store quotient
      c = c0 - n(j) * c2 '             remainder
   NEXT j
   Lftj p, lp
Divint = c
END FUNCTION

SUB Dup (p AS SHORT, q AS SHORT) export
DIM AS INTEGER j, o = i(p) - i(q)
   FOR j = i(q) - 1 TO i(q) + Getl(p)
      n(j) = n(o + j)
   NEXT j
END SUB

SUB EnQ (q AS SHORT, p AS SHORT, k AS SHORT) export
DIM AS INTEGER j, o, lq = Getl(q), lx = lq + Getl(p)
   IF lx < uj THEN
      o = i(p) - lq
      FOR j = lq TO lx - 1
         n(i(q) + j) = n(o + j)
      NEXT j
      IF k = 0 THEN k = 1
      n(i(q) + lx) = -ABS(k) '         end of record
      n(i(q) - 1) = lx + 1
   ELSE
      Errorh "list overflow Sub EnQ"
      ERROR 6
   END IF
END SUB

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

SUB Euclid (p AS SHORT, q AS SHORT, d AS SHORT) export
DIM AS INTEGER s = Absf(q)
   IF Isf(q, 0) THEN Dup p, d: EXIT SUB
   IF Gets(p) THEN Moddiv p, q
   '
   Letf d, 1: Letf t0, 0
   DO UNTIL Isf(p, 0)
      SWAP i(p), i(q)
      SWAP i(d), i(t0)
      Divd p, q, t1 '                  Euclidean division & inversion
      Mult t0, t1, t2: Subt d, t2 '    r_t = r_t-2 - q_t * r_t-1
   LOOP
   SWAP i(p), i(t0)
   SWAP i(d), i(q)
   '
   IF Gets(p) AND NOT Bezsw THEN Add p, q
   Sets q, s
END SUB

FUNCTION ExQ (p AS SHORT, ByRef k AS SHORT, q AS SHORT, _
 ByRef j AS SHORT) AS SHORT export
DIM AS INTEGER s, o, lq = Getl(q)
   IF j < lq THEN
      ExQ = -1: o = i(p) - j - 1
      DO WHILE j < lq
         s = n(i(q) + j): j += 1
         IF s < 0 THEN EXIT DO '       end of record
         n(o + j) = s
      LOOP
      k = ABS(s): n(o + j) = 0
      n(i(p) - 1) = o + j - i(p)
   ELSE
      ExQ = 0
   END IF
END FUNCTION

SUB Fctl (p AS SHORT, a AS SHORT) export
DIM AS INTEGER r, c, c0
   Letf p, 1: r = 1
   IF a < 3 THEN
      IF a = 2 THEN Inc p, 1
      EXIT SUB
   END IF
   DO
      c = r
      DO: r += 1
         c0 = c: c *= r
      LOOP UNTIL c > MB OR r > a
      Letf t1, c0: Lmul p, t1 '        partial product
   LOOP UNTIL r > a
END SUB

SUB Gcd (p AS SHORT, q AS SHORT, d AS SHORT) export
   Dup p, d: Dup q, t0
   DO UNTIL Isf(t0, 0)
      SWAP i(t0), i(d)
      Divd t0, d, t1
   LOOP
   Sets d, 0
END SUB

FUNCTION Getl (p AS SHORT) AS SHORT export
   Getl = n(i(p) - 1) AND M1
END FUNCTION

FUNCTION Gets (p AS SHORT) AS SHORT export
   Gets = n(i(p) - 1) AND SB
END FUNCTION

FUNCTION Getw (p AS SHORT, j AS SHORT) AS SHORT export
   Getw = n(i(p) + j)
END FUNCTION

FUNCTION Hp2 (p AS SHORT) AS SHORT export
DIM AS INTEGER a = n(i(p) + Getl(p) - 1), a0 = MH
   DO UNTIL a AND a0
      a0 SHR= 1: LOOP '                bitscan L=>R
Hp2 = a0
END FUNCTION

SUB Inc (p AS SHORT, a AS SHORT) export
DIM AS INTEGER j, b = SGN(Gets(p) OR 1) * a
DIM AS INTEGER lp, c = CINT(n(i(p))) + b
   IF c >= 0 AND c < MB THEN
      n(i(p)) = c
   ELSE
      c = MB + b: lp = Getl(p)
      FOR j = i(p) TO i(p) + lp '      add single word
         c += M2 + n(j)
         n(j) = c AND M1: c SHR= LB
      NEXT j
      IF c = MB THEN
         Lftj p, lp
      ELSE
         n(i(p)) = MB - n(i(p))
         n(i(p) + 1) = 0: Chs p
      END IF
   END IF
END SUB

FUNCTION Isf (p AS SHORT, a AS SHORT) AS SHORT export
DIM AS INTEGER s, r = n(i(p) - 1)
   IF a = 0 THEN r = r AND M1 '        mask off sign(0)
   s = (SGN(a) AND SB) OR 1
Isf = r = s AND n(i(p)) = ABS(a)
END FUNCTION

FUNCTION Ismin1 (p AS SHORT, m AS SHORT) AS SHORT export
   Inc p, 1
   Ismin1 = Cmp(p, m) = 0
   Inc p, -1
END FUNCTION

FUNCTION IsPPrm (p AS SHORT, d AS SHORT, w AS SHORT) AS SHORT export
DIM AS INTEGER j, k, a, x = 0, lp = Getl(p)
DIM AS SHORT r = 2
IsPPrm = 0
   IF Gets(p) OR Isf(p, 1) THEN EXIT FUNCTION
   '
   IF lp = 1 OR w >= 0 THEN '          trial divisions
      k = lp * 54
      IF k > 669 THEN k = 669 '        pr(669) = 4999
      n(i(d) - 1) = 1: n(i(d) + 1) = 0
      FOR j = 1 TO k
         Dup p, t1: n(i(d)) = Nxtprm(r)
         IF Divint(t1, n(i(d))) = 0 THEN
            x = Isf(t1, 1): EXIT FOR ' d | N
         END IF
         x = Cmp(d, t1) = 1
         IF x THEN EXIT FOR '          d > sqrt(N)
      NEXT j
      r = Nxtprm(-2)
      IF j <= k THEN
         IF x THEN x = 1
         IsPPrm = x: EXIT FUNCTION
      END IF
   END IF
   '
   a = 2 + INT(RND * 255) '            1-byte random witness
   Dup p, d
   Inc d, -1: k = Odd(d) '             determine N - 1 = (2^ k) * odd(m)
   '
   FOR r = 1 TO ABS(w) '               Miller-Rabin test
      Letf t3, 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
            x = Ismin1(t3, p)
            IF x THEN EXIT FOR
         NEXT j
         IF x = 0 THEN EXIT FOR
      END IF
      a += 1
   NEXT r
IsPPrm = x
END FUNCTION

SUB Isqrt (p AS SHORT, q AS SHORT) export
DIM AS INTEGER j, c, sw, t, d, k, lq, lp = n(i(p) - 1)
   IF lp AND SB THEN
      Errorh "illegal argument Sub Isqrt"
      ERROR 5: EXIT SUB
   END IF
   '
   j = i(p) + lp - 1: c = n(j)
   k = lp AND 1: lq = (lp + 1) \ 2 '   length(sqrt(p))
   IF k = 0 THEN c = c SHL LB + n(j - 1)
   n(i(q) - 1) = 1: n(i(q) + 1) = 0
   n(i(q)) = INT(SQR(c)) '             q' = sqrt(MSB(p))
   IF lq = 1 THEN EXIT SUB
   '
   sw = 0: d = 1
   DO
      IF d SHL 1 < lq THEN
         Lsft q, -d: d SHL= 1 '        q:= doublen(q')
         '
         t = d SHL 1 - k
         c = i(p) - i(t0) + lp - t
         FOR j = i(t0) TO i(t0) + t '  t0:= ShRcopy(p)
           n(j) = n(c + j): NEXT j
         n(i(t0) - 1) = t
      ELSEIF NOT sw THEN
         Lsft q, d - lq: sw = -1
      END IF
      '
      IF sw THEN Dup p, t0
      Divd t0, q, t1 '                 Heron's algorithm
      Add q, t1: Rsft q, 1 '           q' = (q + t0 / q) / 2
      '
      IF sw THEN
         Subt t1, q
         t = n(i(t1) - 1) = 1 AND n(i(t1)) < 2
         IF t THEN EXIT DO '          (t0 / q) - q' = 0 or 1
      END IF
   LOOP
END SUB

FUNCTION IsSqr (p AS SHORT, q AS SHORT) AS SHORT STATIC export
'Square test from Henri Cohen's `Computational Number Theory', 1.7.2
DIM AS INTEGER q11(10), q63(62), q64(63), q65(64), c, k, sw
IF NOT sw THEN
   FOR k = 0 TO 32
      c = k * k
      q11(c MOD 11) = -1
      q63(c MOD 63) = -1
      q64(c AND 63) = -1
      q65(c MOD 65) = -1
   NEXT k: sw = -1
END IF
   IsSqr = 0
   IF q64(n(i(p)) AND 63) THEN
      Letf q, 45045
      Dup p, t0: Divd t0, q, t1
      c = CINT(n(i(t0) + 1)) SHL LB + n(i(t0))
      IF q63(c MOD 63) THEN
         IF q65(c MOD 65) THEN
            IF q11(c MOD 11) THEN
               IF Gets(p) THEN EXIT FUNCTION
               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 SHORT, q AS SHORT) AS SHORT export
DIM AS INTEGER r, x = 1
   '
   Dup p, t0: Dup q, t1
   IF (n(i(t1)) AND 1) = 0 THEN '      even(N)
      IF Isf(t1, 0) THEN '             (a/0)
         Sets t0, 0
         IF NOT Isf(t0, 1) THEN x = 0
         Kronec = x: EXIT FUNCTION
      END IF
      IF (n(i(t0)) AND 1) = 0 THEN '   both even
         Kronec = 0: EXIT FUNCTION
      END IF
      '
      IF Odd(t1) AND 1 THEN '          make odd(N)
         r = n(i(t0)) AND 7
         IF r = 3 OR r = 5 THEN x = -x
      END IF
   END IF
   '
   IF Gets(t0) THEN '                  negative(a)
      IF Mp2(t1, 4) = 3 THEN x = -x '  N mod 4 = 3
      Sets t0, 0 '                     Abs(a)
   END IF
   Sets t1, 0 '                        Abs(N)
   '
   DO
      IF (n(i(t0)) AND 1) = 0 THEN '   even(a)
         IF Isf(t0, 0) THEN
            IF NOT Isf(t1, 1) THEN x = 0
            EXIT DO
         END IF
         '
         IF Odd(t0) AND 1 THEN '       make odd(a)
            r = n(i(t1)) AND 7
            IF r = 3 OR r = 5 THEN x = -x
         END IF
      END IF
      r = n(i(t0)) AND n(i(t1))
      IF r AND 2 THEN x = -x '         reciprocity
      SWAP i(t0), i(t1)
      Divd t0, t1, t3
   LOOP
Kronec = x
END FUNCTION

FUNCTION LargeInit (ix AS SHORT, g AS ZSTRING PTR) AS SHORT export
DIM s AS STRING, t AS INTEGER
LargeInit = 0
   '
   t = Asize \ (5 + ix) '       "5" for four swap-rows + row zero
   IF t < 5 OR t > M1 THEN
      Errorh "out of memory LargeInit fail"
      ERROR 7: EXIT FUNCTION
   END IF
   uj = t
   REDIM i(-4 TO ix) AS INTEGER
   FOR t = -4 TO ix
      i(t) = 1 + (t + 4) * uj '        initialize indices
      Letf t, 0 '                      set numbers to zero
   NEXT t
   '
   s = SPACE$(64): Work 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
   '
   RANDOMIZE TIMER
LargeInit = uj: uj = uj - 2
END FUNCTION

SUB Lcm (p AS SHORT, q AS SHORT, d AS SHORT) export
   Gcd p, q, d
   IF Isf(d, 1) THEN
      Mult p, q, d
   ELSE
      IF Isf(d, 0) THEN EXIT SUB
      Dup p, t0: Divd t0, d, t1
      Mult t1, q, d
   END IF
END SUB

SUB Ldiv (p AS SHORT, d AS SHORT) export
   Divd p, d, t1: SWAP i(p), i(t1)
END SUB

SUB Letf (p AS SHORT, c AS INTEGER) export
DIM AS INTEGER j = 0, c0 = ABS(c)
   DO
      n(i(p) + j) = c0 AND M1
      c0 SHR= LB: j += 1 '             split DWord c
   LOOP WHILE c0: n(i(p) + j) = 0
   n(i(p) - 1) = (SGN(c) AND SB) OR j ' sign, length
END SUB

SUB Lftj (p AS SHORT, k AS SHORT) export
DIM AS INTEGER j, ip = i(p), t = 1
   FOR j = k TO 0 STEP -1
      IF n(ip + j) THEN t = j + 1: EXIT FOR
   NEXT j: n(ip + t) = 0
   n(ip - 1) = (n(ip - 1) AND SB) OR t
END SUB

SUB Lmul (p AS SHORT, q AS SHORT) export
   Mult p, q, t2: SWAP i(p), i(t2)
END SUB

FUNCTION Logf (p AS SHORT) AS DOUBLE export
STATIC lMB AS DOUBLE, sw AS INTEGER
DIM AS INTEGER j, t = Getl(p) - 1, k = t - 3
DIM AS LONGINT nr = n(i(p) + t)
IF NOT sw THEN
   lMB = LOG(MB): sw = -1
END IF
   IF k < 0 THEN k = 0
   FOR j = t - 1 TO k STEP -1
      nr = nr SHL LB + n(i(p) + j)
   NEXT j
Logf = LOG(nr) + k * lMB
END FUNCTION

SUB Lsft (p AS SHORT, r AS SHORT) export
DIM AS INTEGER c, j, k, m, ip = i(p), lp = Getl(p)
   IF n(ip) = 0 AND lp = 1 THEN EXIT SUB
   '
   IF r > 0 THEN
      k = r \ LB: m = r - k * LB
      IF m THEN
         c = 0
         FOR j = ip TO ip + lp
            c SHR= LB
            c OR= CINT(n(j)) SHL m
            n(j) = c AND M1 '          store product
         NEXT j
         IF n(ip + lp) THEN lp += 1
         n(ip + lp) = 0
      END IF
   ELSE
      k = ABS(r)
   END IF
   '
   IF k THEN
      lp += k
      IF lp > uj THEN
         Errorh "overflow Sub Lsft"
         ERROR 6: EXIT SUB
      END IF
      FOR j = ip + lp TO ip + k STEP -1
        n(j) = n(j - k): NEXT j '      ShL words
      FOR k = j TO ip STEP -1
        n(k) = 0: NEXT k
   END IF
   Setl p, lp
END SUB

SUB Modbal (p AS SHORT, m AS SHORT) export
DIM AS INTEGER r, s, z
   Divd p, m, t1
   IF NOT Isf(p, 0) THEN
      s = Absf(p): z = Absf(m)
      Dup p, t1: Lsft t1, 1
      r = Cmp(t1, m)
      IF r = 1 THEN '                  Abs(2p) > m
	 Subt p, m: s XOR= SB '        balance p mod m
      ELSEIF r = 0 THEN
	 s = 0
      END IF
      Sets p, s: Sets m, z
   ELSE
      Sets p, 0
   END IF
END SUB

SUB Moddiv (p AS SHORT, m AS SHORT) export
DIM s AS INTEGER
   Divd p, m, t1
   IF Gets(p) THEN '                   make positive residue
      IF NOT Isf(p, 0) THEN
         s = n(i(m) - 1)
         n(i(m) - 1) = s OR SB '      -Abs(m)
         Subt p, m: Sets m, s
      ELSE
         Sets p, 0
      END IF
   END IF
END SUB

SUB Modmlt (p AS SHORT, q AS SHORT, m AS SHORT) export
   Mult p, q, t2: SWAP i(p), i(t2)
   Divd p, m, t1
END SUB

SUB Modpwr (p AS SHORT, q AS SHORT, m AS SHORT) export
DIM AS INTEGER j, k, s, a, a0, sw = NOT Isf(m, 0)
   IF Isf(q, 0) THEN Letf p, 1: EXIT SUB
   '
   IF sw THEN '                        enable reduction mod m
      s = Absf(m)
      IF Gets(q) THEN
         Euclid p, m, t3 '             handle negative exponent
         IF NOT Isf(t3, 1) THEN
            Errorh "impossible inverse Sub Modpwr"
            ERROR 11: EXIT SUB
         END IF
      ELSE
         Moddiv p, m '                 initial reduction
      END IF
   ELSE
      s = Absf(p)
      IF (n(i(q)) AND 1) = 0 THEN s = 0
   END IF
   '
   k = i(q) + Getl(q) - 1
   FOR j = k TO i(q) STEP -1 '         L=>R binary exponentiation
      a = n(j): a0 = MH '              unsigned bitvector q
      IF j = k THEN
         Dup p, t0
         a0 = Hp2(q) SHR 1 '           handle highest set bit(q)
      END IF
      WHILE a0
         Squ t0, t1
         SWAP i(t0), i(t1) '           square t0
         IF sw THEN Divd t0, m, t1 '   reduce Mod m
         IF a AND a0 THEN '            read bit
            Mult t0, p, t1 '           t0 times base p
            SWAP i(t0), i(t1)
            IF sw THEN Divd t0, m, t1
         END IF
         a0 SHR= 1
      WEND
   NEXT j
   SWAP i(p), i(t0)
   IF sw THEN
      Sets m, s
   ELSE
      Sets p, s
   END IF
END SUB

SUB Modsqu (p AS SHORT, m AS SHORT) export
   Squ p, t2: SWAP i(p), i(t2)
   Divd p, m, t1
END SUB

FUNCTION Mp2 (p AS SHORT, a AS SHORT) AS SHORT export
DIM AS INTEGER m = n(i(p)) AND a - 1
   IF Gets(p) THEN m = a - m AND a - 1
Mp2 = m
END FUNCTION

SUB Mult (p AS SHORT, q AS SHORT, r AS SHORT) export
DIM AS INTEGER ip = i(p), lp = n(ip - 1) AND M1
DIM AS INTEGER iq = i(q), lq = n(iq - 1) AND M1
DIM AS INTEGER j, k, c, c0, ir = i(r), lx = lp + lq - 1
   IF lx > uj THEN
      Errorh "overflow Sub Mult"
      ERROR 6: EXIT SUB
   END IF
   IF lp < lq THEN SWAP ip, iq: SWAP lp, lq
   n(ir - 1) = n(ip - 1) XOR n(iq - 1)
   '
   ip -= ir: c0 = n(iq)
   c = 0: FOR j = ir TO ir + lp '      initialize destination
      c += c0 * n(ip + j)
      n(j) = c AND M1: c SHR= LB
   NEXT j
   FOR k = j TO ir + lx
     n(k) = 0: NEXT k
   '
   iq -= ir
   FOR k = ir + 1 TO ir + lq - 1 '     multiply and add
      ip -= 1: c0 = n(iq + k)
      c = 0: FOR j = k TO k + lp
         c += c0 * n(ip + j) + n(j)
         n(j) = c AND M1: c SHR= LB
      NEXT j
   NEXT k
   Lftj r, lx
END SUB

FUNCTION Nxtprm (ByRef sw AS SHORT) AS INTEGER STATIC export
DIM AS INTEGER c(2), cp(2), cv(2), d(2), b(2)
DIM AS INTEGER fl, r, clp, t, ini '    3-row stack
DIM b4 AS STRING * 4
IF NOT ini THEN
   clp = LOF(Prmnr): t = 0: ini = -1
END IF
IF sw AND 1 THEN
   IF fl THEN
      DO
         d(t) XOR= 6: c(t) += d(t) '   skip multiples of 2 and 3
         IF b(t) = 0 THEN
            b(t) = 30
            IF cp(t) < clp THEN
               GET #Prmnr, cp(t), b4 ' next bitvector in PrimFlgs.bin
               cv(t) = CVL(b4)
               cp(t) += 4
            ELSE
               cv(t) = &H1BF6FDBF '    5-folds excluded
            END IF
         END IF
         r = cv(t) AND 1 '             read bit
         b(t) -= 1: cv(t) SHR= 1
      LOOP UNTIL r
   ELSE
      c(t) += d(t)
      d(t) SHL= 1: fl = c(t) = 5
   END IF
   Nxtprm = c(t)
ELSE
   '
   IF sw < 0 THEN '                    pop
      IF t THEN t -= 1: fl = c(t) > 4
   ELSE
      IF sw AND t < 2 THEN t += 1 '    push
      c(t) = 2: cp(t) = 1: fl = 0
      d(t) = 1: b(t) = 0: sw = -1 '    initialize
   Nxtprm = 2
   END IF
END IF
END FUNCTION

FUNCTION Odd (p AS SHORT) AS SHORT export
DIM AS INTEGER a, lp, ip = i(p), t = 0
Odd = 0
   lp = n(ip - 1) AND M1
   IF lp = 1 THEN
      a = n(ip)
      IF a = 0 THEN a = 1
   ELSE
      DO UNTIL n(ip + t)
         t += 1
      LOOP
      a = n(ip + t): t *= LB
   END IF
   DO UNTIL a AND 1
      a SHR= 1: t += 1
   LOOP
   IF t = 0 THEN EXIT FUNCTION
   IF lp = 1 THEN
      n(ip) = a
   ELSE
      Rsft p, t
   END IF
Odd = t
END FUNCTION

SUB Powr (p AS SHORT, k AS SHORT) export
   Letf t3, ABS(k)
   Letf t2, 0: Modpwr p, t3, t2
END SUB

FUNCTION PrimCeil () AS INTEGER export
DIM b4 AS STRING * 4
DIM AS INTEGER c = 0
   IF LOF(Prmnr) > 3 THEN
      GET #Prmnr, 1, b4
      IF CVL(b4) = &HB76BDBF THEN '    valid primelist
         c = 5 + (LOF(Prmnr) \ 4) * 90
      END IF
   END IF
PrimCeil = c
END FUNCTION

SUB Printf (ByRef f AS STRING, ByRef g AS STRING, ByRef h AS STRING, sw AS SHORT)
DIM s AS STRING, k AS INTEGER
   SELECT CASE sw
   CASE 0
      PRINT f; g; h;
      IF Lognr THEN PRINT #Lognr, f; g; h;
   CASE 1
      PRINT f; g; h
      IF Lognr THEN PRINT #Lognr, f; g; h
   CASE ELSE
      k = LEN(g)
      IF LEFT$(g, 1) = "-" THEN k -= 1
      s = "  [" + LTRIM$(STR$(k)) + "]"
      PRINT f; g; h; s
      IF Lognr THEN PRINT #Lognr, f; g; h; s
   END SELECT
END SUB

SUB Printn (p AS SHORT, ByRef f AS STRING, ByRef h AS STRING, sw AS SHORT) export
DIM AS INTEGER k = Decf(p)
DIM AS STRING g = STRING$(k, "0")
   CnvSt g: Printf f, g, h, sw
END SUB

SUB Printr (p AS SHORT, q AS SHORT, ByRef f AS STRING, e AS SHORT) export
DIM k AS INTEGER, d AS INTEGER
DIM g AS STRING, s AS STRING
   Divd p, q, t1: k = Decf(t1)
   g = STRING$(k, "0"): CnvSt g
   k -= e: d = 0
   IF k > 0 THEN
      g = LEFT$(g, e) '                truncate
      s = " E" + LTRIM$(STR$(k))
   END IF
   IF Isf(p, 0) OR k > 0 THEN '        integer quotient
      Printf f, g, s, 1
   ELSE
      IF g = "0" THEN
	 IF Gets(t1) THEN g = "-0"
	 d = -INT((Logf(q) - Logf(p)) * .43429)
	 IF d < -6 THEN k += d '       enlarge precision
      END IF
      Printf f, g, "", 0
      g = STRING$(ABS(k), "0")
      Ratdec g, p, q
      IF d < -6 THEN
	 g = RIGHT$(g, e - 1) '        shift left
	 s = " E" + LTRIM$(STR$(d))
      END IF
      Printf ".", g, s, 1
   END IF
END SUB

SUB Prints (ByRef f AS STRING, sw AS SHORT) export
   SELECT CASE sw
   CASE 0
      PRINT f;
      IF Lognr THEN PRINT #Lognr, f;
   CASE 1
      PRINT f
      IF Lognr THEN PRINT #Lognr, f
   CASE ELSE
      PRINT f: PRINT
      IF Lognr THEN PRINT #Lognr, f: PRINT #Lognr,
   END SELECT
END SUB

SUB Ratdec (g AS ZSTRING PTR, p AS SHORT, q AS SHORT) export
DIM s AS STRING: DIM AS INTEGER r, j, c = LEN(*g)
DIM AS INTEGER lp = c \ LD, m = c - lp * LD
   IF m = 0 THEN m = LD: lp -= 1
   Letf t1, MD: c = 1
   FOR j = 0 TO lp
      Lmul p, t1 '                         p * MD
      Divd p, q, t0 '                  p:=        mod q
      r = Divint(t0, MD) '             r:=       / q mod MD
      s = LTRIM$(STR$(r))
      IF j < lp THEN
         c += LD
      ELSE
         s = STRING$(LD - LEN(s), "0") + s
         s = LEFT$(s, m): c += m '     last digit
      END IF
      MID$(*g, c - LEN(s), LD) = s '    stuff
   NEXT j
END SUB

SUB Readst (p AS SHORT, gp AS ZSTRING PTR) export
DIM AS INTEGER j, t, sw, k, lp, m
DIM AS STRING g = TRIM$(*gp)
   Letf p, 0: k = LEN(g)
   IF g = "0" OR k = 0 THEN EXIT SUB
   sw = 0: t = 0
   IF LEFT$(g, 1) = "-" THEN
      sw = -1: t = 1: k -= 1
   END IF
   lp = k \ LD: m = k - lp * LD
   IF m = 0 THEN m = LD: lp -= 1
   '
   Letf t1, MD: k = t + 1: t = m
   FOR j = 0 TO lp
      IF j THEN Lmul p, t1: t = LD '   p:= p * MD
      Inc p, VAL(MID$(g, k, t)) '      add base MD digit
      k += t
   NEXT j
   t = Getl(p)
   IF sw THEN n(i(p) - 1) = SB
   Lftj p, t - 1
END SUB

SUB Rndf (p AS SHORT, k AS SHORT) export
DIM AS INTEGER j, r = ABS(k)
DIM AS INTEGER t = r \ LB, m = r - t * LB
   IF m THEN
      t += 1
   ELSE
      m = LB
   END IF
   IF t > uj THEN
      Errorh "overflow Sub Rndf"
      ERROR 6: EXIT SUB
   END IF
   FOR j = i(p) TO i(p) + t - 2
      n(j) = INT(RND * MB)
   NEXT j: r = 1
   FOR j = 2 TO m: r SHL= 1: NEXT
   n(i(p) + t - 1) = r + INT(RND * r)
   IF t = 0 THEN t = 1: n(i(p)) = 0
   n(i(p) - 1) = t: n(i(p) + t) = 0
END SUB

SUB Rsft (p AS SHORT, r AS SHORT) export
DIM AS INTEGER ip = i(p), lp = Getl(p) - 1
DIM AS INTEGER c0, c, j, k, m
   '
   IF r > 0 THEN
      k = r \ LB: m = r - k * LB
      IF m THEN
         c = 0
         FOR j = ip + lp TO ip STEP -1
            c0 = c SHL LB OR n(j)
            c = c0 AND M1
            n(j) = c0 SHR m AND M1 '   store quotient
         NEXT j
         IF n(ip + lp) = 0 AND lp THEN lp -= 1
      END IF
   ELSE
      k = ABS(r)
   END IF
   '
   IF k THEN
      lp -= k
      IF lp < 0 THEN Letf p, 0: EXIT SUB
      FOR j = ip TO ip + lp
        n(j) = n(j + k): NEXT j '      ShR words
   END IF
   n(ip + lp + 1) = 0
   Setl p, lp + 1
END SUB

SUB Setl (p AS SHORT, a AS SHORT) export
   n(i(p) - 1) = (n(i(p) - 1) AND SB) OR ABS(a)
END SUB

SUB Sets (p AS SHORT, a AS SHORT) export
   n(i(p) - 1) = (n(i(p) - 1) AND M1) OR (SGN(a) AND SB)
END SUB

SUB Setw (p AS SHORT, j AS SHORT, a AS SHORT) export
   n(i(p) + j) = ABS(a)
END SUB

SUB Squ (p AS SHORT, q AS SHORT) export
DIM AS INTEGER j, o, k, c, c2, ip = i(p), iq = i(q)
DIM AS INTEGER lp = n(ip - 1) AND M1, lx = lp + lp - 1
   IF lx > uj THEN
      Errorh "overflow Sub Squ"
      ERROR 6: EXIT SUB
   END IF
   '
   j = iq: n(iq - 1) = 0
   FOR k = ip TO ip + lp - 1 '         initialize destination
       c = CINT(n(k)) * n(k)
       n(j) = c AND M1: c SHR= LB
       n(j + 1) = c: j += 2
   NEXT k
   '
   o = ip - iq
   FOR k = 1 TO lp - 1 '               add mixed terms
      o -= 1: iq += 1
      c2 = CINT(n(ip + k)) SHL 1
      c = 0: FOR j = iq TO iq + k - 1
         c += c2 * n(o + j) + n(j)
         n(j) = c AND M1: c SHR= LB
      NEXT j
      c += n(j)
      n(j) = c AND M1
      n(j + 1) += c SHR LB
   NEXT k
   Lftj q, lx
END SUB

SUB Subt (p AS SHORT, q AS SHORT) export
DIM AS INTEGER ix = i(p), lx = n(ix - 1) AND M1
DIM AS INTEGER im = i(q), lm = n(im - 1) AND M1
DIM AS INTEGER j, o, c, s = n(ix - 1) XOR n(im - 1)
   '
   IF s AND SB THEN '                  distinct signs: add
      IF lx < lm THEN SWAP lm, lx: im = ix
      FOR j = im + lm + 1 TO im + lx
        n(j) = 0: NEXT j
      '
      o = i(q) - i(p)
      c = 0: FOR j = i(p) TO i(p) + lx
         c += CINT(n(o + j)) + n(j)
         n(j) = c AND M1: c SHR= LB
      NEXT j
      '
   ELSE '                              subtract
      s = lx - lm
      IF s = 0 THEN '                  equal lengths
         o = im - ix
         FOR j = ix + lx - 1 TO ix STEP -1
            s = n(j) - n(o + j)
            IF s THEN EXIT FOR
         NEXT j
         IF s = 0 THEN
            n(ix - 1) = 1: n(ix) = 0 ' unsigned zero p
            n(ix + 1) = 0: EXIT SUB
         END IF
      END IF
      '
      IF s < 0 THEN '                  p:= -(q - p)
         n(ix - 1) XOR= SB
         SWAP lm, lx: SWAP im, ix
      END IF
      FOR j = im + lm + 1 TO im + lx
        n(j) = 0: NEXT j
      '
      im -= i(p): ix -= i(p)
      c = MB: FOR j = i(p) TO i(p) + lx
         c += M2 - n(im + j) + n(ix + j)
         n(j) = c AND M1: c SHR= LB
      NEXT j
   END IF
   Lftj p, lx
END SUB

SUB Swp (p AS SHORT, q AS SHORT) export
   SWAP i(p), i(q)
END SUB

SUB Term () export
   CLOSE 'all files, then terminate
END SUB

SUB Triald (q AS SHORT, p AS SHORT, c AS INTEGER) export
DIM AS INTEGER cp, k, r, s = n(i(p) - 1)
DIM AS SHORT sw = 2
   n(i(q) - 1) = 0 '                   clear list {q}
   n(i(p) - 1) = s AND M1 '            Abs(p)
   r = n(i(p)) < 2 AND n(i(p) - 1) = 1
   IF r THEN EXIT SUB
   DO
      cp = Nxtprm(sw) '                read primes
      IF cp > c THEN EXIT DO
      Letf t3, cp
      FOR k = 0 TO 32766
         Dup p, t0: Divd t0, t3, t1 '  trial divide p
         IF NOT Isf(t0, 0) THEN EXIT FOR ' nonzero remainder
         SWAP i(p), i(t1) '            running quotient p:= p/dr
      NEXT k
      IF k THEN
         IF k < 32767 THEN
            EnQ q, t3, k '             store pfactor
         ELSE
            Errorh "exp. overflow Sub Triald"
            ERROR 6: EXIT SUB
         END IF
      END IF
      r = Cmp(t3, t1) = 1 '            divisor > Isqrt(p)
   LOOP UNTIL r
   IF r AND NOT Isf(p, 1) THEN
      EnQ q, p, 1: Letf p, 1
   END IF
   cp = Nxtprm(-2)
   Sets p, s
END SUB

SUB Work (ByRef f AS STRING) export
   f = TRIM$(ENVIRON$("LARGEINT"))
   IF f = "" THEN f = CURDIR$
   IF RIGHT$(f, 1) <> "\" THEN f += "\"
END SUB

