# HG changeset patch
# User Steve Losh <steve@stevelosh.com>
# Date 1541115866 14400
# Node ID 5e5dd3fd5ae005c36c39e5d0cb96cf0579e886e2
# Parent  99fa06464617ebc7d2cad78be1ba63fb4907f86c
Problems 71 and 73

diff -r 99fa06464617 -r 5e5dd3fd5ae0 src/problems.lisp
--- a/src/problems.lisp	Tue Jul 31 16:18:33 2018 +0000
+++ b/src/problems.lisp	Thu Nov 01 19:44:26 2018 -0400
@@ -1592,6 +1592,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:
@@ -2365,6 +2419,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 99fa06464617 -r 5e5dd3fd5ae0 src/utils.lisp
--- a/src/utils.lisp	Tue Jul 31 16:18:33 2018 +0000
+++ b/src/utils.lisp	Thu Nov 01 19:44:26 2018 -0400
@@ -167,6 +167,30 @@
                          (cl-strings:split l ,delimiter%))))
            (,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
@@ -345,6 +369,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))
@@ -781,8 +815,6 @@
 (define-modify-macro adjoinf (item &rest keyword-args) adjoin%)
 
 
-
-
 (defun pythagorean-triplet-p (a b c)
   (math a ^ 2 + b ^ 2 = c ^ 2))
 
@@ -818,6 +850,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."
@@ -974,7 +1040,6 @@
   (* precision (round number precision)))
 
 
-
 ;;;; A* Search ----------------------------------------------------------------
 (defstruct path
   state