Programming Praxis

Here is our version of the elliptic curve factorization function:

(define ec-factor ; Crandall/Pomerance Algorithm 7.4.4
  (let* ((b1-saved 0) (b2-saved 0) (ps '()) (as '()) (ds '()) (d 210)
         (sv (make-vector (+ d 1) #f)) (bv (make-vector (+ d 1) #f)))
    (lambda (n b1 b2 s)
      (when (not (and (= b1 b1-saved) (= b2 b2-saved)))
        (set! b1-saved b1) (set! b2-saved b2) (set! ps (primes b2))
        (set! as (let loop ((ps ps) (as (list)))
                   (let ((a (ilog (car ps) b1)))
                     (if (= a 1) (append (reverse as) (cycle 1))
                       (loop (cdr ps) (cons a as))))))
        (set! ds (cons 0 (map (lambda (x) (/ x 2))
                   (map - (cdr ps) (reverse (cdr (reverse ps))))))))
      (let-values (((an ad q) (curve12 n s)))
        (let stage1 ((ps ps) (as as) (ds ds) (q q))
          (if (< (car ps) b1)
              (stage1 (cdr ps) (cdr as) (cdr ds)
                      (multiply (expt (car ps) (car as)) q an ad n))
              (let ((g (gcd (cdr q) n))) (if (< 1 g n) (list 1 g)
                (let* ((v (- b1 1)) (r (multiply v q an ad n))
                       (t (multiply (- b1 1 d d) q an ad n))
                       (alpha (modulo (* (car r) (cdr r)) n)))
                  (vector-set! sv 1 (double q an ad n))
                  (vector-set! sv 2 (double (vector-ref sv 1) an ad n))
                  (do ((i 3 (+ i 1))) ((< d i))
                    (vector-set! sv i
                      (add (vector-ref sv (- i 1))
                           (vector-ref sv 1)
                           (vector-ref sv (- i 2)) n)))
                  (do ((i 1 (+ i 1))) ((< d i))
                    (vector-set! bv i
                      (modulo (* (car (vector-ref sv i))
                                 (cdr (vector-ref sv i))) n)))
                  (let stage2 ((v v) (ps ps) (ds ds) (r r) (t t) (g 1))
                    (if (null? ps) (let ((g (gcd g n))) (if (< 1 g n) (list 2 g) #f))
                      (let ((v2d (+ v d d)) (alpha (modulo (* (car r) (cdr r)) n)))
                        (let loop ((ps ps) (ds ds) (g g))
                          (if (and (pair? ps) (< (car ps) v2d))
                              (loop (cdr ps) (cdr ds) (modulo (* g (+ (*
                                    (- (car r) (car (vector-ref sv (car ds))))
                                    (+ (cdr r) (cdr (vector-ref sv (car ds)))))
                                    (- alpha) (vector-ref bv (car ds)))) n))
                              (let* ((old-r r) (r (add r (vector-ref sv d) t n)))
                              (let ((g (gcd g n))) (if (< 1 g n) (list 2 g)
                                  (stage2 v2d ps ds r old-r g))))))))))))))))))

The only change we made was to replace variable d with i and variable r with v; since some implementations of Scheme don’t recognize differences in case, we couldn’t retain the original variable names. We reused the elliptic arithmetic functions add, double and multiply from the previous exercises.

Our ec-factor function is a drop-in replacement for the earlier factors function. An updated version of the function, with improved Pollard rho (the Brent variation) and Pollard p−1 (the two-stage improved standard continuation) functions, is shown on the next page. In this form, complicated factorizations are ten to fifty times faster than the old function.

The complete factors function is shown on the next page and at http://programmingpraxis.codepad.org/SdrzlWbh. Two sample factorizations are shown below:

> (factors (repunit 66))
Trial division: bound=1000
  Found factors (3 7 11 11 13 23 37 67)
  Remaining co-factor 58993628694531296147719957524688808978821367471004660551
Pollard rho: bound=500000, constant=1
  Found factor 4093, remaining co-factor 144132979952434146464011623563862225699
53913381628307
Pollard rho: bound=500000, constant=2
  Found non-prime factor 190056571
Trial division: bound=1000
Pollard rho: bound=500000, constant=1
  Found factor 21649, remaining co-factor 8779
  Factorization complete
Pollard rho: bound=500000, constant=3
  Found factor 513239, remaining co-factor 147761341791192537461208953746267179103
Pollard rho: bound=500000, constant=4
  Found factor 599144041, remaining co-factor 246620731710144034397913595783
Pollard rho: bound=500000, constant=5
Pollard p-1: b1=100000, b2=5000000
Elliptic curve 1: b1=50000, b2=2500000, s=1325338914
Elliptic curve 2: b1=50000, b2=2500000, s=1196705985
  In stage 2, found factor 183411838171, remaining co-factor 1344628210313298373
  Factorization complete
(3 7 11 11 13 23 37 67 4093 8779 21649 513239 599144041 183411838171 1344628210313298373)
> (factors (repunit 69))
Trial division: bound=1000
  Found factors (3 37 277)
  Remaining co-factor 3613722025274371844768956682314083036104696754516249101086646213
Pollard rho: bound=500000, constant=1
Pollard p-1: b1=100000, b2=5000000
  Found factor 203864078068831, remaining co-factor 17726134292546940492766254883101518975039516780123
Pollard p-1: b1=100000, b2=5000000
  Found factor 11111111111111111111111, remaining co-factor 1595352086329224644348978893
  Factorization complete
(3 37 277 203864078068831 11111111111111111111111 1595352086329224644348978893)

Pages: 1 2 3