``````NB. Elliptic curve arithmetic and factorization.
NB. factor_ecj.ijs
NB. Cliff Reiter
NB. June 2003
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

``````