'Subject: Pollard's double-step P-1 method finds factor p of N if
'         p-1 is a product of small primes and a single larger prime.
'Author : Sjoerd.J.Schaper
'Date   : 24-05-2004
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DEFLNG B, F
CONST Fx = 2147483629 '                fixed limit
CONST N = 0, d = 1, a = 2, p = 3, qfl = 4
CONST dq = 76, qcl = qfl + dq
t = LargeInit(qcl, "pollarp1")
DIM LnB AS SINGLE, i(dq) '             difference table
FOR q = qfl TO qcl: i(q - qfl) = q: NEXT
CLS : tim# = TIMER: lin = 8
RANDOMIZE TIMER

begin:
LOCATE 1, 1
PRINT "   +---------------------------+"
PRINT "   |  input                    |"
PRINT "   |  N to split and bound B,  |"
PRINT "   |  0 to exit.               |"
PRINT "   |                           |"
PRINT "   +---------------------------+"
LOCATE 5, 7: INPUT ; "N"; Ni$
IF INSTR(Ni$, "'") THEN GOTO begin
IF VAL(Ni$) < 2 THEN GOTO eind
INPUT ; " B"; Bi$: B = 0 '             smoothness bound
IF ABS(VAL(Bi$)) < Fx THEN B = ABS(VAL(Bi$))
IF B = 0 THEN B = Fx \ 1000 '          default B
B1 = B * 100 - dq
IF B1 > Fx THEN B1 = Fx '              second bound
LnB = LOG(B)

CLS : LOCATE 8, 1
Readst N, Ni$: t = 0
Printn N, "  N = ", "", 2

IF IsPPrm(N, d, 5) THEN '              preliminary tests
   Prints " N is probably prime", 2
   GOTO begin
ELSE
   Swp p, d: Gcd N, p, d
   IF NOT Isf(d, 1) THEN GOTO jump
   IF IsSqr(N, d) THEN GOTO jump
END IF

Step1:
q = 2 + INT(RND * 98): t = 0
Letf a, CLNG(q): Gcd a, N, d '         random base a
IF NOT Isf(d, 1) THEN GOTO jump
sw = 0: fr = 1
r = 0: Kr = 4: t = 1
DO
   IF sw THEN
      c = fr ^ INT(LnB / LOG(fr)) '    prime-power
   ELSE
      c = 43545600 '                   initial exponent 10! * 12
   END IF
   Letf d, c
   Modpwr a, d, N '                    cumulative a-power
   IF Isf(a, 1) THEN GOTO Step1
   IF r = Kr THEN '                    periodic gcd check
      Inc a, -1
      IF Kr < 16 THEN Kr = Kr * 2
      Gcd a, N, d
      IF NOT Isf(d, 1) THEN GOTO jump
      Inc a, 1: r = 0
   END IF
   fr = Nxtprm(sw): r = r + 1
LOOP UNTIL fr > B

Squ a, i(0) '                          standard continuation:
Divd i(0), N, p
FOR q = 1 TO dq '                      store even powers a^ 2,... 154
   Mult i(q - 1), i(0), i(q)
   Divd i(q), N, p
NEXT q

Letf d, fr '                           make search for one last big factor
Modpwr a, d, N
Letf p, 1: t = 2
DO
   Inc a, -1
   Modmlt p, a, N: Inc a, 1 '          cumulative product
   IF r = dq THEN
      Gcd p, N, d
      IF NOT Isf(d, 1) THEN EXIT DO
      r = 0
   END IF
   f0 = fr: fr = Nxtprm(1)
   q = (fr - f0) \ 2 - 1
   Modmlt a, i(q), N
   r = r + 1
LOOP UNTIL fr > B1

jump:
IF t > 0 THEN  '                       print phase, final exponent
   Prints "  t =" + STR$(t) + "  fr =" + STR$(fr), 1
END IF
Divd N, d, p '                         gcd in d
IF Isf(d, 1) OR Isf(p, 1) THEN
   Prints "failure", 1
ELSE
   Printn d, "  f = ", "", 2
   Printn p, "N/f = ", "", 2
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

