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

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 )

Google photo

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

Twitter picture

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

Facebook photo

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

Connecting to %s

%d bloggers like this: