# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1545001352 18000
# Node ID 2a837c783fe744d3bd3cd2e563cfc11761398310
# Parent  5e5dd3fd5ae005c36c39e5d0cb96cf0579e886e2# Parent  c19da8761e5752072bfbf856d95e855c6e938f51
Merge.

diff -r c19da8761e57 -r 2a837c783fe7 src/problems.lisp
--- a/src/problems.lisp	Sat Dec 15 17:25:12 2018 -0500
+++ b/src/problems.lisp	Sun Dec 16 18:02:32 2018 -0500
@@ -1503,6 +1503,60 @@
     (while (<= n 1000000))
     (finally (return result))))
 
+(defun problem-71 ()
+  ;; Consider the fraction, n/d, where n and d are positive integers. If n<d and
+  ;; HCF(n,d)=1, it is called a reduced proper fraction.
+  ;;
+  ;; If we list the set of reduced proper fractions for d ≤ 8 in ascending order
+  ;; of size, we get:
+  ;;
+  ;; 1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3,
+  ;; 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
+  ;;
+  ;; It can be seen that 2/5 is the fraction immediately to the left of 3/7.
+  ;;
+  ;; By listing the set of reduced proper fractions for d ≤ 1,000,000 in
+  ;; ascending order of size, find the numerator of the fraction immediately to
+  ;; the left of 3/7.
+  (iterate
+    ;; We could theoretically iterate through F_1000000 til we find 3/7, but
+    ;; we'd have about 130000000000 elements to churn though.  We can do it much
+    ;; faster.
+    ;;
+    ;; Instead, we start by computing F₇ and finding 3/7 and whatever's to the
+    ;; left of it.  We can find what's in between these two in the next
+    ;; iteration by taking their mediant, then just repeat until we hit the
+    ;; limit.
+    (with limit = 1000000)
+    (with target = 3/7)
+    (with guess = (iterate (for x :in-farey-sequence (denominator target))
+                           (for px :previous x)
+                           (finding px :such-that (= x target))))
+    (for next = (mediant guess target))
+    (if (> (denominator next) limit)
+      (return (numerator guess))
+      (setf guess next))))
+
+(defun problem-73 ()
+  ;; Consider the fraction, n/d, where n and d are positive integers. If n<d and
+  ;; HCF(n,d)=1, it is called a reduced proper fraction.
+  ;;
+  ;; If we list the set of reduced proper fractions for d ≤ 8 in ascending order
+  ;; of size, we get:
+  ;;
+  ;; 1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3,
+  ;; 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
+  ;;
+  ;; It can be seen that there are 3 fractions between 1/3 and 1/2.
+  ;;
+  ;; How many fractions lie between 1/3 and 1/2 in the sorted set of reduced
+  ;; proper fractions for d ≤ 12,000?
+  (iterate
+    ;; Just brute force this one with our fancy Farey sequence iterator.
+    (for x :in-farey-sequence 12000)
+    (while (< x 1/2))
+    (counting (< 1/3 x))))
+
 (defun problem-74 ()
   ;; The number 145 is well known for the property that the sum of the factorial
   ;; of its digits is equal to 145:
@@ -2270,6 +2324,8 @@
 (test p63 (is (= 49 (problem-63))))
 (test p67 (is (= 7273 (problem-67))))
 (test p69 (is (= 510510 (problem-69))))
+(test p71 (is (= 428570 (problem-71))))
+(test p73 (is (= 7295372 (problem-73))))
 (test p74 (is (= 402 (problem-74))))
 (test p79 (is (= 73162890 (problem-79))))
 (test p81 (is (= 427337 (problem-81))))
diff -r c19da8761e57 -r 2a837c783fe7 src/utils.lisp
--- a/src/utils.lisp	Sat Dec 15 17:25:12 2018 -0500
+++ b/src/utils.lisp	Sun Dec 16 18:02:32 2018 -0500
@@ -181,6 +181,30 @@
                          (str:split ,delimiter% l))))
            (,kwd ,var :next (parse-line (next-line))))))))
 
+(defmacro-driver (FOR var IN-FAREY-SEQUENCE n)
+  "Iterate `var` through the Farey sequence Fₙ in ascending order.
+
+  The Farey sequence Fₙ is the sorted sequence of all reduced fractions between
+  0 and 1 (inclusive) with denominators ≤ `n`.
+
+  Example:
+
+    (iterate (for x :in-farey-sequence 5) (collect x))
+    ; =>
+    (0 1/5 1/4 1/3 2/5 1/2 3/5 2/3 3/4 4/5 1)
+
+  "
+  (with-gensyms (n% prev)
+    (let ((kwd (if generate 'generate 'for)))
+      `(progn
+         (with ,n% = ,n)
+         (with ,prev = 0)
+         (,kwd ,var :do-next
+          (if-first-time
+            (setf ,var 0)
+            (psetf ,var (or (next-farey-term ,prev ,var ,n%) (terminate))
+                   ,prev ,var)))))))
+
 
 (defmacro when-first-time (&body body)
   `(if-first-time
@@ -359,6 +383,16 @@
          (u (- 1 v w)))
     (values u v w)))
 
+(defun mediant (x y)
+  "Return the mediant of `x` and `y`."
+  ;; The mediant operation is:
+  ;;
+  ;; A   C   A + C
+  ;; - ⊕ - = -----
+  ;; B   D   B + D
+  (/ (+ (numerator x) (numerator y))
+     (+ (denominator x) (denominator y))))
+
 
 ;;;; Digits -------------------------------------------------------------------
 (defun digits-length (n &optional (radix 10))
@@ -830,6 +864,40 @@
         (recur (mv* d triple))))))
 
 
+;;;; Farey Sequences ----------------------------------------------------------
+(defun next-farey-term (x y n)
+  "Return the next term in the Farey sequence Fₙ after `x` and `y`, or `nil`.
+
+  For convenience, if `y` is zero `x` is not examined.
+
+  If `y` is `1` the sequence is complete and `nil` is returned.
+
+  "
+  (declare (type (rational 0 1) x y)
+           (type (integer 1) n))
+  (case y
+    (0 (/ 1 n))
+    (1 nil)
+    (t (let* ((a (numerator x))
+              (b (denominator x))
+              (c (numerator y))
+              (d (denominator y))
+              (k (floor (+ n b) d))
+              (p (- (* k c) a))
+              (q (- (* k d) b)))
+         (/ p q)))))
+
+(defun farey (n)
+  "Return the Farey sequence Fₙ as a fresh list."
+  (iterate (for x :in-farey-sequence n)
+           (collect x)))
+
+(defun approximate-farey-length (n)
+  "Return an approximation of the length of the Farey sequence Fₙ."
+  (/ (* 3 (expt n 2))
+     (expt pi 2)))
+
+
 ;;;; Geometric Numbers --------------------------------------------------------
 (defun squarep (n)
   "Return whether `n` is a perfect square."