'Subject: Lenstra's elliptic curve method finds factor p of N if
'         some random elliptic curve group has a smooth order mod p.
'         Uses homogeneous Montgomery form BY^2Z = X^3 + AX^2Z + XZ^2.
'         Recursive double-step version.
'Author : Sjoerd.J.Schaper
'Date   : 03-12-2005
'Code   : QuickBasic, PDS, VBdos
REM $LANG: "qb"
REM $INCLUDE: 'largeint.bi'
DEFLNG B, F
DECLARE SUB Factor (d, N)
DECLARE FUNCTION Splitf (d, N)
DECLARE SUB Eadd (xm, zm, xn, zn, xd, zd)
'Add points P_m and P_n, using P_(m-n) = P_d
DECLARE SUB ECM (d, N)
'Elliptic curve engine attempts to
'split N in N/ d and d with bound B
DECLARE SUB Edbl (xn, zn)
'Double point P_n
DECLARE SUB Ellmlt (xv, zv, c)
'Calculate point c*P in (xv: :zv)
DECLARE SUB Prntfr (k, fr)
'Print curve number and current factor

CONST Fx = 2147483629, lim = 131071 '  prime trial limit
CONST a = 2, p = 3, qu = 4, x0 = 5, z0 = 6
CONST xs = 7, zs = 8, xt = 9, zt = 10
CONST xu = 11, zu = 12, qfl = zu + 1 ' constant pointers
CONST dq = 153, qcl = qfl + dq, sfl = qcl
CONST stax = 6, scl = sfl + stax
t = LargeInit(scl, "ellcrvfr")
DIM SHARED B, N, i(dq), St(stax) '     difference table, stack
FOR q = qfl TO qcl: i(q - qfl) = q: NEXT
FOR t = sfl TO scl: St(t - sfl) = t: NEXT
CLS : tim# = TIMER: lin = 8
RANDOMIZE TIMER
N = 0: d = 1 '                         variable pointers

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
Readst N, Ni$
'
LOCATE 5, 7: PRINT SPACE$(25) + "|"
LOCATE 5, 7: INPUT "B"; Ni$: B = 0
IF ABS(VAL(Ni$)) < Fx \ 1000 THEN
   B = ABS(VAL(Ni$)) '                 smoothness bound
END IF
IF B < 10 THEN B = 1000 '              default B

LOCATE lin, 1
Printn N, " N = ", "", 2
Factor d, N
Prints "", 1
lin = CSRLIN
GOTO begin

eind:
LOCATE lin, 1
Ni$ = STR$(CSNG(TIMER - tim#))
Prints "Timer:" + Ni$ + " s", 1
CALL Term
SYSTEM

SUB Eadd (xm, zm, xn, zn, xd, zd)
   Mult xm, xn, t0
   Mult zm, zn, t1: Subt t0, t1
   Divd t0, N, t1: Modsqu t0, N '      t0:= (XmXn - ZmZn)^ 2 (mod N)
   Mult xm, zn, t3
   Mult zm, xn, t1: Subt t3, t1
   Divd t3, N, t1: Modsqu t3, N '      t3:= (XmZn - ZmXn)^ 2 (mod N)
   Mult t0, zd, xm: Moddiv xm, N '     xm:= t0 * Z|m-n| (mod N)
   Mult t3, xd, zm: Moddiv zm, N '     zm:= t3 * X|m-n| (mod N)
END SUB

SUB ECM (d, N)
DIM LnB AS SINGLE, g AS STRING
DIM c, fr, t, x
   B1 = B * 64
   IF B1 > Fx THEN B1 = Fx '           second bound
   LnB = LOG(B)
   '
   IF IsSqr(N, d) THEN
      Dup d, N: EXIT SUB '             N is square
   END IF
   '
Stp1:
   Letf a, 2: k = 0
   DO
      k = k + 1: t = 0
      Sets a, 1
      Inc a, 1 + INT(RND * 5)
      Sets a, SGN(RND - .5)
      x = 2 + INT(RND * 10) '          draw elements A, X
      Letf xs, CLNG(x): Letf zs, 1 '   P0 = (X: :1)
      Gcd xs, N, d
      IF NOT Isf(d, 1) THEN GOTO jump
      '
      sw = 0: fr = 1
      r = 0: Kr = 4: t = 1
      DO
         IF sw THEN '                  prime-power
            c = fr ^ INT(LnB / LOG(fr))
         ELSE
            c = 43545600 '             initial factor 10! * 12
         END IF
         Ellmlt xs, zs, c
         IF Isf(zs, 0) THEN GOTO Stp1
         IF r = Kr THEN
            IF Kr < 32 THEN Kr = Kr * 2
            Gcd zs, N, d
            IF NOT Isf(d, 1) THEN GOTO jump
            Prntfr k, fr: r = 0
         END IF
         fr = Nxtprm(sw): r = r + 1
      LOOP UNTIL fr > B
      '
      Dup xs, i(0): Dup zs, i(1)
      Edbl i(0), i(1) '                2*P
      Dup i(0), i(2): Dup i(1), i(3)
      Edbl i(2), i(3) '                4*P
      FOR q = 4 TO dq - 1 STEP 2 '     store even points 2*P,.., 154*P
         Dup i(q - 2), i(q)
         Dup i(q - 1), i(q + 1)
         Eadd i(q), i(q + 1), i(0), i(1), i(q - 4), i(q - 3)
      NEXT q
      '
      Dup xs, x0: Dup zs, z0
      Letf p, 1: t = 2
      DO
         Dup x0, xs: Dup z0, zs '      make search for one last big factor
         f0 = fr: Ellmlt xs, zs, fr '  giant step
         Dup zs, d
         DO
            Modmlt p, d, N
            fr = Nxtprm(1): q = fr - f0 - 2
            IF q >= dq THEN EXIT DO
            Mult xs, i(q + 1), d
            Mult zs, i(q), t0: Subt d, t0
            Sets d, 0: Divd d, N, t0 ' baby step: |XmZn - ZmXn| mod N
         LOOP
         Gcd p, N, d
         IF NOT Isf(d, 1) THEN GOTO jump
         Prntfr k, fr
      LOOP UNTIL fr > B1
   LOOP
   '
jump:
   Divd N, d, p '                      gcd in d
   IF Isf(d, 1) OR Isf(p, 1) THEN
      IF Isf(p, 1) THEN Swp p, d
      Swp p, N: GOTO Stp1
   END IF
   '
   IF t > 0 THEN
      g = "  x =" + STR$(x) '          log curve elements,
      c = SGN(Gets(a) OR 1) * Getw(a, 0)
      Prints g + "  a =" + STR$(c), 0
      c = x * (x * (x + c) + 1)
      Prints "  byż=" + STR$(c), 1
      g = "  k =" + STR$(k) '          curve number, final factor
      Prints g + "  fr =" + STR$(fr), 1
   END IF
   Swp p, N
END SUB

SUB Edbl (xn, zn)
   Mult xn, zn, t0: Divd t0, N, t1 '   t0:= Xn * Zn (mod N)
   Modsqu xn, N: Modsqu zn, N '        xn:= Xn^ 2 (mod N), zn:= Zn^ 2 (mod N)
   Mult a, t0, t3
   Add t3, xn: Add t3, zn '            t3:= Xn^ 2 + A * XnZn + Zn^ 2
   Subt xn, zn: Modsqu xn, N '         xn:= (Xn^ 2 - Zn^ 2)^ 2 (mod N)
   Lsft t0, 2: Mult t0, t3, zn '       zn:= 4 * XnZn * t3 (mod N)
   Moddiv zn, N
END SUB

SUB Ellmlt (xv, zv, c)
   c0 = ABS(c)
   DO UNTIL c0 AND 1
      Edbl xv, zv: c0 = c0 \ 2
   LOOP '                              c0 = (2^ i) * odd(k)
   Dup xv, xt: Dup zv, zt
   Dup xv, xu: Dup zv, zu '            V = T = -U = (2^ i) * P
   DO UNTIL c0 = 1 '                   R=>L Montgomery ladder
      Edbl xt, zt: c0 = c0 \ 2 '       T:= 2T
      IF c0 AND 1 THEN
         Eadd xv, zv, xt, zt, xu, zu ' V:= V + T, using V - T = U
      ELSE
         Eadd xu, zu, xt, zt, xv, zv ' U:= U - T, using U + T = V
      END IF
   LOOP
END SUB

SUB Factor (d, N)
   pr = 0: sq = 0
   IF Getl(N) > 2 THEN '               preliminary tests
      sq = IsSqr(N, d)
      IF sq THEN Swp N, d '            N is square
      pr = IsPPrm(N, d, 5) '           N or sqrt(N) is pprime
   END IF

   IF pr THEN
      Setl qu, 0
      EnQ qu, N, 1
   ELSE
      Triald qu, N, lim '              trial divide N
      IF NOT Isf(N, 1) THEN
         IF NOT Splitf(d, N) THEN EXIT SUB
      END IF
   END IF

   k = 1: t = 0
   WHILE ExQ(N, k, qu, t) '            read factors
      g$ = "": IF sq THEN k = k * 2
      IF k > 1 THEN g$ = "^ " + LTRIM$(STR$(k))
      Printn N, SPACE$(5), g$, 1
   WEND
END SUB

SUB Prntfr (k, fr)
STATIC Ton AS SINGLE
   IF TIMER - Ton > 1 THEN
      lin = CSRLIN: Ton = TIMER
      PRINT "  k ="; k; "  fr ="; fr; "     ";
      LOCATE lin, 1
   END IF
END SUB

FUNCTION Splitf (d, N)
Splitf = 0
   t = 0
   DO
      IF IsPPrm(N, d, -5) THEN '           N is probably prime
         k = 1
         FOR j = t TO 1 STEP -1
            DO
               Dup St(j), t0 '             trace possible common
               Divd t0, N, t1 '            factors in the stack
               IF NOT Isf(t0, 0) THEN EXIT DO
               Swp St(j), t1: k = k + 1
            LOOP
            IF Isf(t0, 1) AND Isf(t1, 0) THEN
               Swp St(t), St(j): t = t - 1 ' pop empty entry
            END IF
         NEXT j
         EnQ qu, N, k
         Swp St(t), N: t = t - 1 '         pop pprime
         '
      ELSE
         ECM d, N: t = t + 1 '             make search for divisor
         Printn N, "N/f = ", "", 2
         IF t > stax THEN
            Prints "  out of stack space", 1
            EXIT FUNCTION
         END IF
         Swp N, St(t): Swp d, N '          push N/ d, check d
      END IF
   LOOP UNTIL t < 0
Splitf = -1
END FUNCTION

