Lenstra’s Algorithm

August 4, 2009

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)))
      (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

4 Responses to “Lenstra’s Algorithm”

  1. robert said

    You say this is of immense practical value, would you mind expanding on what practical applications this has? I understand the mathematical problem it solves, but what kinds of real-world problems require this kind of computation?

  2. Oisín said

    It’s used in cryptography, for one thing.

  3. Joel said

    Primes are the basis for most keys in modern cryptography algorithms, so if you want to break those algorithms, one way to do so is to factor the primes and get the keys back.

  4. nvcytfh said

    Then you can check if you were lucky; if gcd(d,n) > 1, d is a factor
    I think you meant gcd(d, n) is the factor

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: