'Subject: Williams's double-step, Lucas series P + 1 method finds factor
'         p of N if all prime factors but one of p+/-1 are at most B,
'         and a single larger factor is at most B1.
'Author : Sjoerd.J.Schaper
'Date   : 26-05-2004
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DEFLNG B, F
DECLARE SUB Lucadd (a, c)
'Calculate Lucas function V_c(P,1) in a
CONST Fx = 2147483629 '                fixed limit
CONST N = 0, d = 1, p = 2, dc = 3
CONST x0 = 4, x = 5, y = 6, z = 7, qfl = 8
CONST dq = 76, qcl = qfl + dq
t = LargeInit(qcl, "williap1")
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:
PRINT " k =";
FOR k = 1 TO 4
   PRINT k; : t = 0
   DO '                                initialize Lucas sequence V_1(P,1)
      q = 3 + INT(RND * 97)
      Letf x, CLNG(q): Gcd x, N, d
      IF NOT Isf(d, 1) THEN GOTO jump
      Squ x, dc: Inc dc, -4
      Moddiv dc, N '                   P^ 2 - 4 (mod N)
   LOOP WHILE Isf(dc, 0) '             discriminant <> 0
   '
   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 index 10! * 12
      END IF
      Lucadd x, c
      IF Isf(x, 2) THEN GOTO Step1
      IF r = Kr THEN
         Inc x, -2
         IF Kr < 16 THEN Kr = Kr * 2
         Gcd x, N, d
         IF NOT Isf(d, 1) THEN GOTO jump
         Inc x, 2: r = 0
      END IF
      fr = Nxtprm(sw): r = r + 1
   LOOP UNTIL fr > B
   '
   Squ x, i(0): Inc i(0), -2
   Moddiv i(0), N
   FOR q = 1 TO dq '                   store even indices V_2,.., _154
      Mult i(q - 1), i(0), i(q)
      IF q > 1 THEN
         Subt i(q), i(q - 2)
      ELSE
         Inc i(q), -2
      END IF
      Moddiv i(q), N
   NEXT q
   '
   Dup x, x0
   Letf p, 1: t = 2 '                  make search for one last big factor
   DO
      Dup x0, x
      f0 = fr: Lucadd x, (fr) '        giant step
      Dup x, d: Inc d, -2
      DO
         Sets d, 1: Modmlt p, d, N
         fr = Nxtprm(1)
         q = (fr - f0) \ 2 - 1
         IF q > dq THEN EXIT DO '      baby step:
         Dup x, d: Subt d, i(q) '      sqrt(V_m+n - 2)(V_m-n - 2)
      LOOP
      Gcd p, N, d
      IF NOT Isf(d, 1) THEN GOTO jump
   LOOP UNTIL fr > B1
NEXT k

jump:
PRINT
IF t > 0 THEN
   g$ = "(D/f)=" + STR$(Kronec(dc, d)) ' order
   g$ = g$ + "  t =" + STR$(t) '        phase,
   Prints g$ + "  fr =" + STR$(fr), 1 ' final index
END IF
Divd N, d, p '                         gcd in d
IF Isf(d, 1) OR Isf(p, 1) THEN
   Prints "fail", 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

SUB Lucadd (a, c)
   DO UNTIL c AND 1
      c = c \ 2
      Squ a, t0: Swp a, t0
      Inc a, -2: Divd a, N, t0 '       V_2n = (V_n)^ 2 - 2
   LOOP '                              c = (2^ i) * odd(k)
   Dup a, y: Dup a, z
   DO UNTIL c = 1 '                   R=>L Montgomery ladder
      c = c \ 2
      Squ y, t0: Swp y, t0
      Inc y, -2: Divd y, N, t0 '       double...
      IF c AND 1 THEN
         Lmul a, y: Subt a, z '        V_m+n = V_m * V_n - V_m-n
         Divd a, N, t0
      ELSE '                           ...and add
         Lmul z, y: Subt z, a
         Divd z, N, t0
      END IF
   LOOP
END SUB

