NB. Elliptic curve arithmetic and factorization.
NB. factor_ecj.ijs
NB. Cliff Reiter
NB. June 2003
NB. Corrections made February 2005
NB. (thanks to Roger Hui)

NB. Elliptic curves are E=.a,b where (y^2)=(x^3)+(a*x)+b
NB. Moduli are n where n-:0 corresponds to rational arithmetic

NB. ****************************************************************
NB. Main functions:
NB. fac_smpr n NB. get small prime factors in n 
NB. fac_p_1 n  NB. get a factor via Pollard p-1
NB. fac_rho n  NB. get a factor via Pollard rho
NB. fac_ecmj n NB. get facor via random elliptic curve
NB. primeq n   NB. use msr to determine whether n is prime
NB. factor n   NB. use the above to find the full factorization of n
NB. ****************************************************************


NB. *******************************
NB. elliptic curve arithmetic in
NB. modified projective coordinates
NB. *******************************

NB. Elliptic curve double point (mod proj coord)
NB. ab ecdj n Q
ecdj=: 2 : 0
'a b'=.m.
'x y z'=.3{.y.,1
if. 0 e. y,z do. 'x2 y2 z2'=.0 1 0x else.
  m=.(3**:x)+a**:*:z
  s=.4*x*yy=.*:y
  x2=.(*:m)+_2*s
  y2=.(m*s-x2)+_8**:yy
  z2=.2*y*z
  end.
n.|x2,y2,z2
)

NB. Ellitpic curve add (mod proj coord)
NB. P ab ecaj n Q
ecaj=: 2 : 0
:
'a b'=.m.
'x1 y1 z1'=.3{.x.,1
'x2 y2 z2'=.3{.y.,1
if. z1=0 do. 'x3 y3 z3'=.x2,y2,z2 else.
  if. z2=0 do. 'x3 y3 z3'=.x1,y1,z1 else.
  NB. genl case
  u1=.x2*z12=.*:z1
  u2=.x1*z22=.*:z2
  s1=.y2*z1*z12
  s2=.y1*z2*z22
  w=.u1-u2
  r=.s1-s2
  if. w=0 do. 'x3 y3 z3'=.m. ecdj n. x1,y1,z1 else.
    NB. continue genl
    t=.u1+u2
    m=.s1+s2
    x3=.(*:r)-t*w2=.*:w
    y3=.-:(r*(_2*x3)+t*w2)-m*w*w2
    z3=.z1*z2*w
    end. end. end.
n.|x3,y3,z3
)

NB. Scalar mult ladder (mod proj coord)
NB. m ab ecmj n P
ecmj=: 2 : 0
:
if. x.=0 do. q=.0 1 0x else.
  q=.y.=.3{.y.,1
  for_mn. #.}.}:|:#:3 1*x. do.
    q=.m. ecdj n. q
    select. mn
    case. 2 do. q=.q m. ecaj n. y.
    case. 1 do. q=.q m. ecaj n. 1 _1 1* y.
    end.
  end.
end.
q
)


NB. obtain a random EC mod y. (smaller than 2^31)
ranEC=:3 : 0
whilst. g=y. do.
  'x y a'=.x:?3#y.<.<:2^31
  b=.y.|(*:y)-x*(a+*:x)
  g=.y. +. +/ 4 27* (a,b) ^ 3 2
  end.
g,a,b,x,y
)

NB. ***********************************************
NB. Functions for main elliptic curve factorization
NB. ***********************************************

NB. [list of B1 B2s] fac_ecmj n
fac_ecmj=: 3 : 0"_ 0
(20 d_B1_B2 y.) fac_ecmj y.
:
max=.#x.
for_ecc. i.max do.
  'B1 B2'=.ecc{x.
  Last_random_ec=:ranEC y.
  'g a b x y'=.Last_random_ec
  if. -.g e. 1,y. do. g,0,ecc return. end.
  Q=. x,y,1
  ab=.a,b
  Q=.B1 ab fac_ecmj_s1 y. Q
  g=.y. +. {:Q
  if. -. g e. 1,y. do. g,1,ecc return. end.
  if. g=1 do.
    Q=.(B1,B2) ab fac_ecmj_s2 y. Q
    g=.y. +. {:Q
    if. -. g e. 1,y. do. g,2,ecc return. end.
    end.
end.
(g,_1,max)
)

NB. Gives a range of B1 B2 values
NB. [num_steps] d_B1_B2 m
NB. log to "usual" B1
NB. Handy if relatively small factors appear
d_B1_B2=:3 : 0
20 d_B1_B2 y.
:
lg=.^. y.
rt=.^(%:%2)+%:(*^.)^. %:y.
<. (,.(*^.)) lg*(rt%lg)^ (i.%<:) x.
)

NB. B1 (ab) fac_ecmj_s1 (n) xy1
fac_ecmj_s1=:2 : 0
:
Q=.y.
for_p. p: i. p:^:_1 x. do.
  Q=.(p^<.p ^. x.) m. ecmj n. Q
  end.
Q
)

NB. (B1,B2) ab fac_ecmj_s2 n Q
fac_ecmj_s2=:2 : 0
:
'B1 B2'=.x.
Q=.y.
Sd=.,:S1=.m. ecdj n. Q
D=.100
for_d. 1+i.D-1 do.
  Sd=.Sd,S1 m. ecaj n. {:Sd
  end.
'piB1 piB2'=.  p:^:_1 B1,B2
ps=.p: piB1+i.piB2-piB1
Q=.({.ps) m. ecmj n. Q
for_di. <:-:(}.-}:)ps do.
  for. i.<.(di%D) do. Q=.Q m. ecaj n. {:Sd end.
  Q=.Q m. ecaj n. (D|di){Sd
  end.
Q
)

NB. ***********************************************
NB. Functions for auxiliary factorization tecniques
NB. ***********************************************


NB. [B] fac_p_1 n
fac_p_1=: 3 : 0"0
(d_p_1_B y.) fac_p_1 y.
:
c=.2x
exp=.y.&|@^
for_k. i. p:^:_1 x. do.
  p=.p: k
  c=.c exp (p^<.p ^. x.)
  g=.y. +. c-1
  if. -. g e. 1,y. do. g,1 return. end.
end.
g,_1
)

d_p_1_B=: 10000"_ <. >.@:(^&0.5 )

NB. Pollard rho
NB. [maxiter] fac_rho n
fac_rho=:3 : 0"0
(d_rho_maxj y.) fac_rho y.
:
n=.y.
f=.1x&+@*:
x1=.f 2x
x2=.f f 2x
g=.n +. n|x2-x1
j=.1
while. (g e. 1,n)*. j<x. do.
  x1=.n|f x1
  x2=.n|f f x2
  g=.n +. n|x2-x1
  j=.j+1
end.
g,j
)

d_rho_maxj=: 10000&<.@<.@(^&(%3))

NB. Essentially trial division by small primes
NB. [num_primes] fac_rho n
fac_smpr=: 3 : 0
24 fac_smpr y.
:
p=.p:i.x.
z=.0
whilst. 0<+/peg do.
  peg=.p=p +. y.  NB. prime equals gcd
  z=.z+peg
  y.=.y.%*/peg#p
end.
(z#p);y.
)

smpr=:>@{.@fac_smpr

rmsmpr=: >@{:@fac_smpr

randi=: (#.?)@:(#&10x)"0

NB. ***********************************************
NB. Functions for probablistic testing of primality
NB. ***********************************************

NB. for prob prime testing
pwr2in=:3 :0
j=.1
while. 0=(2x ^j)|y. do.
  j=.j+1
end.
j-1
)

msr=:4 : 0"0"1 0
b=.x.
n=.x: y.
k=.pwr2in n-1
m=.(n-1)%2x ^k
x=.n&|@(b&^)m
if. 1=x do. 1 else.
j=.0
while. (x~:n-1)*.j<k do.
  x=.n|x^2
  j=.j+1
end.
if.x = n-1 do.1 else.0 end.
end.
)

primeq=:3 : 0"0
100 primeq y.
:
b=.p:i.x.>.3
if. y. e. b do. 1 return. end.
if. -.*/(2{.b) msr y. do. 0 return. end.
*/(2}.b) msr y.
)

NB. ****************************
NB. Main functions for factoring
NB. ****************************

factor=: 3 : 0"0
'ps fs'=.fac_smpr y.
if. fs=1 do. fs=.i.0 end.
q=.primeq fs
ps=.ps,q#fs
fs=.(-.q)#fs
while. 0<#fs do.
  f0=.{.fs
  'g v'=.fac_p_1 f0
  if. g e. 1,f0 do.
    'g v'=.fac_rho f0
    while. g e. 1,f0 do.
      'g v'=.2{.>{.fac_ecmj f0
      end.
    end.
  fs=.(}.fs),g,f0%g
  q=.primeq fs
  ps=.ps,q#fs
  fs=.(-.q)#fs
  end.
/:~ps
)

NB. Timing utilities
ts=:6!:2 , 7!:2@]

time=:6!:2

fac_time=:3 : 0"0
t=.time 'z=:factor y.'
t;z
)

NB. Test numbers and trials
test_numbers_and_trials=:0 : 0
]i1=:1001
]i2=:10001
]i3=:1000001

]p1=:*/x: p: 1000 * 1 10
]p2=:*/x: p: 1000000 * 1 10
]p3=:*/x: p: 1000000 * 1 2 3 4 5

]m29=:_1+2^29x
]m57=:_1+2^57x
]m71=:_1+2^71x
]m83=:_1+2^83x
]m97=:_1+2^97x
]m101=:_1+2^101x
]m61m89=:*/x: _1+2^x: 61 89

mm=:152415787533657061564561727x

NB. For Vector Notes 
(9!:1)1710440016
,.x=:rmsmpr@randi 10#25
time 'x=:rmsmpr@randi 10#25'
fac_p_1 x
time 'fac_p_1 x'
fac_rho x
time 'fac_rho x'

(9!:1)1710440016
,.x=:rmsmpr@randi 10#25
fac_ecmj x
time 'fac_ecmj x'
)

NB.  Sample usage:
NB.  factor m97