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

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

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?

It’s used in cryptography, for one thing.

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.