## 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, )).

divisors([], Divisors) -> Divisors;
divisors([{P, R}|Factors], Divisors) ->
divisors(Factors, [D * imath:pow(P, E) || D <- Divisors, E <- seq(0, R)]).
```