Mathematicians are an odd bunch. Names for just about everyhing. There are amicable numbers and perfect number, sociable number and betrothed numbers. There are sublime number, frugal number, and quasiperfect number. Heck, there are powerful number, smooth number, and even sphenic numbers. Rather a lot to deal with all told… So let’s just focus on two of them: perfect numbers and amicable numbers.

The two are rather related, in that they’re both defined in terms of their sum of proper divisors. Start with divisors: the set of numbers that evenly divide a given number. So say the divisors of 25 are 1, 5, and 25. Proper adds the constraint that a number is not a divisor of itself (so only 1 and 5). Sum them and you have spd(25) = 6 .

If spd(n) = n , a number is considered perfect. If spd(n) = m and spd(m) = n , the m and n are amicable. Simple¹ as that. The same context can be extended to a sequence:

So that’s what we’re dealing with today, given Tuesday’s post on Programming Praxis. In a single sentence: find all perfect numbers, amicable pairs, and amicable chains less than one million.

Let’s get started.

Perfect numbers and amicable pairs

First we need to write a function that can quickly calculate the sum-of-divisors of a number. We’ll be using the math/number-theory module to calculate divisors and jbclements’ memoize module to save us a bit of time calculating a pile of these over and over again:

(require math/number-theory
         memoize)

; Calculate sum of proper divisors (proper, thus subtracting i)
(define/memo (sum-of-divisors n)
  (- (apply + (divisors n)) n))

Testing it out with the value we already know:

> (sum-of-divisors 25)
6

From here and with the definitions above, we can directly define perfect numbers:

; A number is perfect if its divisors sum to itself
(define (perfect? n)
  (= n (sum-of-divisors n)))

And amicable pairs:

; Two numbers are amicable if each numbers's sum of divisors is the other
(define (amicable? m n)
  (and (= m (sum-of-divisors n))
       (= n (sum-of-divisors m))))

As an example:

> (perfect? 6)
#t
> (perfect? 25) ; (spd 25) -> 6
#f
> (amicable? 25 6)
#f
> (amicable? 220 284) ; (spd 220) -> 284; (spd 284) -> 220
#t

Looks good. But how about we think a little bit bigger, eh?

Amicable chains

For the next part of the problem, we want to calculate amicable chains. Sequences of numbers of any length such that the sum of proper divisors of each number is the next and the chain wraps around at the end. On the smaller end, all perfect numbers are amicable chains of length 1. All amicable pairs are chains of length 2. One example amicable chain of length 5 would be (12496 14264 14536 15472 14288)².

Given all of that, we should be able to write a nice recursive function to determine if a number is in an amicable chain and, if so, return it. To do so, we’re going to essentially need to recur down a chain of numbers until one of three/four conditions is met:

• We find the original number in the sequence; return the amicable chain
• The chain terminates at 1 (any prime number will have 1 as the sum of divisors after which the chain will terminate)
• We find a number other than the original for the second time; there is a chain, but the number in question is not part of it (example 25 -> 6 which then loops)
• optional: The current number is larger than some bound; it may be a chain, but for our purposes return #f

Translating this relatively directly into code, we have:

; An amicable chain is a sequence of numbers where each's sum of divisors is the next
; Ignore chains that leave the given bounds
(define (amicable-chain n [bound +inf.0])
  (let/ec return
    (let loop ([x (sum-of-divisors n)] [prev '()])
      (cond
        [(= x n)
         (cons x prev)]
        [(or (> x bound) (= x 1) (member x prev))
         (return #f)]
        [else
         (loop (sum-of-divisors x)
               (cons x prev))]))))

As given earlier:

> (amicable-chain 12496)
'(12496 14264 14536 15472 14288)
> (amicable-chain 15472)
'(15472 14288 12496 14264 14536)

That gives us most of our framework. All that’s left is extended it to a bunch of numbers at a time (and as a bonus: making a pretty picture).

Generating a graph

Basic idea: Iterate through all numbers 2 through some bound. For each in turn, if it’s an amicable chain, add the numbers in the chain to a graph. Skip any number added as we iterate upwards.

; Create a graph of all amicable chains within the bounds
(define (amicable-chains-graph bound)
  (define graph (unweighted-graph/directed '()))
  (define colors (make-hash))

  (for ([i (in-range 2 bound)] #:unless (has-vertex? graph i))
    (cond 
      [(amicable-chain i bound) 
       => (λ (chain)
            (displayln chain)
            (for ([from (in-list chain)] 
                  [to (in-list (snoc (car chain) (cdr chain)))])
              (hash-set! colors from
                         (case (length chain)
                           [(1) 0] [(2) 1] [else 2]))
              (add-directed-edge! graph from to)))]))

  (values graph colors))

That’s actually all there is to it. One interesting trick is iterating over the list chain and (snoc (car chain) (cdr chain)). Basically, that’s a sneaky way to iterate over all sequential numbers in the list. snoc is the reversed form of cons. Rather than adding to the head of a list, add to the end. It’s terribly inefficient, but at the very least it serves our purposes.

With a helper wrapping the graphviz function to actually generate the files, and we can actually write out our results:

Looks like we have four perfect numbers (in red; 6, 28, 496, and 8128), a whole pile of amicable pairs, and two nice chains (the larger has almost 30 numbers). I’m actually a little surprised that such a large chain already appears under a million. It sort of makes me wonder just how large of chains there might be out there if you look further…

And that’s it for today. As always, the code is on GitHub: amicable-chain.rkt. Enjoy!