Clojure Companion Cube

My solution to Project Euler Problem 72 is rather simple. We are asked to find the length of a Farey Sequence for all fractions with denominators less than one million.

Wikipedia’s page on the Farey Sequence gives us the property we need to solve this problem: a relationship between Euler’s Totient function and the length of a Farey Sequence. Reusing the euler-totient function from problem 69 and the solution was just a few lines:

(defn farey-seq-length [n]
  (reduce + (map euler-totient (range 1 (inc n)))))

(time (farey-seq-length 1000000)) ;; "Elapsed time: 36900.410736 msecs"

There is doubtless a much, much faster way to accomplish this, but I haven’t the time for a better analysis of the problem today. Besides, I was looking for an excuse to use that euler-totient function from problem 69…

Problem 70 was another problem that is nearly impossible to accomplish without spending time to understanding the math. A brute force approach will be far to computationally expensive to work, although it took me a few tries to discover this.

Rather than show the pages of failed attempts, I’ll just show the argument that brought me the correct answer. Each step after the first took roughly an hour to figure out, making this one of the more challenging problems.

1. A prime number will be relatively prime to all numbers less than it, and so primes will minimize n/phi(n).
2. However, there are no prime numbers p where p-1 is a permutation of p, so a simple prime number will not be the answer.
3. The next type of number that would give a large totient would be a product of two primes. Therefore, we need to find two prime numbers that when multiplied are slightly less than 1e7.
4. Careful! Not just any pair of prime numbers will work. The two prime numbers should also be near sqrt(1e7) so that the totient function will be maximized. Picking a low prime like 2 or 3 will result in a much lower totient function because many numbers will not be relatively prime.

My algorithm became the following:

1. Starting at sqrt(1e7), go up and down grabbing prime numbers.
2. Form pairs of these numbers. (ie, make combinations lazily)
3. Filter to only the pairs less than 1e7.
4. Filter to only pairs whose product and totient are permutations of each other.
5. Take a bunch (say, 50) of these numbers, and pick the one with the lowest n/phi(n).

In code that looks like:

;; Finally, this is how I did it:
(use '[clojure.contrib.lazy-seqs :only (primes)])

(defn prime? [n]
  (and (< 1 n)
       (not-any? #(zero? (rem n %)) (take-while #(<= (* % %) n) primes))))

(def candidates
  (let [center 3162] ;; Rough estimate of sqrt(1e7)
    (interleave (take-while #(> % 1) (filter prime? (iterate dec center)))
                (take-while #(> (* 2 center) %)
                            (filter prime? (iterate inc center))))))

(defn lazy-combinations
  "Returns lazy combinations of the list, possibly infinitely long."
  ([as] (lazy-combinations [(first as)] (rest as)))
  ([as bs]
     (let [b (first bs)]
       (lazy-cat (map #(vector % b) as)
                 (lazy-combinations (conj as b) (rest bs))))))

(defn totient [a b]
  (* a b (reduce * (map #(- 1 (/ 1 %)) [a b]))))

(defn are-perms? [a b] (= (sort (str a)) (sort (str b))))

(defn euler-70-works []
  (->> (lazy-combinations candidates)
       (filter (fn [[a b]] (and (< (* a b) 1e7)
                                (are-perms? (* a b) (totient a b)))))
       (take 50) ;; arbitrary, but it works
       (map (fn [[a b]] [(* a b) (/ (* a b) (totient a b))]))
       (sort-by second)
       first
       first))

(time (euler-70-works)) ;;"Elapsed time: 1658.035361 msecs"

Far from actually needing to take the first 50 candidates after filtering, the above algorithm will give the answer in only 4 candidates. However, during the exploration of the problem, I didn’t actually know which pair of prime numbers would work, so I was a little conservative and took more than was needed to be sure it would give the correct answer. Since it only takes a second to run, performance isn’t an issue.

Whew! That was harder than I expected!

After a long and wonderful three week vacation (and another week of just getting back into the flow), it’s time to return to the daily little bit of practice to stay sharp. Let’s get started!

Problem 69 wasn’t a programming problem at all, it was a math problem. One could solve it with a few minutes of logic, a pencil, and paper.

The key insight is that multiplying small primes together will result in products with the fewest number of relative primes. Because a small phi will maximize n/phi(n), our answer is in fact simply the product of the first few primes.

(use '[clojure.contrib.lazy-seqs :only (primes)])

(defn euler-69 []
  (last (take-while #(< % 1e6) (reductions * primes))))

(time (euler-69)) ;; "Elapsed time: 0.263869 msecs"

I really was pleased when I realized how reductions made this code a one-liner. Problem solved!

Although not really necessary, if you need some help or insight into the totient function, you might enjoy playing with these functions, as I did.

(defn coprime? [p q]
  {:pre [(< 0 p) (< 0 q)]}
  (not-any? #(and (zero? (rem p %))
                  (zero? (rem q %)))
            (take-while #(<= (* 2 %) p) primes)))

(defn brute-totient [n]
  (count (filter #(coprime? n %) (range 1 n))))

Of course, it’s much wiser to use an analytic equation for the totient of larger numbers (read up on Euler Products):

;; Prime-factors-of from problem 12
(defn prime-factors-of
  "Returns a sorted list of prime factors of num, including multiplicity."
  [num]
  (let [q (Math/sqrt num)
        factor? (fn [nom den] (zero? (rem nom den)))]
    (loop [n num
           d 2
           r []]
      (cond
       (> d q) (concat r [n])
       (= n d) (concat r [n])
       (factor? n d) (recur (/ n d) d (conj r d))
       true          (recur n (inc d) r)))))

(defn unique-prime-factors-of [n]
  (vec (into #{} (prime-factors-of n))))

(defn euler-totient [n]
  (* n (reduce * (map #(- 1 (/ 1 %)) (unique-prime-factors-of n)))))

Ciao for today!