Lehman’s Factoring Algorithm

August 22, 2017

The explicit algorithms given by Crandall and Pomerance make this task easy:

```(define (wheel n) ; factoring by 2,3,5-wheel
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (f 2) (fs (list)) (w 0))
(if (< n (* f f)) (reverse (cons n fs))
(if (zero? (modulo n f))
(loop (/ n f) f (cons f fs) w)
(loop n (+ f (vector-ref wheel w)) fs
(if (= w 10) 3 (+ w 1))))))))```
```(define (fermat n) ; fermat's factoring algorithm
(call-with-current-continuation (lambda (return)
(do ((a (inexact->exact (ceiling (sqrt n))) (+ a 1)))
((square? (- (* a a) n)) =>
(lambda (b) (return (- a b)))))))))```
```(define (factors n) ; lehman's factoring algorithm
(define (wheel n limit)
(let ((wheel '#(1 2 2 4 2 4 2 4 6 2 6)))
(let loop ((n n) (f 2) (fs (list)) (w 0))
(if (< limit f) (values n fs)
(if (< n (* f f)) (values 1 (cons n fs))
(if (zero? (modulo n f)) (loop (/ n f) f (cons f fs) w)
(loop n (+ f (vector-ref wheel w)) fs (if (= w 10) 3 (+ w 1)))))))))
(define (lehman n)
(call-with-current-continuation (lambda (return)
(do ((k 1 (+ k 1))) ((exact (ceiling sqrt-4kn)) (+ a 1)))
((< (floor (+ sqrt-4kn (/ (expt n 1/6) (* 4 (sqrt k))))) a))
(let ((b (square? (- (* a a) 4kn))))
(if b (return (gcd (+ a b) n)))))))
(return n))))
(let ((limit (iroot 3 n)))
(call-with-values (lambda () (wheel n limit))
(lambda (n fs)
(if (= n 1) (reverse fs)
(let* ((p (lehman n)) (q (/ n p)))
(if (= q 1) (reverse (cons p fs))
(reverse (cons (max p q) (cons (min p q) fs))))))))))```

In both Fermat’s and Lehman’s algorithm, we use `call-with-current-continuation` to make an explicit return, so the relation to the pseudocode is clear. Because of the trial division, both the wheel and Lehman’s algorithm give a complete factorization, with all of the returned factors proven primes. Fermat’s algorithm, by contrast, finds prime factors, but the co-factor is not necessarily prime; for instance, `(fermat 165)` returns 11, which is prime, but the co-factor 15 is composite.

We did the following timing experiments:

```> (time (wheel 9428602644157043))
(time (wheel 9428602644157043))
no collections
0.677028356s elapsed cpu time
0.677727926s elapsed real time
64 bytes allocated
(48392053 194837831)
> (time (fermat 9428602644157043))
(time (fermat 9428602644157043))
1 collection
3.094562790s elapsed cpu time, including 0.000125445s collecting
3.123183707s elapsed real time, including 0.000129115s collecting
972992 bytes allocated, including 8354000 bytes reclaimed
48392053
> (time (factors 9428602644157043))
(time (factors 9428602644157043))
1 collection
0.041777502s elapsed cpu time, including 0.000139381s collecting
0.041849967s elapsed real time, including 0.000142982s collecting
7188208 bytes allocated, including 8600032 bytes reclaimed
(48392053 194837831)```

The first algorithm made 12,904,551 trial divisions of 2,3,5-pseudoprimes and took two-thirds of a second. Fermat’s algorithm started with a = 97100992 and made 24,513,950 passes through an inner loop more complicated than the trial division algorithm, ending with a = 121614942 and b = 73222889, taking more than three seconds. Lehman’s algorithm tested 56,339 trial divisors and tried multipliers from 1 to 5814 and computed a total of 17,424 as, an average of about three a per k, and took less than a twentieth of a second. Both trial division and Fermat’s algorithm have an order of growth proportional to the square root of n, but trial division performs its inner loop half as many times as Fermat’s algorithm, and its inner loop is much less complicated, giving it a huge speed advantage.

As a further test, I calculated the factors of repunits Rn = (10n − 1) / 9 in succession, starting from n= 2. After about four hours, it reached R31. That’s good, but Pollard’s rho algorithm made the same computation in less than a second, which explains why Lehman’s algorithm is not used in practice.

You can run the program at http://ideone.com/b3AN9T.

Pages: 1 2

2 Responses to “Lehman’s Factoring Algorithm”

1. Paul said

In Python.

```def td_factor(n, limit=10**6):
'stop after first factor found'
for f in accumulate(chain((2, 1, 2, 2), cycle((4, 2, 4, 2, 4, 6, 2, 6)))):
if f*f > n or f > limit:
return 0
if n % f == 0:
return f

def fermat(n):
for a in range(isqrt(n)+1, (n + 9) // 6 + 1):
b = a * a - n
if is_square(b):
return a - isqrt(b)
return 0

def lehman(n):
s3n = ceil(n ** (1/3))
f = td_factor(n, limit=s3n)
if f:
return f
for k in range(1, s3n + 1):
sk = 2*sqrt(k*n)
for a in range(ceil(sk), floor(sk + n**(1/6) / (4*sqrt(k)))+1):
b = a*a - 4*k*n
if is_square(b):
return gcd(a + isqrt(b), n)
return 0
```
2. John said

I shall show this to my father, he celebrated his 89th birthday yesterday. He will be pleased that the algorithm that carries his name is still in use.