Project Euler: Problem 21

The last two Euler problems were rather anticlimactic, so let’s just keep going.

Let d(n) be defined as the sum of proper divisors of n (numbers less than n which divide evenly into n).
If d(a) = b and d(b) = a, where a ≠ b, then a and b are an amicable pair and each of a and b are called amicable numbers.

For example, the proper divisors of 220 are 1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284. The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

Evaluate the sum of all the amicable numbers under 10000.

We’ll solve this problem in F#. The goal is to identify amicable numbers, but we start off with a few simple helper functions we’ll need:

let Divisors n = [1..n / 2] |> List.filter (fun x -> n % x = 0)

let SumDiv n = Divisors n |> List.sum

These two provide us with a simple implementation of d(n). The next step is using this to identify the amicable numbers. The definition provided in the problem description depends on amicable pairs. From it, we’ll derive the definition of an amicable number. First, we have d(a) = b and d(b) = a. Rewriting, we get d(d(a)) = a. We also have a ≠ b. Using the previous definition of d(a) = b, we can rewrite this to a ≠ d(a).

Hence, a number n is an amicable number if d(d(n)) = n and d(n) ≠ n. Using this definition, we can go through all the numbers under 10000 and filter out the amicable pairs. We then take the sum to arrive at the answer.

let result = [1..9999] |>
                List.filter (fun x -> SumDiv x <> x && SumDiv(SumDiv x) = x) |>
                List.sum

This solution is relatively straightforward, but also inefficient. It takes several seconds to return the answer, which is acceptable but also clearly non-optimal. Indeed, our Divisors function uses an inefficient brute-force technique which is slowing down response times. Let’s see if we can optimize this.

First, rather than filtering an integer list to get our divisors, we use a bitmap to represent them: a boolean array m for integer n where each position m[i] indicates whether or not i is a proper divisor for n. We compute the bitmap as follows:

let DivMap n =
    if n = 1 then
        [| false; true |]
    else
        let m = Array.create((n / 2) + 1) false;
        m.[1] <- true
        for i in 2 .. m.Length - 1 do
            if m.[i] = false && n % i = 0 then
                m.[i] <- true
                m.[n / i] <- true
        m

Although more complicated, this function computes the divisors of a number significantly faster than our previous naive approach. Next, we modify our SumDiv function to use it:

let SumDiv n =
    let mutable s = 0
    let m = DivMap n
    for i in 1 .. m.Length - 1 do
        if m.[i] = true then
            s <- s + i
    s

Overall these two changes get our response time down to less than a second. Still not instant, but good enough for me.