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.

Pages: 1 2

3 Responses to “A Divisor Apology”

  1. matthew said

    How does your program perform with, say, 278914005382139703576000?

  2. matthew said

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

    > (time (length (divisors4 278914005382139703576000)))
    cpu time: 471 real time: 472 gc time: 359
    1032192
    > (time (length (divisors4 40729680599249024150621323470)))
    cpu time: 1052 real time: 1054 gc time: 814
    2097152
    
  3. David said

    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.

    % Divisors(n)
    % - Return a list of all of the divisors of n
    %
    divisors(N) when is_integer(N) -> divisors(factorize(N));
    
    divisors(L) when is_list(L) -> sort(divisors(L, [1])).
    
    divisors([], Divisors) -> Divisors;
    divisors([{P, R}|Factors], Divisors) ->
        divisors(Factors, [D * imath:pow(P, E) || D <- Divisors, E <- seq(0, R)]).
    

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: