Programming Praxis

Our version of Lenstra’s algorithm follows the pseudo-code exactly, searching a single curve for factors up to the limit:

(define (lenstra n limit)
  (call-with-current-continuation (lambda (return)
    (let* ((x (randint n)) (y (randint n)) (a (randint n))
           (b (- (* y y) (* x x x) (* a x)))
           (d (+ (* 4 a a a) (* 27 b b))) (g (gcd d n))
           (ps (primes limit)))
      (define (inverse x)
        (let loop ((x (modulo x n)) (a 1))
          (cond ((zero? x) (return #f))
                ((= x 1) a)
                (else (let ((q (- (quotient n x))))
                        (loop (+ n (* q x)) (modulo (* q a) n)))))))
      (define (add p1 p2)
        (let* ((x1 (car p1)) (y1 (cdr p1))
               (x2 (car p2)) (y2 (cdr p2))
               (g (gcd (- x1 x2) n)))
          (if (zero? g) (return #f)) (if (< 1 g n) (return g))
          (let* ((slope (* (inverse (- x1 x2)) (- y1 y2)))
                 (newx (modulo (- (* slope slope) x1 x2) n))
                 (newy (modulo (- (* slope (- x1 newx)) y1) n)))
            (cons newx newy))))
      (define (double p)
        (let* ((x (car p)) (y (cdr p)) (g (gcd (+ y y) n)))
          (if (zero? g) (return #f)) (if (< 1 g n) (return g))
          (let* ((slope (* (inverse (+ y y)) (+ (* 3 x x) b)))
                 (newx (modulo (- (* slope slope) x x) n))
                 (newy (modulo (- (* slope (- x newx)) y) n)))
            (cons newx newy))))
      (define (mult k p)
        (cond ((= k 1) p)
              ((even? k) (mult (quotient k 2) (double p)))
              (else (add (mult (quotient k 2) (double p)) p))))
      (if (zero? g) (return #f)) (if (< 1 g n) (return g))
      (let loop ((p (car ps)) (q (car ps)) (ps (cdr ps)) (xy (cons x y)))
        (cond ((null? ps) (return #f))
              ((< limit q) (loop (car ps) (car ps) (cdr ps) xy))
              (else (loop p (* p q) ps (mult p xy)))))))))

The only thing that might be unfamiliar is call-with-current-continuation, which sets up an early exit when a factor is found; calling (return x) causes lenstra to immediately stop what it is doing and return x. Factor counts curves until it finds a factor:

(define (factor n limit curves)
  (let loop ((curves curves))
    (if (zero? curves) #f
      (let ((f (lenstra n limit)))
        (if f f (loop (- curves 1)))))))

Factors calls factor repeatedly until it finds all the factors of the original number, using trial-factors to “prime” the pump using trial division:

(define (trial-factors n ps)
  (let loop ((ps ps) (facts '()) (cofact n))
    (cond ((or (null? ps) (= cofact 1))
            (values (reverse facts) cofact))
          ((< cofact (* (car ps) (car ps)))
            (values (reverse (cons cofact facts)) 1))
          ((zero? (modulo cofact (car ps)))
            (loop ps (cons (car ps) facts)
                  (quotient cofact (car ps))))
          (else (loop (cdr ps) facts cofact)))))

(define (factors n . args) ; factors n limit curves trial-limit
  (let ((limit (if (< 1 (length args)) (car args) 10000))
        (curves (if (< 2 (length args)) (cadr args) 1000))
        (trial-limit (if (< 3 (length args)) (caddr args) #e1e6)))
    (call-with-values
      (lambda () (trial-factors n (primes trial-limit)))
      (lambda (facts cofact)
        (let fact ((n cofact) (facts facts))
          (cond ((= n 1) (sort < facts))
                ((prime? n) (sort < (cons n facts)))
                (else (let ((f (factor n limit curves)))
                        (cond ((not f)
                                (display "failed with factors")
                                (for-each (lambda (f)
                                            (display " ") (display f))
                                          (sort < facts))
                                (display " and co-factor ")
                                (display n) (newline))
                             ((prime? f) (fact (/ n f) (cons f facts)))
                             (else (let ((fs (fact (/ n f) (list f))))
                                     (fact (apply / n fs) (append fs facts)))))))))))))

We use sort, rand, and randint from the Standard Prelude, primes from the exercise on the Sieve of Eratosthenes and prime? from the exercise on Miller/Rabin primality checking. You can run the code at http://programmingpraxis.codepad.org/S6l8kId7.

Pages: 1 2