Project Euler Problem 23

Ivar Thorson bio photo By Ivar Thorson

Trying to solve this problem the naive way quickly led me into computational problems…

(defn proper-divisors [n]
  "Returns a list of the proper divisors of n."
  (filter #(zero? (mod n %)) (range 1 (+ 1 (/ n 2)))))

(defn abundant? [n]
  (< n (reduce + (proper-divisors n))))

(def abundant-numbers (filter abundant? (range 1 20000)))

(defn all-pairs [coll]
  (for [i (range (count coll))
        j (range (inc i) (count coll))]
    [(nth coll i) (nth coll j)]))

(defn euler-23 []
  (let [sums (into #{} (map (fn [[a b]] (+ a b))
                            (all-pairs abundant-numbers)))
        not-sums (filter #(nil? (sums %)) (range 1 2000))]
    (reduce + not-sums)))

(euler-23)

After two minutes passed and no signs of stopping, I killed the process and decided to check how much time it was taking just to generate the abundant numbers.

(time (reduce + (filter abundant? (range 1 28123))))
"Elapsed time: 37798.750415 msecs"

Aha! It looks like it’s finally time to implement a more efficient version of proper-divisors, since that function has come up twice before, and is clearly one of the bottlenecks in this problem.

Back in problem 12 we used a trick to find the number of proper divisors from the number of primes, but stopped short of actually computing all of the proper divisors. It’s time for that to change.

The hypothesis that I am working under is that if we generate the list of proper divisors directly from prime numbers, it will be more efficient than finding the proper divisors through trial and error.

Here is the code, explanation to follow:

(defn all-pairs [coll]
  (for [i (range (count coll))
        j (range (inc i) (count coll))]
    [(nth coll i) (nth coll j)]))

(defn choose-k [coll k]
  "Returns all possible unique groups of size k in coll."
  (let [n (count coll)]
    (if (<= k 1)
      (map vector coll)
      (reduce into []
              (for [i (range 1 (inc n))]
                (map #(into [(nth coll (dec i))] %)
                     (choose-k (drop i coll) (dec k))))))))

(defn all-combos [coll]
  "Returns all possible unique combos of elements of coll."
  (reduce into []
          (map #(choose-k coll %) (range 1 (inc (count coll))))))

(defn proper-divisors-from-primes [num]
  (into #{1} (map #(apply * %) (all-combos (prime-factors-of num)))))

(defn abundant? [n]
  (< n (reduce + (proper-divisors-from-primes n))))

(def abundant-numbers (filter abundant? (range 1 28123)))

(time (reduce + abundant-numbers)) 
"Elapsed time: 7342.117583 msecs"

It was difficult for me to come up with the algorithm behind choose-k, and it took more than an hour of concentration to figure out exactly how to make it properly recursive. It’s an interesting exercise which I solved by first writing by-pairs, then writing by-triples, then by-quads, etc until the pattern became obvious enough to refactor the code into a recursive definition.

Once choose-k worked smoothly, everything else flowed out smoothly, and the performance of this way of generating primes is much better! Nothing like a 5-fold improvement to brighten my day.

Alas, this was not the only performance problem of the naive algorithm: generating all of those sums of abundant numbers had too heavy a performance toll. What is a better way?

This time, my curiosity to know the answer overwhelmed my discipline to always solve the problem myself first, and I admit that I peeked at the interesting solution that hircus wrote on clojure-euler. Putting the abundant numbers into a hashmap, and checking for pairs the way he did is very clever, and I reprint it here:

;; This is a clever one-liner! Thanks "hircus"
(defn sum-of-abundants? [i abundants]
  (some (fn [a] (abundants (- i a))) abundants))
 
(defn euler-23 []
  (let [abundants (into (sorted-set) (filter abundant? (range 12 28112)))]
    (reduce + (for [i (range 1 28124)
                    :when (not (sum-of-abundants? i abundants))]
                i))))

In the end, the performance of my approach was half as fast as the solution posted at clojure-euler (7.68secs vs 3.93secs), although creating a function able to generate all possible combinations of elements in a list – as I did in this exercise – may prove useful later. Hircus’s solution uses a tight loop and finds two divisors at a time, which gives good performance that is pretty hard to beat.