; elliptic curves
(define (ecm-plus ecm p1 p2)
(define a (car ecm)) (define b (cadr ecm)) (define m (caddr ecm))
(define x car) (define y cadr) (define z caddr)
(define (sq x) (* x x)) (define infinity (list 0 1 0))
(define (inv x) ; modular inverse
(let loop ((x (modulo x m)) (a 1))
(cond ((zero? x) #f) ((= x 1) a)
(else (let ((q (- (quotient m x))))
(loop (+ m (* q x)) (modulo (* q a) m)))))))
(define (return n d)
(let ((d (inv d)))
(if (not d) #f
(let ((x3 (modulo (- (* n d n d) (x p1) (x p2)) m)))
(list x3 (modulo (- (* n d (- (x p1) x3)) (y p1)) m) 1)))))
(cond ((or (not (pair? p1)) (not (pair? p2))) #f)
((zero? (z p1)) p2) ((zero? (z p2)) p1)
((= (x p1) (x p2))
(if (zero? (modulo (+ (y p1) (y p2)) m)) infinity
(return (+ (* 3 (sq (x p1))) a) (* 2 (y p1)))))
(else (return (- (y p2) (y p1)) (- (x p2) (x p1))))))
(define (ecm-negate ecm p)
(define x car) (define y cadr) (define z caddr)
(list (x p) (modulo (- (y p)) (caddr ecm)) (z p)))
(define (ecm-double ecm p) (ecm-plus ecm p p))
(define (ecm-minus ecm p1 p2) (ecm-plus ecm p1 (ecm-negate ecm p2)))
(define (ecm-times ecm n p)
(cond ((= n 1) p)
((odd? n) (ecm-plus ecm p (ecm-times ecm (quotient n 2) (ecm-double ecm p))))
(else (ecm-times ecm (quotient n 2) (ecm-double ecm p)))))
(display (ecm-plus '(4 4 5) '(1 2 1) '(4 3 1))) (newline)
(display (ecm-negate '(4 4 5) '(1 2 1))) (newline)
(display (ecm-double '(4 4 5) '(1 2 1))) (newline)
(display (ecm-minus '(4 4 5) '(4 2 1) '(4 3 1))) (newline)
(display (ecm-times '(4 4 5) 3 '(1 2 1))) (newline)
(display (ecm-times '(4 4 5) 4 '(1 2 1))) (newline)
(display (ecm-times '(4 4 5) 5 '(1 2 1))) (newline)