Elliptic curve computations mod N in projective coordinates

On Jul 10, 2007 9:08 AM Alasdair McAndrew? wrote:

I have written some highly unoptimized code for factoring integers using Lenstra's elliptic curve method, with the "birthday paradox" phase two developed by Richard Brent. Even at this stage, it can factor the seventh Fermat number 2^2^7+1 in 352 seconds, as opposed to 1877 seconds by the in-built factoring method. If anybody is interesting in developing this code further, do let me know.

All I really want to do is to extend the methods in intfact.spad to include ecm - which the file itself recommends.

-Alasdair

reference

(1) -> mod(a:Integer,b:PositiveInteger):Integer == (a>0 => a rem b::Integer; b + (a rem b))

   Function declaration mod : (Integer, PositiveInteger) -> Integer has
      been added to workspace.

a is a parameter which describes the elliptic curve, and which is the one used for addition and point doubling.

ecAdd(a:Integer,X:List Integer,Y:List Integer,N:Integer):List Integer ==
  local u,v,A,x3,y3,z3
  if (X=Y)::Boolean then return ecDouble(a,X,N)
  else
    u := Y.2*X.3-X.2*Y.3
    v := Y.1*X.3-X.1*Y.3
    A := u^2*X.3*Y.3-v^3-2*v^2*X.1*Y.3
    x3 := (v*A) rem N
    y3 := (u*(v^2*X.1*Y.3-A)-v^3*X.2*Y.3) rem N
    z3 := (v^3*X.3*Y.3) rem N
    return [x3,y3,z3]

   Function declaration ecAdd : (Integer, List(Integer), List(Integer)
      , Integer) -> List(Integer) has been added to workspace.

ecDouble(a:Integer,X:List Integer,N:Integer):List Integer ==
  local w,s,B,h,x3,y3,z3
  w := a*X.3^2+3*X.1^2
  s := X.2*X.3
  B := X.1*X.2*s
  h := w^2 - 8*B
  x3 := (2*h*s) rem N
  y3 := (w*(4*B-h)-8*X.2^2*s^2) rem N
  z3 := (8*s^3) rem N
  return [x3,y3,z3]

   Function declaration ecDouble : (Integer, List(Integer), Integer)
       -> List(Integer) has been added to workspace.

ecMultiply(a:Integer,X:List Integer,n:Integer,N:Integer):List Integer ==
  local lst,P
  lst := wholeRagits(n::RadixExpansion 2)
  P:=X
  for i in rest lst repeat
    P:=ecDouble(a,P,N)
    if i = 1 then P:=ecAdd(a,P,X,N)
  return P

   Function declaration ecMultiply : (Integer, List(Integer), Integer, 
      Integer) -> List(Integer) has been added to workspace.

Here is the main routine; you have to enter both the bounds B and r; for example:

   ecFactor(2^101-1,10000,700)

ecFactor(N:Integer,B:Integer,r:Integer):List Integer ==
  local a:Integer, x:Integer,y:Integer,b,g,curve,found,Q,f,i,tmp
  found:=false
  i:=0
  while not(found) repeat
    curve:=false
    while not(curve) repeat
      a:=random(N)
      x:=random(N)
      y:=random(N)
      b:=y^2-x^3-a*x
      g:=gcd(4*a^3+27*b^2,N)
      if g = 1 then 
        curve:=true
        i:=i+1
      if (i rem 20)=0 then output i
    Q:=ecPhaseOne(a,[x,y,1],B,N)
    f := gcd(Q.3,N)
    if f>1 and f<N then
      found:=true
    else
      tmp:=ecPhaseTwo(a,Q,N,r)
      if tmp>1 and tmp<N then
        found:=true
  return [f,N/f]

   Function declaration ecFactor : (Integer, Integer, Integer) -> List(
      Integer) has been added to workspace.

ecPhaseOne(a:Integer,X:List Integer,B:PositiveInteger,N:PositiveInteger):List Integer ==
  local P,i,e
  P := X
  for i in primes(2,B) repeat
    e := wholePart(log2(B::DoubleFloat)/log2(i::DoubleFloat))::PositiveInteger
    P := ecMultiply(a,P,i^e,N)
  return P

   Function declaration ecPhaseOne : (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> List(Integer) has been added
      to workspace.

ecPhaseTwo(a:Integer,Q:List Integer,N:PositiveInteger,r:PositiveInteger):Integer ==
  local randr,Qs,tmp,i,j,k,d,result
  result := 1
  randr := rest wholeRagits(2^r+random(2^r)::RadixExpansion 2)
  Qs := [Q]
  for i in 2..r repeat
    tmp := ecDouble(a,Qs.(i-1),N)
    if randr.i = 1 then tmp := ecAdd(a,tmp,Q,N)
    Qs := append(Qs,[tmp])
  for i in 2..r repeat
    if member?(Qs.i,[Qs.k for k in 1..i-1]) then
      --j := position(Qs.i,[Qs.k for k in 1..i-1])
      d := reduce(*,[reduce(*,[Qs.i.3-Qs.j.3 for j in i+1..r]) for i in _
  1..r-1])
      if gcd(N,d)>1 then result:=gcd(N,d)
  return result

   Function declaration ecPhaseTwo : (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> Integer has been added to 
      workspace.

See http://www.alpertron.com.ar/ECM.HTM for values of B.

For example:

)set message time on

factor 2^101-1

(1)

Time: 1.32 (EV) + 0.03 (OT) = 1.35 sec

ecFactor(2^101-1,10000,700)

Compiling function ecDouble with type (Integer, List(Integer), 
      Integer) -> List(Integer)

Compiling function ecAdd with type (Integer, List(Integer), List(
      Integer), Integer) -> List(Integer)

Compiling function ecMultiply with type (Integer, List(Integer), 
      Integer, Integer) -> List(Integer)

Compiling function ecPhaseOne with type (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> List(Integer)

Compiling function ecPhaseTwo with type (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> Integer

Compiling function ecFactor with type (Integer, Integer, Integer)
       -> List(Integer)

Compiling function G27 with type Integer -> Boolean

Compiling function G29 with type NonNegativeInteger -> Boolean 
   20

(2)

Time: 0.02 (IN) + 3.79 (EV) + 0.05 (OT) + 0.01 (GC) = 3.88 sec

test( reduce(*,%)=2^101-1 )

(3)

Time: 0 sec

)abbrev package EC EllipticCurve
--#pile
--#include "axiom"
--import from Integer
EllipticCurve(): with
    ecAdd:(Integer,List Integer,List Integer,Integer)->List Integer
    ecDouble:(Integer,List Integer,Integer) -> List Integer
    ecMultiply:(Integer,List Integer,Integer,Integer) -> List Integer
  == add
    ecAdd(a:Integer,X:List Integer,Y:List Integer,N:Integer):List Integer ==
      if X=Y then
        ecDouble(a,X,N)
      else
        u := Y.2*X.3-X.2*Y.3
        v := Y.1*X.3-X.1*Y.3
        A := u^2*X.3*Y.3-v^3-2*v^2*X.1*Y.3
        x3 := (v*A) rem N
        y3 := (u*(v^2*X.1*Y.3-A)-v^3*X.2*Y.3) rem N
        z3 := (v^3*X.3*Y.3) rem N
        [x3,y3,z3]

    ecDouble(a:Integer,X:List Integer,N:Integer):List Integer ==
      w := a*X.3^2+3*X.1^2
      s := X.2*X.3
      B := X.1*X.2*s
      h := w^2 - 8*B
      x3 := (2*h*s) rem N
      y3 := (w*(4*B-h)-8*X.2^2*s^2) rem N
      z3 := (8*s^3) rem N
      [x3,y3,z3]

    ecMultiply(a:Integer,X:List Integer,n:Integer,N:Integer):List Integer ==
      -- import from RadixExpansion 2, List Integer
      lst := wholeRagits(n::RadixExpansion 2)
      P:=X
      for i in rest lst repeat
        P:=ecDouble(a,P,N)
        if i = 1 then P:=ecAdd(a,P,X,N)
      P

   Compiling FriCAS source code from file 
      /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/5435117267253907166-25px007.spad
      using old system compiler.
   EC abbreviates package EllipticCurve 
------------------------------------------------------------------------
   initializing NRLIB EC for EllipticCurve 
   compiling into NRLIB EC 
   compiling exported ecAdd : (Integer,List Integer,List Integer,Integer) -> List Integer
Time: 0.34 SEC.

   compiling exported ecDouble : (Integer,List Integer,Integer) -> List Integer
Time: 0.07 SEC.

   compiling exported ecMultiply : (Integer,List Integer,Integer,Integer) -> List Integer
Time: 0 SEC.

(time taken in buildFunctor:  0)

;;;     ***       |EllipticCurve| REDEFINED

;;;     ***       |EllipticCurve| REDEFINED
Time: 0 SEC.

   Cumulative Statistics for Constructor EllipticCurve
      Time: 0.41 seconds

   finalizing NRLIB EC 
   Processing EllipticCurve for Browser database:
--->-->EllipticCurve(constructor): Not documented!!!!
--->-->EllipticCurve((ecAdd ((List (Integer)) (Integer) (List (Integer)) (List (Integer)) (Integer)))): Not documented!!!!
--->-->EllipticCurve((ecDouble ((List (Integer)) (Integer) (List (Integer)) (Integer)))): Not documented!!!!
--->-->EllipticCurve((ecMultiply ((List (Integer)) (Integer) (List (Integer)) (Integer) (Integer)))): Not documented!!!!
--->-->EllipticCurve(): Missing Description
; compiling file "/var/aw/var/LatexWiki/EC.NRLIB/EC.lsp" (written 22 NOV 2024 01:07:36 AM):

; wrote /var/aw/var/LatexWiki/EC.NRLIB/EC.fasl
; compilation finished in 0:00:00.028
------------------------------------------------------------------------
   EllipticCurve is now explicitly exposed in frame initial 
   EllipticCurve will be automatically loaded when needed from 
      /var/aw/var/LatexWiki/EC.NRLIB/EC

factor(N:Integer,B:Integer,r:Integer):List Integer ==
  local a:Integer, x:Integer,y:Integer,b,g,curve,found,Q,f,i,tmp
  found:=false
  i:=0
  while not(found) repeat
    curve:=false
    while not(curve) repeat
      a:=random(N)
      x:=random(N)
      y:=random(N)
      b:=y^2-x^3-a*x
      g:=gcd(4*a^3+27*b^2,N)
      if g = 1 then 
        curve:=true
        i:=i+1
      if (i rem 20)=0 then output i
    Q:=ecPhase1(a,[x,y,1],B,N)
    f := gcd(Q.3,N)
    if f>1 and f<N then
      found:=true
    else
      tmp:=ecPhase2(a,Q,N,r)
      if tmp>1 and tmp<N then
        found:=true
  return [f,N/f]

   Function declaration factor : (Integer, Integer, Integer) -> List(
      Integer) has been added to workspace.

Time: 0 sec

ecPhase1(a:Integer,X:List Integer,B:PositiveInteger,N:PositiveInteger):List Integer ==
  local P,i,e
  P := X
  for i in primes(2,B) repeat
    e := wholePart(log2(B::DoubleFloat)/log2(i::DoubleFloat))::PositiveInteger
    P := ecMultiply(a,P,i^e,N)$EC
  return P

   Function declaration ecPhase1 : (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> List(Integer) has been added
      to workspace.

Time: 0 sec

ecPhase2(a:Integer,Q:List Integer,N:PositiveInteger,r:PositiveInteger):Integer ==
  local randr,Qs,tmp,i,j,k,d,result
  result := 1
  randr := rest wholeRagits(2^r+random(2^r)::RadixExpansion 2)
  Qs := [Q]
  for i in 2..r repeat
    tmp := ecDouble(a,Qs.(i-1),N)$EC
    if randr.i = 1 then tmp := ecAdd(a,tmp,Q,N)$EC
    Qs := append(Qs,[tmp])
  for i in 2..r repeat
    if member?(Qs.i,[Qs.k for k in 1..i-1]) then
      --j := position(Qs.i,[Qs.k for k in 1..i-1])
      d := reduce(*,[reduce(*,[Qs.i.3-Qs.j.3 for j in i+1..r]) for i in _
  1..r-1])
      if gcd(N,d)>1 then result:=gcd(N,d)
  return result

   Function declaration ecPhase2 : (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> Integer has been added to 
      workspace.

Time: 0 sec

See http://www.alpertron.com.ar/ECM.HTM for values of B.

For example:

)set message time on

factor 2^101-1

(4)

Time: 1.33 (EV) = 1.34 sec

factor(2^101-1,10000,700)

Compiling function ecPhase1 with type (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> List(Integer)

Compiling function ecPhase2 with type (Integer, List(Integer), 
      PositiveInteger, PositiveInteger) -> Integer

Compiling function factor with type (Integer, Integer, Integer) -> 
      List(Integer) 
   20

(5)

Time: 0.03 (IN) + 4.17 (EV) + 0.03 (OT) + 0.02 (GC) = 4.25 sec

test( reduce(*,%)=2^101-1 )

(6)

Time: 0 sec