A Divisor Apology
September 9, 2016
Let’s begin with the divisors
function that I have been using up until now:
(define (divisors1 n) ; divisors of n, including 1 and n (let ((divs (list 1))) (do ((fs (factors n) (cdr fs))) ((null? fs) (sort < divs)) (let ((temp (list))) (do ((ds divs (cdr ds))) ((null? ds) (set! divs (append divs temp))) (let ((d (* (car fs) (car ds)))) (when (not (member d divs)) (set! temp (cons d temp))))))))) > (time (length (divisors1 (apply * (primes 73))))) (time (length (divisors1 (...)))) 20 collections 17399.347438981s elapsed cpu time, including 0.274695336s collecting 17629.483049516s elapsed real time, including 0.588428532s collecting 200539040 bytes allocated, including 54712832 bytes reclaimed 2097152
That’s nearly five hours, and simply unacceptable. The problem is the call to member
, which performs a linear search through the growing list of divisors; that’s fine when the number of divisors is small, but is a killer for the highly composite numbers that we have been working with in recent exercises. How could I possibly write code that bad, and not notice it for many months? Mea culpa, mea culpa, mea maxima culpa. A better approach calculates the divisors in a pair of nested loops, one in the do
and one in the map
:
(define (divisors2 n) (let ((divs (list 1))) (do ((fs (factors n) (cdr fs))) ((null? fs) (unique = (sort < divs))) (set! divs (append (map (lambda (d) (* d (car fs))) divs) divs))))) > (time (length (divisors2 (apply * (primes 73))))) (time (length (divisors2 (...)))) 30 collections 1.341846042s elapsed cpu time, including 0.415476502s collecting 1.392317452s elapsed real time, including 0.425363038s collecting 314199616 bytes allocated, including 549322704 bytes reclaimed 2097152
That’s three and a half orders of magnitude better. Moving the calls to unique
and sort
inside the loop improves things even more by reducing the size of the intermediate lists:
(define (divisors3 n) (let ((divs (list 1))) (do ((fs (factors n) (cdr fs))) ((null? fs) divs) (set! divs (unique = (sort < (append (map (lambda (d) (* d (car fs))) divs) divs))))))) > (time (length (divisors3 (apply * (primes 73))))) (time (length (divisors3 (...)))) 48 collections 0.918762104s elapsed cpu time, including 0.554705330s collecting 0.961626818s elapsed real time, including 0.583698212s collecting 500117920 bytes allocated, including 478435328 bytes reclaimed 2097152
Our final version notes that the sort is unnecessary, because the factors are generated in order, and the map
retains that order. That means we can merge, an O(n) operation, instead of sort, an O(n log n) operation. Further, we can save a pass over the data by combining the merge with the unique operation:
(define (merge-uniq xs ys) (let loop ((xs xs) (ys ys) (zs (list))) (cond ((and (null? xs) (null? ys)) (reverse zs)) ((null? xs) (loop xs (cdr ys) (cons (car ys) zs))) ((null? ys) (loop (cdr xs) ys (cons (car xs) zs))) ((< (car xs) (car ys)) (loop (cdr xs) ys (cons (car xs) zs))) ((< (car ys) (car xs)) (loop xs (cdr ys) (cons (car ys) zs))) (else (loop xs (cdr ys) zs))))) (define (divisors4 n) (let loop ((fs (factors n)) (divs (list 1))) (if (null? fs) divs (loop (cdr fs) (merge-uniq (map (lambda (d) (* d (car fs))) divs) divs))))) > (time (length (divisors4 (apply * (primes 73))))) (time (length (divisors5 (...)))) 19 collections 0.441666600s elapsed cpu time, including 0.264594108s collecting 0.465659087s elapsed real time, including 0.278028872s collecting 207327184 bytes allocated, including 117907312 bytes reclaimed 2097152
Less than half a second. I can live with that. You can run the program at http://ideone.com/8HI7m2.
How does your program perform with, say, 278914005382139703576000?
The answer is, not too bad (I was wondering about what would happen with lots of repeated factors resulting in the merge unique having to do more work):
Erlang can calculate the divisors of those rather large numbers within 2 – 4 seconds. (matthew’s smaller number 1.6 seconds; 4 seconds for the bigger number.) I’m not claiming any credit, since the implementation is just a simple and rather naive list comprehension.