; elliptic curve factorization

(define rand
  (let* ((a 3141592653) (c 2718281829)
         (m (expt 2 35)) (x 5772156649)
         (next (lambda ()
                 (let ((x-prime (modulo (+ (* a x) c) m)))
                   (set! x x-prime) x-prime)))
         (k 103)
         (v (list->vector (reverse
              (let loop ((i k) (vs (list x)))
                (if (= i 1) vs
                  (loop (- i 1) (cons (next) vs)))))))
         (y (next))
         (init (lambda (s)
                 (set! x s) (vector-set! v 0 x)
                 (do ((i 1 (+ i 1))) ((= i k))
                   (vector-set! v i (next))))))
    (lambda seed
      (cond ((null? seed)
              (let* ((j (quotient (* k y) m))
                     (q (vector-ref v j)))
                (set! y q)
                (vector-set! v j (next)) (/ y m)))
            ((eq? (car seed) 'get) (list a c m x y k v))
            ((eq? (car seed) 'set)
              (let ((state (cadr seed)))
                (set! a (list-ref state 0))
                (set! c (list-ref state 1))
                (set! m (list-ref state 2))
                (set! x (list-ref state 3))
                (set! y (list-ref state 4))
                (set! k (list-ref state 5))
                (set! v (list-ref state 6))))
            (else (init (modulo (numerator
                    (inexact->exact (car seed))) m))
                  (rand))))))

(define (randint . args)
  (cond ((null? (cdr args))
          (inexact->exact (floor (* (rand) (car args)))))
        ((< (car args) (cadr args))
          (+ (inexact->exact (floor (* (rand) (- (cadr args) (car args))))) (car args)))
        (else (+ (inexact->exact (ceiling (* (rand) (- (cadr args) (car args))))) (car args)))))

(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-factor n)
  (let* ((a (randint n)) (b (- a)) (ecm (list a b n)) (p '(1 1 1)))
    (let loop ((p p) (q (ecm-plus ecm p p)))
      (if (or (not q) (equal? '(0 1 0) q))
          (let ((g (gcd (- (car p) 1) n)))
            (if (= g 1) (ecm-factor n) g))
          (loop q (ecm-plus ecm q '(1 1 1)))))))

(display (ecm-factor 455839))

761