help links subscribe changes refresh edit

	SandBox
SandBoxElemntry ←	↑	→ SandBoxEtc

SandBoxEllipticCurves

Edit detail for SandBoxEllipticCurves revision 1 of 2

	1 2
Editor: Time: 2007/11/18 18:01:14 GMT-8
Note: small optimizations

changed:
-
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":http://grampus.jaist.ac.jp:8080/miyaji-lab/member/PaperPS/ICICS97r1.pdf

\begin{axiom}
mod(a:Integer,b:PositiveInteger):Integer == (a>0 => a rem b::Integer; b + (a rem b))
\end{axiom}

'a' is a parameter which describes the elliptic curve, and which is
the one used for addition and point doubling.
\begin{axiom}
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]

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]

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  
\end{axiom}

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

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

\begin{axiom}
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]

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

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
\end{axiom}

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

For example:
\begin{axiom}
)set message time on
factor 2^101-1
\end{axiom}
\begin{axiom}
ecFactor(2^101-1,10000,700)
\end{axiom}
\begin{axiom}
test( reduce(*,%)=2^101-1 )
\end{axiom}


From BillPage Tue Jul 10 22:45:51 -0500 2007
From: Bill Page
Date: Tue, 10 Jul 2007 22:45:51 -0500
Subject: spad
Message-ID: <20070710224551-0500@wiki.axiom-developer.org>

\begin{spad}
)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  
\end{spad}

\begin{axiom}
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]

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

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
\end{axiom}

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

For example:
\begin{axiom}
)set message time on
factor 2^101-1
\end{axiom}
\begin{axiom}
factor(2^101-1,10000,700)
\end{axiom}
\begin{axiom}
test( reduce(*,%)=2^101-1 )
\end{axiom}

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

axiom

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.

Type: Void

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

axiom

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.

Type: Void

axiom

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.

Type: Void

axiom

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.

Type: Void

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

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

axiom

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.

Type: Void

axiom

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.

Type: Void

axiom

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.

Type: Void

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

For example:

axiom

)set message time on

factor 2^101-1

$\label{eq1}{7432339208719}\ {341117531003194129}$

(1)

Type: Factored(Integer)

axiom

Time: 0.12 (IN) + 13.45 (EV) + 0.16 (OT) = 13.73 sec

axiom

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

axiom

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

axiom

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

axiom

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

axiom

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

axiom

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

axiom

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

axiom

Compiling function G731 with type Integer -> Boolean

axiom

Compiling function G733 with type NonNegativeInteger -> Boolean 
   20

$\label{eq2}\left[{7432339208719}, \:{341117531003194129}\right]$

(2)

Type: List(Integer)

axiom

Time: 0.26 (IN) + 28.01 (EV) + 0.41 (OT) = 28.68 sec

axiom

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

$\label{eq3} \mbox{\rm true}$

(3)

Type: Boolean

axiom

Time: 0.02 (OT) = 0.02 sec

spad --Bill Page, Tue, 10 Jul 2007 22:45:51 -0500 reply

spad

)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

spad

   Compiling FriCAS source code from file 
      /var/zope2/var/LatexWiki/4146989362989036246-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
****** comp fails at level 8 with expression: ******
error in function ecAdd 

(IF (= X Y)
    (|ecDouble| |a| X N)
    (SEQ
     (LET |u|
       (- (* (Y 2) (X 3)) (* (X 2) (Y 3))))
     (LET |v|
       (- (* (Y 1) (X 3)) (* (X 1) (Y 3))))
     (LET A
       (- (- (* (* | << | (** |u| 2) | >> | (X 3)) (Y 3)) (** |v| 3))
          (* (* (* 2 (** |v| 2)) (X 1)) (Y 3))))
     (LET |x3|
       (|rem| (* |v| A) N))
     (LET |y3|
       (|rem|
        (- (* |u| (- (* (* (** |v| 2) (X 1)) (Y 3)) A))
           (* (* (** |v| 3) (X 2)) (Y 3)))
        N))
     (LET |z3|
       (|rem| (* (* (** |v| 3) (X 3)) (Y 3)) N))
     (|exit| 1 (|construct| |x3| |y3| |z3|))))
****** level 8  ******
$x:= (** u 2)
$m:= $EmptyMode
$f:=
((((|v| #) (|u| #) (N # #) (Y # #) ...)))

   >> Apparent user error:
   NoValueMode
    is an unknown mode

axiom

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.

Type: Void

axiom

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.

Type: Void

axiom

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.

Type: Void

axiom

Time: 0 sec

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

For example:

axiom

)set message time on

factor 2^101-1

$\label{eq4}{7432339208719}\ {341117531003194129}$

(4)

Type: Factored(Integer)

axiom

Time: 14.42 (EV) + 0.01 (OT) = 14.43 sec

axiom

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

   You cannot now use EllipticCurve in the context you have it.

axiom

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

   There are 1 exposed and 3 unexposed library operations named reduce 
      having 2 argument(s) but none was determined to be applicable. 
      Use HyperDoc Browse, or issue
                             )display op reduce
      to learn more about the available operations. Perhaps 
      package-calling the operation or using coercions on the arguments
      will allow you to apply the operation.

   Cannot find a definition or applicable library operation named 
      reduce with argument type(s) 
                                 Variable(*)
                              Factored(Integer)

      Perhaps you should use "@" to indicate the required return type, 
      or "$" to specify which version of the function you need.