Programming Praxis

We use a 2,3,5-wheel for the trial division before calling Hart’s algorithm:

(define (wheel-factors n limit)
  (let ((wheel (vector 1 2 2 4 2 4 2 4 6 2 6)))
    (let loop ((n n) (f 2) (next 0) (fs (list)))
      (cond ((< limit f) (values n fs))
            ((< n (* f f)) (values 1 (cons n fs)))
            ((zero? (modulo n f))
              (loop (/ n f) f next (cons f fs)))
            (else (loop n (+ f (vector-ref wheel next))
                        (if (= next 10) 3 (+ next 1)) fs))))))

We run Hart's algorithm without limit, so we complete the factorization no matter what; our variable ni is the product of n and i, since we never need i by itself:

(define (one-line-factor n)
  (let loop ((ni n))
    (let* ((s (isqrt ni))
           (s (if (= ni (* s s)) s (+ s 1)))
           (m (modulo (* s s) n))
           (t (isqrt m)))
      (if (= (* t t) m)
          (gcd (- s t) n)
          (loop (+ ni n))))))

Now we combine the two functions; the max 2 is necessary when n < 8 so n^1/3 < 2:

(define (factors n)
  (call-with-values
    (lambda () (wheel-factors n (max 2 (expt n 1/3))))
    (lambda (n fs)
      (if (= n 1) fs
        (let ((f (one-line-factor n)))
          (if (= f 1) (cons n fs)
            (cons (/ n f) (cons f fs))))))))

Worst-case for the algorithm occurs when n is prime. As a test we (slowly) compute the factors of the mersenne numbers:

> (do ((n 2 (+ n 1))) (#f)
    (display n) (display ":")
    (for-each
      (lambda (f) (display " ") (display f))
      (sort < (factors (- (expt 2 n) 1))))
    (newline))
2: 3
3: 7
4: 3 5
5: 31
6: 3 3 7
7: 127
8: 3 5 17
9: 7 73
10: 3 11 31
11: 23 89
12: 3 3 5 7 13
13: 8191
14: 3 43 127
15: 7 31 151
16: 3 5 17 257
17: 131071
18: 3 3 3 7 19 73
19: 524287
20: 3 5 5 11 31 41
21: 7 7 127 337
22: 3 23 89 683
23: 47 178481
24: 3 3 5 7 13 17 241
25: 31 601 1801
26: 3 2731 8191
27: 7 73 262657
28: 3 5 29 43 113 127
29: 233 1103 2089
30: 3 3 7 11 31 151 331

We used isqrt from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/fuTsFyTD.

Pages: 1 2