Problem 021: Amicable Numbers

It's been a long time since I've posted! That's alright though, we're back into it now. From Project Euler:

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.

$n!$ means $n × (n − 1) × ... × 3 × 2 × 1$

There's a lot of potential optimizations for this one, but I found it works just fine by brute forcing it. We'll start by defining a function that sums the divisors of a number:

(defn sum-divisors
  "Sums the divisors for the provided number."
  [n]
  (assert (> 0 n))
  ; Loop through all the numbers from 2 to sqrt(n). We only need to
  ; go to sqrt(n) because all divisors less than that will have a
  ; corresponding divisor greater than sqrt(n).
  (loop [current 2
         ; Start total from 1 because 1 is always a divisor.
         total 1]
    (if (> current (Math/sqrt n))
      total
      (recur (+ 1 current)
             (if (= 0 (mod n current))
               ; In the case of a square number, we only want to
               ; add the number once.
               (if (= current (quot n current))
                 (+ total current)
                 ; Here we're not square, so we add both the number
                 ; and it's complement on the other
                 ; side of sqrt(n).
                 (+ total current (quot n current)))
               total)))))

After that, our code is pretty straight-forward:

(defn sum-amicable-numbers
  "Sums all amicable numbers up to and including `max`."
  [max]
  (->> (range 1 (+ 1 max))
       (map #(list % (sum-divisors %)))
       ; Prune out numbers that are amicable with themselves.
       (remove #(= (first %) (last %)))
       ; Filter out any non-amicable numbers.
       (filter #(= (first %) (sum-divisors (last %))))
       ; And sum them all up.
       (map first)
       (reduce +)))

This runs very fast:

$ lein run
Processing...
Total is: 41274
"Elapsed time: 153.883368 msecs"

Problem 012: Highly divisible triangle number

From Project Euler:

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Let us list the factors of the first seven triangle numbers:

1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.

What is the value of the first triangle number to have over five hundred divisors?

A way to solve this problem is fairly straight-forward: progress through all the triangle numbers, figure out their divisors, and count them. A simple brute force method is to just enumerate all the numbers up to $\sqrt{n}$, and for each one that divides $n$, add 2 to our count (one for this number, and one for $n$ divided by it). The corner case is if $n$ is square, just count the square root once:

(defn num-divisors [n]
  (loop [current 1
         total 0]
    (cond
      (> current (Math/sqrt n)) total
      (== current (Math/sqrt n)) (+ 1 total)
      :else (recur
              (+ 1 current)
              (+ total
                 (if (== 0 (mod n current)) 2 0))))))

We can then do a simple search through all the triangle numbers to get the first one that has at least 500 divisors:

(defn search-brute-force [limit]
  (loop [n 1
         triangle 1]
    (let [d (num-divisors n)]
     (if (> d limit)
       triangle
       (recur (+ 1 n)
              (+ triangle (+ 1 n)))))))

Unfortunately, this goes well past the one minute limit. Let's see if we can speed this up a bit. Take the formula for a triangle number $T_n$:
$$
T_n = \sum_{i=1}^n = \frac{n(n + 1)}{2}
$$
From this formula we can infer a few other interesting details. Whenever you have two sequential numbers one of them will be even and the other will be odd. Any even number will be divisible by 2, so the 2 in the denominator will get cancelled out by that and we'll be left with a divisor of $T_n$. The odd number left over will also be a divisor.

To write this out more explicitly, if $n$ is even then $\frac{n}{2}$ and $n + 1$ are divisors of $T_n$; if $n$ is odd then $n$ and $\frac{n + 1}{2}$ are divisors of $T_n$.

Next we'll want to dig into the properties of divisors. Any number can be broken down into a product of prime numbers, and some of those primes might be duplicated (for example, 4 = 2 * 2). So we can write any number $n$ like this:
$$
n = p_1^{c_1} * p_2^{c_2} * ...
$$
Then for each prime $p_i$, there's going to be $c_i$ factors of $n$ from that prime: $p_i$, $p_i^2$, $p_i^3$, all the way up to $c_i$. Each one of those can multiply with every other factor of $n$ to produce another factor (for example with $n = 28$ you can have $2 * 7$ to produce 14). So the number of divisors $d(n)$ is $c_i + 1$ multiplied by all the other possibilities:
$$
d(n) = (c_1 + 1) * (c_2 + 1) * ...
$$
Lastly, if you have two numbers $a$ and $b$ that are co-prime with one another, this means that they do not share any divisors other than 1. Therefore none of the $p_i$s are the same, so the number of divisors is just the product of the number of divisors for $a$ and $b$ separately. Since $n$ and $n + 1$ are necessarily coprime, we can just multiply the number of divisors of $n$ with the number of divisors of $n + 1$, dividing the even one by 2 to remove the denominator in the formula for $T_n$:
$$
d(T_n) = \left\{
\begin{array}{ll}
d(\frac{n}{2}) * d(n + 1) & n\ even \\
d(n) * d(\frac{n + 1}{2}) & n\ odd
\end{array}
\right.
$$
These two numbers should be quite a lot smaller than $T_n$, so the search time will be much faster.

Let's write out a little function to calculate the number of divisors of a triangle number:

(defn num-triangle-divisors [n]
  (if (even? n)
    (* (num-divisors (quot n 2))
       (num-divisors (+ 1 n)))
    (* (num-divisors n)
       (num-divisors (quot (+ 1 n) 2)))))

And then tweak our search function to use this instead:

(defn search-triangle [limit]
  (loop [n 1
         triangle 1]
    (let [d (num-triangle-divisors n)]
     (if (> d limit)
       triangle
       (recur (+ 1 n)
              (+ triangle (+ 1 n)))))))

This one is much faster, even with the 3s or so it takes for lein to start up:

$ time lein run
76576500

real 0m3.848s
user 0m4.687s
sys 0m0.355s

Problem 005: Smallest Multiple

From Project Euler:

2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?

This is a pretty straight-forward problem: we're just trying to find the least-common multiple (LCM) of a set of numbers. To find the LCM, it is just the product of the two numbers $a$ and $b$ divided by their greatest common divisor (GCD). If we think for a second about why this works: any number is just a product of prime factors. A prime number is just a product of itself and 1, composite numbers are a more complex sequence: 12 is a product of 2 * 2 * 3, 8 is a product of 2 * 2 * 2. When you multiply the two numbers together, you're a making a long sequence of prime factors that will result in the product of $a$ and $b$.

This product however is not the least common multiple, because there could be factors in common between $a$ and $b$: both 8 and 12 have a common factor of 2 * 2, if we simply multiplied them together we'd get 96 when the LCM is actually 24. We should divide by the GCD (2 * 2) to filter out these common factors. This gives us the correct result: 96 / 4 = 24.

Here's the code, which uses the Euclidean algorithm to find the GCD of two numbers:

(defn gcd [a b]
  (cond 
    (< a b) (gcd b a)
    (== 0 b) a
    :else (gcd b (mod a b))))

(defn lcm [a b]
  (if (or (== a 0) (== b 0))
    0
    (quot (* a b) (gcd a b))))

Then we can just do a reduce to get the LCM of the range:

(println (reduce lcm (range 1 21)))

This runs pretty fast (note that on my machine the JVM takes about 1.5s to start up):

$ time clj 5.clj
232792560

real 0m1.585s
user 0m3.136s
sys 0m0.151s

An alternative approach is to use something like the Sieve of Eratosthenes (SoT). This is an algorithm for finding prime numbers, but we can use the same idea here. The way the SoT works is you take a list of numbers from 2 to $N$. For each number $i$, you loop through the rest of the list, removing any multiple of $i$. After the first pass we'll have removed all numbers divisible by 2, after the second pass all numbers divisible by 3, etc. What we'll be left with is just a list of prime numbers. We'll apply the same idea here, but instead of removing the number if it is divisible, we'll divide out our current factor.

(defn lcm-sieve [xs]
  (if (empty? xs) ()
    (let [head (first xs)
          tail (rest xs)]
      (cons head
            (lcm-sieve (map #(if (== 0 (mod % head))
                               (quot % head)
                               %)
                              tail))))))

This algorithm efficiently gives us the list of prime factors for the LCM:

user=> (lcm-sieve (range 1 11))
(2 3 2 5 1 7 2 3 1)

These can be combined using reduce:

(println (->> (range 1 21)
              (lcm-sieve)
              (reduce *)))

This gives the correct result:

rob@alien ~/code/euler $ time clj 5.clj
232792560

real 0m1.549s
user 0m3.014s
sys 0m0.143s

Both methods seem quite fast, I hit integer overflow on my machine without either of them taking more than half a second. The reason I dove into the second one is that the SoT algorithm is the basis for a lot of my solutions in the future, so it's good to see a version of it now.

Problem 001: Multiples of 3 and 5

From Project Euler:

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

For this one, it's pretty simple just to iterate through all the numbers and add them up as needed (I use threading syntax here because I think it's more readable):

user=> (->> (range 1 1000)
            (filter #(or (== 0 (mod % 3)) (== 0 (mod % 5))))
            (reduce +))
233168

However we don't even need to loop. Since we care only about the sum and not the numbers themselves, we can apply a little bit of number theory. What we want is the sum of all the multiples of 3 or 5. This is equal to the sum of all the multiples of 3, plus the sum of all the multiples of 5, minus the sum of all the multiples of 15 to avoid double-counting.

We can rewrite the sum of the multiples of 3 from 1 to 999 (because we don't want to include 1000) as 3 times the sum of all numbers from 1 to 999 / 3, or 333. There's a formula to get the sum of numbers from 1 to n:
$$
\sum^n_{i=1} i = \frac{n(n+1)}{2}
$$
Which we can write in Clojure like this:

(defn sumto [n] (quot (* n (+ n 1)) 2))

To get the sum of the multiples of 3 from 1 to 1000, we can just write:

(* 3 (sumto (quot 999 3)))

Then to get the final result, it's simply a matter of combining everything together:

user=> (+ (* 3 (sumto (quot 999 3)))
          (* 5 (sumto (quot 999 5)))
          (- (* 15 (sumto (quot 999 15)))))
233168