largeint.bas

一个大整数计算的小程序,可以进行大数的计算,可以对素数进行测试

	    x = 0: EXIT DO '           return N in d
	 END IF
	 r = 0
	 DO
	    r = r + 1: Rsft t0, 1
	 LOOP UNTIL n(0, i(t0)) AND 1
	 IF r AND 1 THEN
	    r = n(0, i(t1)) AND 7
	    IF r = 3 OR r = 5 THEN x = -x
	 END IF
      END IF
      IF Isf(t0, 1) THEN '            (1/N) = 1 for all N
	 SWAP i(d), i(t0)
	 EXIT DO
      END IF
      r = (n(0, i(t0)) AND 3) = 3
      s = (n(0, i(t1)) AND 3) = 3
      IF r AND s THEN x = -x
      SWAP i(t0), i(t1)
   LOOP
jump:
Kronec = x
END FUNCTION

SUB Letf (p AS INTEGER, c AS LONG)
   n(-1, i(p)) = SGN(c + .5)
   j = 0: c0 = ABS(c) '                SGN(0) = 1
   DO
      n(j, i(p)) = c0 AND M1
      SftR c0, LB: j = j + 1 '         split DWord c
   LOOP WHILE c0
   n(-2, i(p)) = j
   n(j, i(p)) = 0
END SUB

SUB Lftj (p AS INTEGER, k AS INTEGER)
   t = 1
   FOR j = k TO 0 STEP -1
      IF n(j, i(p)) THEN t = j + 1: EXIT FOR
   NEXT j
   n(-2, i(p)) = t
   n(t, i(p)) = 0
END SUB

SUB Lsft (p AS INTEGER, r AS INTEGER)
   k = r \ LB: m = r - k * LB
   t = n(-2, i(p))
   IF m > 0 THEN
      c = 0
      FOR j = 0 TO t
	 c0 = n(j, i(p))
	 SftL c0, m: c = c OR c0 '     ShL element, paste carry
	 n(j, i(p)) = c AND M1 '       mask off high bits
	 SftR c, LB
      NEXT j
      IF n(t, i(p)) THEN t = t + 1
      n(t, i(p)) = 0
   END IF
   IF k > 0 THEN
      t = t + k
      IF t >= vj THEN
	 Errorh "overflow in Sub Lsft"
	 Letf p, 1: ERROR 6: EXIT SUB
      END IF
      IF Isf(p, 0) THEN EXIT SUB
      FOR j = t TO k STEP -1
	 n(j, i(p)) = n(j - k, i(p)) ' ShL array
      NEXT j
      FOR j = 0 TO k - 1: n(j, i(p)) = 0: NEXT
   END IF
   n(-2, i(p)) = t
END SUB

SUB Modbal (p AS INTEGER, m AS INTEGER)
   s = n(-1, i(p))
   z = n(-1, i(m))
   n(-1, i(p)) = 1
   n(-1, i(m)) = 1
   Divd p, m, t1
   Copyf m, t1: Rsft t1, 1
   IF NOT Isf(t1, 0) THEN
      r = Cmp(p, t1)
      IF r = 1 THEN '                  Abs(p) > m \ 2
	 Subt p, m: s = -s '           balance p mod m
      ELSEIF r = 0 THEN
	 IF (n(0, i(m)) AND 1) = 0 THEN s = z
      END IF
   END IF
   IF Isf(p, 0) THEN s = 1
   n(-1, i(p)) = s
   n(-1, i(m)) = z
END SUB

SUB Moddiv (p AS INTEGER, m AS INTEGER)
   Divd p, m, t1
   IF n(-1, i(p)) = -1 THEN '          make positive residue
      IF Isf(p, 0) THEN
	 n(-1, i(p)) = 1
      ELSE
	 s = n(-1, i(m))
	 n(-1, i(m)) = -1
	 Subt p, m
	 n(-1, i(m)) = s
      END IF
   END IF
END SUB

SUB Modmlt (p AS INTEGER, q AS INTEGER, m AS INTEGER)
   Mult p, q, t1
   Divd p, m, t1
END SUB

SUB Modpwr (p AS INTEGER, q AS INTEGER, m AS INTEGER)
   IF Isf(q, 0) THEN
      Letf p, 1: EXIT SUB
   END IF
   sw = NOT Isf(m, 0) '                enable reduction mod m
   IF sw THEN
      s = n(-1, i(m)): n(-1, i(m)) = 1
      Moddiv p, m '                    initial reduction
   ELSE
      s = n(-1, i(p)): n(-1, i(p)) = 1
      IF (n(0, i(q)) AND 1) = 0 THEN s = 1
   END IF
   '
   k = n(-2, i(q)) - 1
   FOR j = k TO 0 STEP -1 '           L=>R binary exponentiation
      a = n(j, i(q)): a0 = MH '        unsigned bitvector q
      IF j = k THEN
	 Copyf p, t0: a0 = Hp2(q) '    handle highest set bit(q)
	 a = a AND (a0 - 1): a0 = a0 \ 2
      END IF
      DO 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
	    Mult t0, p, t1 '           t0 times base p
	    IF sw THEN Divd t0, m, t1
	 END IF
	 a = a AND (a0 - 1): a0 = a0 \ 2 ' read bit
      LOOP
   NEXT j
   SWAP i(p), i(t0)
   IF sw THEN
      n(-1, i(m)) = s
   ELSE
      n(-1, i(p)) = s
   END IF
END SUB

SUB Modsqu (p AS INTEGER, m AS INTEGER)
   Squ p, t1
   SWAP i(p), i(t1)
   Divd p, m, t1
END SUB

SUB Mult (p AS INTEGER, q AS INTEGER, r AS INTEGER)
   l0 = n(-2, i(p)): l1 = n(-2, i(q))
   lx = l0 + l1 - 1
   IF lx >= vj THEN
      Errorh "overflow in Sub Mult"
      Letf r, 1: ERROR 6: EXIT SUB
   END IF
   sw = l0 < l1
   IF sw THEN
      SWAP i(p), i(q): SWAP l0, l1
   END IF
   '
   c = 0: c2 = n(0, i(q))
   FOR j = 0 TO l0
      c = c + c2 * n(j, i(p)) '        initialize destination
      ln(j) = c AND M1: SftR c, LB
   NEXT j
   FOR j = l0 + 1 TO lx: ln(j) = 0: NEXT
   '
   FOR m = 1 TO l1 - 1 STEP 2
      FOR t = 0 TO 1
	 k = t + m: c = n(k, i(q))
	 FOR j = 0 TO l0 - 1 '         multiply,
	    ln(j + k) = c * n(j, i(p)) + ln(j + k)
	 NEXT j
	 IF k = l1 - 1 THEN EXIT FOR
      NEXT t
      c = 0
      FOR j = m TO k + l0
	 c = c + ln(j) '               normalize,
	 ln(j) = c AND M1: SftR c, LB
      NEXT j
   NEXT m
   '
   FOR j = 0 TO lx '                   copy into n(),
      n(j, i(r)) = ln(j)
   NEXT j
   Lftj r, lx '                        and resize
   n(-1, i(r)) = n(-1, i(p)) * n(-1, i(q))
   IF sw THEN SWAP i(p), i(q)
   SWAP i(p), i(r)
END SUB

FUNCTION Nxtprm& (sw AS INTEGER) STATIC
DIM b4 AS STRING * 4
IF sw = 0 THEN
   c = 2: dc = 1: fl = 0: j = 0 '     initialize
   cpl = LOF(Prmnr): cp = 1
ELSE
   IF fl THEN
      DO
	 dc = 6 - dc: c = c + dc '    skip multiples of 2 and 3
	 IF j = 0 THEN
	    j = 30
	    IF cp < cpl THEN
	       GET #Prmnr, cp, b4 '   next bitvector in PrimFlgs.bin
	       cb = CVL(b4): cp = cp + 4
	    ELSE
	       cb = 469171647 '       5-folds excluded
	    END IF
	 END IF
	 j = j - 1
	 r = cb AND 1: cb = cb \ 2 '  read bit
      LOOP UNTIL r
   ELSE
      c = c + dc: dc = dc * 2: fl = c = 5
   END IF
END IF
Nxtprm = c
END FUNCTION

FUNCTION PrimCeil&
DIM b4 AS STRING * 4
   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 (f AS STRING, g AS STRING, h AS STRING, sw 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 = 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 INTEGER, f AS STRING, h AS STRING, sw AS INTEGER)
   k = Decf(p)
   g$ = STRING$(k, "0")
   CnvSt g$
   Printf f, g$, h, sw
END SUB

SUB Printr (p AS INTEGER, q AS INTEGER, r AS INTEGER, f AS STRING, h AS STRING)
   lx = n(-2, i(r))
   Divd p, q, r: k = Decf(r)
   g$ = STRING$(k, "0"): CnvSt g$
   k = 1 - (LEN(h) = 0)
   IF Isf(p, 0) THEN
      Printf f, g$, "", k
   ELSE
      Printf f, g$, "", 0
      g$ = STRING$(lx, "0")
      Ratdec g$, p, q
      Printf ".", g$, h, k
   END IF
END SUB

SUB Prints (f AS STRING, sw AS INTEGER)
   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 IF
   END SELECT
END SUB

SUB Pwr10 (p AS INTEGER, k AS INTEGER)
   Letf p, 10
   Letf t3, CLNG(k)
   Letf t2, 0
   Modpwr p, t3, t2
END SUB

SUB Pwr2 (p AS INTEGER, k AS INTEGER)
   Letf p, 1
   Lsft p, ABS(k)
END SUB

SUB Ratdec (g AS STRING, p AS INTEGER, q AS INTEGER)
   c = LEN(g): lp = c \ LD: m = c - lp * LD
   IF m = 0 THEN m = LD: lp = lp - 1
   Letf t3, MY: c = 1
   FOR j = 0 TO lp
      c0 = 0
      FOR t = 0 TO 1
	 Mult p, t3, t1 '              dividend * MY
	 Divd p, q, t0 '               partial quotient t0
	 Divd t0, t3, t1 '             base MY digit
	 c0 = c0 * MY + n(0, i(t0))
      NEXT t
      s$ = LTRIM$(STR$(c0))
      IF j < lp THEN
	 c = c + LD
      ELSE
	 s$ = STRING$(LD - LEN(s$), "0") + s$
	 s$ = LEFT$(s$, m): c = c + m 'last digit
      END IF
      MID$(g, c - LEN(s$), LD) = s$ '  stuff
   NEXT j
END SUB

SUB Readst (p AS INTEGER, g AS STRING)
   g = LTRIM$(RTRIM$(g))
   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 = k - 1
   END IF
   lp = k \ LD: m = k - lp * LD
   IF m = 0 THEN m = LD: lp = lp - 1
   k = t + 1: t = m
   Letf t3, MD
   FOR j = 0 TO lp
      IF j THEN
	 Mult p, t3, t1: t = LD '      decimal base t3
      END IF
      Letf t1, VAL(MID$(g, k, t))
      Add p, t1: k = k + t '           convert to base MB
   NEXT j
   Lftj p, n(-2, i(p)) - 1
   IF sw THEN n(-1, i(p)) = -1
END SUB

SUB Rndf (p AS INTEGER, k AS INTEGER)
   k = ABS(k): f! = k * L10
   lp = INT(f!): m = INT((f! - lp) * LB)
   IF m THEN
      lp = lp + 1
   ELSE
      m = LB
   END IF
   IF lp >= vj THEN
      Errorh "overflow in Sub Rndf"
      Letf p, 1: ERROR 6: EXIT SUB
   END IF
   FOR j = 0 TO lp - 2
      n(j, i(p)) = INT(RND * MB)
   NEXT j
   m = m - 1: r = 1
   FOR j = 1 TO m: r = r * 2: NEXT
   m = r * 2: r = r + INT(RND * (m - r))
   n(-2, i(p)) = lp
   n(-1, i(p)) = 1
   n(lp - 1, i(p)) = r
   n(lp, i(p)) = 0
END SUB

SUB Rsft (p AS INTEGER, r AS INTEGER)
   k = r \ LB: m = r - k * LB
   t = n(-2, i(p)) - 1
   IF m > 0 THEN
      c = 0
      FOR j = t TO 0 STEP -1
	 c0 = c OR n(j, i(p)) '        paste carry
	 c = c0 AND M1: SftR c0, m '   save next, ShR element
	 n(j, i(p)) = c0 AND M1 '      mask off high bits
	 SftL c, LB
      NEXT j
      IF n(t, i(p)) = 0 AND t > 0 THEN t = t - 1
   END IF
   IF k > 0 THEN
      t = t - k
      IF t < 0 THEN Letf p, 0: EXIT SUB
      FOR j = 0 TO t
	 n(j, i(p)) = n(j + k, i(p)) ' ShR array
      NEXT j
   END IF
   n(-2, i(p)) = t + 1
   n(t + 1, i(p)) = 0
END SUB

SUB Sete (p AS INTEGER, j AS INTEGER, a AS INTEGER)
   n(j, i(p)) = ABS(a)
END SUB

SUB Setl (p AS INTEGER, a AS INTEGER)
   n(-2, i(p)) = ABS(a)
END SUB

SUB Sets (p AS INTEGER, a AS INTEGER)
   n(-1, i(p)) = SGN(a + .5)
END SUB

SUB Squ (p AS INTEGER, q AS INTEGER)
   l0 = n(-2, i(p))
   lx = l0 + l0 - 1
   IF lx >= vj THEN
      Errorh "overflow in Sub Squ"
      Letf q, 1: ERROR 6: EXIT SUB
   END IF
   IF l0 = 1 THEN
      c = CLNG(n(0, i(p))) * n(0, i(p))
      Letf q, c
   ELSE
      j = 0: n(-1, i(q)) = 1
      FOR k = 0 TO l0 - 1 '            initialize destination
	 c = CLNG(n(k, i(p))) * n(k, i(p))
	 n(j, i(q)) = c AND M1: SftR c, LB
	 n(j + 1, i(q)) = c: j = j + 2
      NEXT k
      '
      FOR k = 1 TO l0 - 1
	 c = 0: c2 = n(k, i(p)) * 2& ' add mixed terms
	 FOR j = 0 TO k - 1
	    c = c + c2 * n(j, i(p)) + n(j + k, i(q))
	    n(j + k, i(q)) = c AND M1: SftR c, LB
	 NEXT j
	 m = k + k
	 FOR j = m TO m + 1
	    c = c + n(j, i(q))
	    n(j, i(q)) = c AND M1: SftR c, LB
	 NEXT j
      NEXT k
      Lftj q, lx
   END IF
END SUB

SUB Subt (p AS INTEGER, q AS INTEGER)
   IF n(-1, i(p)) = n(-1, i(q)) THEN ' subtract
      ix = i(p): im = i(q)
      s = n(-2, ix) - n(-2, im)
      IF s = 0 THEN '                  equal lengths
	 FOR j = n(-2, ix) - 1 TO 0 STEP -1
	    s = n(j, ix) - n(j, im)
	    IF s THEN EXIT FOR
	 NEXT j
	 IF s = 0 THEN Letf p, 0: EXIT SUB
      END IF
      IF s < 0 THEN '                    p:= -(q - p)
	 SWAP ix, im
	 n(-1, i(p)) = -n(-1, i(p))
      END IF
      lx = n(-2, ix): lm = n(-2, im)
      FOR j = lm + 1 TO lx: n(j, im) = 0: NEXT
      c = MB
      FOR j = 0 TO lx
	 c = c + n(j, ix) + M2 - n(j, im)
	 n(j, i(p)) = c AND M1: SftR c, LB
      NEXT j
   ELSE '                                  add
      lx = n(-2, i(p)): lm = n(-2, i(q)): im = i(q)
      IF lx < lm THEN SWAP lx, lm: im = i(p)
      FOR j = lm + 1 TO lx: n(j, im) = 0: NEXT
      c = 0
      FOR j = 0 TO lx
	 c = c + n(j, i(p)) + n(j, i(q))
	 n(j, i(p)) = c AND M1: SftR c, LB
      NEXT j
   END IF
   Lftj p, lx
END SUB

SUB Swp (p AS INTEGER, q AS INTEGER)
   SWAP i(p), i(q)
END SUB

SUB Term
   CLOSE 'all files, then terminate
END SUB