Problem 71 asks us to find a proper fraction just slightly less than 3/7, with a denominator less than one million.

Step 1, as usual, was to just code something brute-force to prove by example that it won’t work.

;; Idea #1: maybe I can exploit the HCF(n,d) = 1 property (coprime) 

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

(defn gcd [a b] (if (zero? b) a (recur b (mod a b)))) 

(defn reduced-proper-fractions [num]
  (for [d (range 1 (inc num))
        n (range 1 d)
        :when (= 1 (gcd n d))]
    (/ n d)))

(time (doall (sort (reduced-proper-fractions 1000)))) 

;; Performance is O(n^2)... much too slow for 1,000,000 elements!

Then I got a bit smarter, and did some reading. Apparently, a sorted list of rational fractions is called a Farey Sequence, and is closely related to a Stern-Brocot tree. More interesting to me was the geometric interpretation of the Farey Sequence, which are called Ford Circles.

My second idea was to play around with the Farey Sequences a little and see if they would be a suitable way of solving the problem.

;; Idea #2: Generate the reduced proper factors in order

(defn next-farey-element [[n a b c d]]
  (let [k (int (/ (+ n b) d))
        p (- (* k c) a)
        q (- (* k d) b)]
    [n c d p q]))

(defn make-farey-seq [n]
  (take-while #(< % 1)
              (map (fn [[n a b c d]] (/ a b))
                   (iterate next-farey-element [n 0 1 1 n]))))

(defn make-decending-farey-seq [n]
  (take-while #(> % 0)
              (map (fn [[n a b c d]] (/ a b))
                   (iterate next-farey-element [n 1 1 (dec n) n]))))

(make-farey-seq 8)
(make-decending-farey-seq 8)
</span>

Although I’m sure there is a way to efficiently solve this problem using a Farey Sequence, it was unclear to me how to start a sequence at anywhere but the top (1) or bottom (0).

Changing strategies, I thought about the continued fractions I had learned about earlier in problem 64 and 65. Perhaps there was a way to solve the problem by using very similar continued fractions?

;; Idea #3: Use the properties of continued fractions

(defn eval-continued-frac [s]
  (reduce #(+ %2 (/ 1 %1)) (last s) (reverse (butlast s))))

(defn frac-as-continued 
  "Returns a sequence of the coefficicents of the continuing fraction
  that is equivilent to ratio frac."
  [frac]
  {:pre [(rational? frac)]}
  (->> (rest (iterate (fn [[x i]]
                   (if (= i (int i))
                     [i 0]
                     [(int i) (/ 1 (- i (int i)))]))
                      [0 frac]))
       (take-while (fn [[x i]] (not (and (zero? x) (zero? i)))))
       (map first)))

;; Remembering that (eval-continued-frac [0 2 2 1]))  -> 3/7
;; So the entry just to the left of it will be almost the same!
;; This will look for a continued fraction slightly less than 3/7
(defn euler-71-weird []
  (->> (map #(eval-continued-frac [0 2 2 1 %]) (iterate inc 1000))
       (take-while #(> 1000000 (denominator %)))
       last))

(time (euler-71-weird)) ;; "Elapsed time: 1442.474092 msecs"

Aha! I did it! Actually, there was quite a bit of other code I wrote to test that would actually work, but let’s ignore that and just consider the above, which concentrates on the core idea: fractions close to each other will have very similar continued fraction expansions. Chopping off digits at the end of the expansion makes the fraction bigger, and adding digits at the end makes them smaller.

If you don’t understand why that worked, I really suggest playing around with eval-continued-frac and frac-as-continued until you intuitively understand the continued fractions a bit better. It will also help to remember that rational fractions have two representations as continued fractions. For example, 3/7 can be written [0 2 3] or [0 2 2 1], which may be helpful.

Addendum: After reading the page at Project Euler regarding how others solved the problem, I smacked my head in disbelief at how stupid I had been to ignore the obvious algorithm. If I had it to do over again, I would have….ah well, let’s do it over again right now!

;; Idea #4: We want the number with the largest numerator n such that
;;                               n/d < 3/7
;;          So we find the biggest n where:    n < (3/7)*d           

(defn euler-71-analysis []
  (loop [bn 0
         bd 1
         d 2]
    (let [n (int (/ (- (* 3 d) 1) 7))]
      (if (>= d 1000000)
        (/ bn bd)
        (if (< (/ bn bd) (/ n d))
          (recur n d (inc d))
          (recur bn bd (inc d)))))))

(time (euler-71-analysis)) ;; "Elapsed time: 739.874679 msecs"

Clear, concise, efficient…why didn’t I think of that the first time!?

Still, I think the continued fractions technique worked surprisingly well, considering I was fumbling around in the dark with just my intuition.