Modern Elliptic Curve Factorization, Part 1
April 23, 2010
(define (add P1 P2 P0 N)
(define (square n) (* n n))
(let* ((x0 (car P0)) (x1 (car P1)) (x2 (car P2))
(z0 (cdr P0)) (z1 (cdr P1)) (z2 (cdr P2))
(t1 (modulo (* (+ x1 z1) (- x2 z2)) n))
(t2 (modulo (* (- x1 z1) (+ x2 z2)) n)))
(cons (modulo (* (square (+ t2 t1)) z0) n)
(modulo (* (square (- t2 t1)) x0) n))))
(define (double P An Ad N)
(define (square n) (* n n))
(let* ((x (car P)) (z (cdr P))
(t1 (modulo (square (+ x z)) N))
(t2 (modulo (square (- x z)) N))
(t3 (- t1 t2)))
(cons (modulo (* t1 t2 4 Ad) N)
(modulo (* (+ (* t3 An) (* t2 Ad 4)) t3) N))))
(define (multiply K P An Ad N)
(cond ((zero? K) (cons 0 0)) ((= K 1) P) ((= K 2) (double P An Ad N))
(else (let loop ((ks (cdr (digits K 2))) (Q (double P An Ad N)) (R P))
(cond ((null? ks) R)
((odd? (car ks))
(loop (cdr ks) (double Q An Ad N) (add Q R P N)))
(else (loop (cdr ks) (add R Q P N) (double R An Ad N))))))))
(define (curve12 N s)
(let* ((u (modulo (- (* s s) 5) N))
(v (modulo (* 4 s) N)) (v-u (- v u)))
(values (modulo (* (* v-u v-u v-u) (+ u u u v)) N)
(modulo (* 4 u u u v) N)
(cons (modulo (* u u u) N)
(modulo (* v v v) N)))))
For example:
> (define n 5569379)
> (define s 4921134)
> (define An #f)
> (define Ad #f)
> (define P #f)
> (call-with-values
(lambda () (curve12 n s))
(lambda (a b c)
(set! An a)
(set! Ad b)
(set! P c)))
> An
3660795
> Ad
1539615
> P
(4087911 . 2770287)
> (define P2 (double P An Ad n))
> P2
(3539269 . 4477480)
> (define P3 (add P2 P P n))
> P3
(2517168 . 4993956)
> (define P4 (double P2 An Ad n))
> P4
(4683984 . 2280932)
> (define P6 (double P3 An Ad n))
> P6
(3440206 . 1480034)
> (define P7 (add P4 P3 P n))
> P7
(2544042 . 2445346)
> (define P13 (add P7 P6 P n))
> P13
(577485 . 4434222)
> (multiply 13 P An Ad n)
(577485 . 4434222)
We use digits from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/t9gygODg.
Programming Praxis 
[…] today’s Programming Praxis exercise we have to implement an algorithm for curve factorization. The Scheme […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/04/23/programming-praxis-modern-elliptic-curve-factorization-part-1/ for a version with comments):
add :: Integral a => (a, a) -> (a, a) -> (a, a) -> a -> (a, a) add (x1, z1) (x2, z2) (x0, z0) n = (f (+) z0, f (-) x0) where f g m = (g ((x1-z1)*(x2+z2)) ((x1+z1)*(x2-z2)))^2 * m `mod` n double :: Integral a => (a, a) -> a -> a -> a -> (a, a) double (x,z) an ad n = (mod x' n, mod z' n) where x' = 4 * ad * (x-z)^2 * (x+z)^2 z' = (4 * ad * (x-z)^2 + t * an) * t t = (x+z)^2 - (x-z)^2 multiply :: Integral a => Int -> (a, a) -> a -> a -> a -> (a, a) multiply 0 _ _ _ _ = (0,0) multiply 1 p _ _ _ = p multiply k p an ad n = let half = multiply (div k 2) p an ad n in if odd k then add half (multiply (div k 2 + 1) p an ad n) p n else double half an ad n curve12 :: Integral a => a -> a -> (a, a, (a, a)) curve12 n s = ((v-u)^3 * (3 * u + v) `mod` n, 4 * u^3 * v `mod` n, (u^3 `mod` n, v^3 `mod` n)) where u = s * s - 5 `mod` n v = 4 * s `mod` nPython version.
def init( n, seed ): u = ( seed*seed - 5 ) % n v = 4*seed % n num = pow(v-u, 3, n) * (3*u + v) % n den = 4 * pow(u,3,n) * v % n A = num, den P = (pow(u, 3, n), pow( v, 3, n) ) return A, P def add( p1, p2, p0, n ): x1,z1 = p1 x2,z2 = p2 x0,z0 = p0 t1 = (x1-z1)*(x2+z2) t2 = (x1+z1)*(x2-z2) newx = z0 * pow( t1 + t2, 2, n ) % n newz = x0 * pow( t1 - t2, 2, n ) % n return newx, newz def double( p, A, n ): x,z = p An, Ad = A t1 = pow(x + z, 2, n) t2 = pow(x - z, 2, n) t = t1 - t2 newx = t1 * t2 * 4 * Ad % n newz = (4 * Ad * t2 + t * An) * t % n return newx, newz def multiply( m, p, A, n ): if m == 0: return 0,0 elif m == 1: return p else: q = double( p, A, n ) if m == 2: return q b = 1 while b < m: b <<= 1 b >>= 2 r = p while b: if m&b: q, r = double(q, A, n), add(q, r, p, n) else: q, r = add(r, q, p, n), double(r, A, n) b >>= 1 return r n,s = 5569379,4921134 A, P1 = init(n, s) P2 = double(P1,A,n) P3 = add(P2,P1,P1,n) P4 = double(P2, A, n) P6 = double(P3, A, n) P7 = add(P4,P3,P1,n) P13 = add(P6,P7,P1,n) P13a = multiply(13,P1,A,n) print P13, P13a, P13 == P13aHi Programmingpraxis
Could you tell me if my implementation is correct or not specially for multiPoint function as i was trying to reduce one more call of mulitpOint function.
--y^2 = x^3+ax^2+b data Elliptic = Conelliptic Integer Integer Integer deriving (Eq,Show) data Point = Conpoint Integer Integer deriving (Eq,Show) addPoints::Point->Point->Point->Integer->Point addPoints (Conpoint x_1 z_1) (Conpoint x_2 z_2) (Conpoint x_0 z_0) n = let x_3= mod ((((x_1-z_1)*(x_2+z_2)+(x_1+z_1)*(x_2-z_2))^2)*z_0) n z_3= mod ((((x_1-z_1)*(x_2+z_2)-(x_1+z_1)*(x_2-z_2))^2)*x_0) n in (Conpoint x_3 z_3) doublePoint::Point->Elliptic->Point doublePoint (Conpoint x z) (Conelliptic an ad n) = (Conpoint (mod x' n) (mod y' n)) where x'= mod (4*ad*(x+z)^2*(x-z)^2) n y'= mod ((4*ad*(x-z)^2+t*an)*t) n where t=(x+z)^2-(x-z)^2 multiPoint::Integer->Integer->Point->Elliptic->Point multiPoint _ 0 _ _ = (Conpoint 0 0) multiPoint _ 1 pOint _ = pOint multiPoint n k pOint eCurve =let p=multiPoint n (div k 2) pOint eCurve in case even k of True -> doublePoint p eCurve _ -> addPoints p (addPoints p pOint pOint n) pOint n -- for 7p i am using 6p+pOint montgomeryCurve::Integer->Integer->[Integer] montgomeryCurve s n = let u=mod (s*s-5) n v=mod (4*s) n an=mod (((v-u)^3)*(3*u+v)) n ad=mod (4*u^3*v) n in [an,ad,n,mod (u*u*u) n,mod (v*v*v) n][…] previous example except using Montgomery curve rather than Lenstra curve. Programmingpraxis used Montgomery curve to find the prime factorisation. I personally feel that Lenstra algorithm worked fine when […]