Theory Farey_Ford

section ‹Farey Sequences and Ford Circles›

theory Farey_Ford
  imports  "HOL-Analysis.Analysis" "HOL-Number_Theory.Totient" "HOL-Library.Sublist"

begin

subsection ‹Library material›

(*added to distribution 2024-04-11*)
lemma sum_squared_le_sum_of_squares:
  fixes f :: "'a  real"
  assumes "finite I"
  shows "(iI. f i)2  (yI. (f y)2) * card I"
proof (cases "finite I  I  {}")
  case True
  then have "(iI. f i / real (card I))2  (iI. (f i)2 / real (card I))"
    using assms convex_on_sum [OF _ _ convex_power2, where a = "λx. 1 / real(card I)" and S=I]
    by simp
  then show ?thesis
    using assms  
    by (simp add: divide_simps power2_eq_square split: if_split_asm flip: sum_divide_distrib)
qed auto

(*added to distribution 2024-04-11*)
lemma sum_squared_le_sum_of_squares_2:
  "(x+y)/2  sqrt ((x2 + y2) / 2)"
proof -
  have "(x + y)2 / 2^2  (x2 + y2) / 2"
    using sum_squared_le_sum_of_squares [of UNIV "λb. if b then x else y"]
    by (simp add: UNIV_bool add.commute)
  then show ?thesis
    by (metis power_divide real_le_rsqrt)
qed

(*added to distribution 2024-04-10*)
lemma sphere_scale:
  assumes "a  0"
  shows   "(λx. a *R x) ` sphere c r = sphere (a *R c :: 'a :: real_normed_vector) (¦a¦ * r)"
proof -
  have *: "(λx. a *R x) ` sphere c r  sphere (a *R c) (¦a¦ * r)" for a r and c :: 'a
    by (metis (no_types, opaque_lifting) scaleR_right_diff_distrib dist_norm image_subsetI mem_sphere norm_scaleR)
  have "sphere (a *R c) (¦a¦ * r) = (λx. a *R x) ` (λx. inverse a *R x) ` sphere (a *R c) (¦a¦ * r)"
    unfolding image_image using assms by simp
  also have "  (λx. a *R x) ` sphere (inverse a *R (a *R c)) (¦inverse a¦ * (¦a¦ * r))"
    using "*" by blast
  also have " = (λx. a *R x) ` sphere c r"
    using assms by (simp add: algebra_simps)
  finally have "sphere (a *R c) (¦a¦ * r)  (λx. a *R x) ` sphere c r" .
  moreover have "(λx. a *R x) ` sphere c r  sphere (a *R c) (¦a¦ * r)"
    using "*" by blast
  ultimately show ?thesis by blast
qed

(*added to distribution 2024-04-10*)
lemma sphere_cscale:
  assumes "a  0"
  shows   "(λx. a * x) ` sphere c r = sphere (a * c :: complex) (cmod a * r)"
proof -
  have *: "(λx. a * x) ` sphere c r  sphere (a * c) (cmod a * r)" for a r c
    by (metis (no_types, lifting) dist_complex_def image_subsetI mem_sphere norm_mult
    right_diff_distrib')
  have "sphere (a * c) (cmod a * r) = (λx. a * x) ` (λx. inverse a * x) ` sphere (a * c) (cmod a * r)"
    by (simp add: image_image inverse_eq_divide)
  also have "  (λx. a * x) ` sphere (inverse a * (a * c)) (cmod (inverse a) * (cmod a * r))"
    using "*" by blast
  also have " = (λx. a * x) ` sphere c r"
    using assms by (simp add: field_simps flip: norm_mult)
  finally have "sphere (a * c) (cmod a * r)  (λx. a * x) ` sphere c r" .
  moreover have "(λx. a * x) ` sphere c r  sphere (a * c) (cmod a * r)"
    using "*" by blast
  ultimately show ?thesis by blast
qed

(*added to distribution 2024-05-15*)
lemma Complex_divide_complex_of_real: "Complex x y / of_real r = Complex (x/r) (y/r)"
  by (metis complex_of_real_mult_Complex divide_inverse mult.commute of_real_inverse)
lemma cmod_neg_real: "cmod (Complex (-x) y) = cmod (Complex x y)"
  by (metis complex_cnj complex_minus complex_mod_cnj norm_minus_cancel)

subsection ‹Farey sequences›

lemma sorted_two_sublist:
  fixes x:: "'a::order"
  assumes sorted: "sorted_wrt (<) l"
  shows "sublist [x, y] l  x < y  x  set l  y  set l  (z  set l. z  x  z  y)"
proof (cases "x < y  x  set l  y  set l")
  case True
  then obtain xs us where us: "l = xs @ [x] @ us"
    by (metis append_Cons append_Nil in_set_conv_decomp_first)
  with True assms have "y  set us"
    by (fastforce simp add: sorted_wrt_append)
  then obtain ys zs where yz: "l = xs @ [x] @ ys @ [y] @ zs"
    by (metis split_list us append_Cons append_Nil)
  have "sublist [x, y] l  ys = []"
    using sorted yz 
    apply (simp add: sublist_def sorted_wrt_append)
    by (metis (mono_tags, opaque_lifting) append_Cons_eq_iff append_Nil assms sorted_wrt.simps(2)
        sorted_wrt_append less_irrefl)
  also have " = (z  set l. z  x  z  y)"
    using sorted yz
    apply (simp add: sublist_def sorted_wrt_append)
    by (metis Un_iff empty_iff less_le_not_le list.exhaust list.set(1) list.set_intros(1))
  finally show ?thesis
    using True by blast 
next
  case False
  then show ?thesis
    by (metis list.set_intros(1) set_subset_Cons sorted_wrt.simps(2) sorted_wrt_append sublist_def
        set_mono_sublist sorted subset_iff)
qed

(* added to distribution 2025-03-20*)
lemma quotient_of_rat_of_int [simp]: "quotient_of (rat_of_int i) = (i, 1)"
  using Rat.of_int_def quotient_of_int by force

(* added to distribution 2025-03-20*)
lemma quotient_of_rat_of_nat [simp]: "quotient_of (rat_of_nat i) = (int i, 1)"
  by (metis of_int_of_nat_eq quotient_of_rat_of_int)

(* added to distribution 2025-03-20*)
lemma int_div_le_self: 
  x div k  x if 0 < x  for x k :: int
  by (metis div_by_1 int_div_less_self less_le_not_le nle_le nonneg1_imp_zdiv_pos_iff order.trans that)

(* not clear to do with these transp lemmas, are they of general use?*)
lemma transp_add1_int:
  assumes "n::int. R (f n) (f (1 + n))"
    and "n < n'"
    and "transp R"
  shows "R (f n) (f n')"
proof -
  have "R (f n) (f (1 + n + int k))" for k
    by (induction k) (use assms in auto elim!: transpE)
  then show ?thesis
    by (metis add.commute assms(2) zle_iff_zadd zless_imp_add1_zle)
qed

lemma refl_transp_add1_int:
  assumes "n::int. R (f n) (f (1 + n))"
      and "n  n'"
      and "reflp R" "transp R"
    shows "R (f n) (f n')"
  by (metis assms order_le_less reflpE transp_add1_int)

lemma transp_Suc:
  assumes "n. R (f n) (f (Suc n))"
      and "n < n'"
      and "transp R"
    shows "R (f n) (f n')"
proof -
  have "R (f n) (f (1 + n + k))" for k
  by (induction k) (use assms in auto elim!: transpE)
  then show ?thesis
    by (metis add_Suc_right add_Suc_shift assms(2) less_natE plus_1_eq_Suc)
qed

lemma refl_transp_Suc:
  assumes "n. R (f n) (f (Suc n))"
      and "n  n'"
      and "reflp R" "transp R"
    shows "R (f n) (f n')"
  by (metis assms dual_order.order_iff_strict reflpE transp_Suc)

(* added to distribution 2025-03-20*)
lemma sorted_subset_imp_subseq:
  fixes xs :: "'a::order list"
  assumes "set xs  set ys" "sorted_wrt (<) xs" "sorted_wrt (≤) ys"
  shows "subseq xs ys"
  using assms
proof (induction xs arbitrary: ys)
  case Nil
  then show ?case
    by auto
next
  case (Cons x xs)
  then have "x  set ys"
    by auto
  then obtain us vs where §: "ys = us @ [x] @ vs"
    by (metis append.left_neutral append_eq_Cons_conv split_list) 
  moreover 
  have "set xs  set vs"
    using Cons.prems by (fastforce simp: § sorted_wrt_append)
  with Cons have "subseq xs vs"
    by (metis § sorted_wrt.simps(2) sorted_wrt_append)
  ultimately show ?case
    by auto
qed

lemma coprime_unimodular_int:
  fixes a b::int
  assumes "coprime a b" "a>1" "b>1"
  obtains x y where "a*x - b*y = 1" "0 < x" "x < b" "0 < y" "y < a"
proof -
  obtain u v where 1: "a * u + b * v = 1"
    by (metis coprime a b cong_iff_lin coprime_iff_invertible_int)
  define k where "k  u div b"
  define x where "x  u - k*b"
  define y where "y  -(v + k*a)"
  show thesis
  proof
    show *: "a * x - b * y = 1" 
      using 1 by (simp add: x_def y_def algebra_simps)
    have "u  k*b" "b>0"
      using assms "*"  by (auto simp: k_def x_def y_def zmult_eq_neg1_iff) 
    moreover have "u - k*b  0"
      by (simp add: k_def b>0 minus_div_mult_eq_mod)
    ultimately show "x > 0"
      by (fastforce simp: x_def)
    show "x < b"
      by (simp add: 0 < b k_def minus_div_mult_eq_mod x_def)
    have "a*x > 1"
      by (metis 0 < x a>1 int_one_le_iff_zero_less less_1_mult linorder_not_less
          mult.right_neutral nle_le)
    then have "¬ b * y  0"
      using "*" by linarith
    then show "y > 0"
      by (simp add: 0 < b mult_less_0_iff order_le_less)
    show "y < a"
      using "*" 0 < x x < b
      by (smt (verit, best) a>1 mult.commute mult_less_le_imp_less)
  qed
qed

subsection ‹Farey Fractions›

type_synonym farey = rat

definition num_farey :: "farey  int" 
  where "num_farey  λx. fst (quotient_of x)" 

definition denom_farey :: "farey  int"
  where "denom_farey  λx. snd (quotient_of x)" 

definition farey :: "[int,int]  farey" 
  where "farey  λa b. max 0 (min 1 (Fract a b))"

lemma farey01 [simp]: "0  farey a b" "farey a b  1"
  by (auto simp: min_def max_def farey_def)

lemma farey_0 [simp]: "farey 0 n = 0"
  by (simp add: farey_def rat_number_collapse)

lemma farey_1 [simp]: "farey 1 1 = 1"
  by (auto simp: farey_def rat_number_expand)

lemma num_farey_nonneg: "x  {0..1}  num_farey x  0"
  by (cases x) (simp add: num_farey_def quotient_of_Fract zero_le_Fract_iff)

lemma num_farey_le_denom: "x  {0..1}  num_farey x  denom_farey x"
  by (cases x) (simp add: Fract_le_one_iff denom_farey_def num_farey_def quotient_of_Fract)

lemma denom_farey_pos: "denom_farey x > 0"
  by (simp add: denom_farey_def quotient_of_denom_pos')

lemma coprime_num_denom_farey [intro]: "coprime (num_farey x) (denom_farey x)"
  by (simp add: denom_farey_def num_farey_def quotient_of_coprime)

lemma rat_of_farey_conv_num_denom:
  "x = rat_of_int (num_farey x) / rat_of_int (denom_farey x)"
  by (simp add: denom_farey_def num_farey_def quotient_of_div)

lemma num_denom_farey_eqI:
  assumes "x = of_int a / of_int b" "b > 0" "coprime a b"
  shows   "num_farey x = a" "denom_farey x = b"
  using Fract_of_int_quotient assms quotient_of_Fract
  by (auto simp: num_farey_def denom_farey_def)

lemma farey_cases [cases type, case_names farey]:
  assumes "x  {0..1}"
  obtains a b where "0a" "ab" "coprime a b" "x = Fract a b"
  by (metis Fract_of_int_quotient Rat_cases assms num_denom_farey_eqI num_farey_le_denom num_farey_nonneg)

lemma rat_of_farey: "x = of_int a / of_int b; x  {0..1}  x = farey a b"
  by (simp add: Fract_of_int_quotient farey_def max_def)

lemma farey_num_denom_eq [simp]: "x  {0..1}  farey (num_farey x) (denom_farey x) = x"
  using rat_of_farey rat_of_farey_conv_num_denom by fastforce

lemma farey_eqI:
  assumes "num_farey x = num_farey y" "denom_farey x = denom_farey y"
  shows   "x=y"
  by (metis Rat.of_int_def assms rat_of_farey_conv_num_denom)

lemma
  assumes "coprime a b" "0a" "a<b"
  shows num_farey_eq [simp]: "num_farey (farey a b) = a"
  and denom_farey_eq [simp]: "denom_farey (farey a b) = b"
  using Fract_less_one_iff quotient_of_Fract zero_le_Fract_iff
  using assms num_denom_farey_eqI rat_of_farey by force+

lemma
  assumes "0  a" "a  b" "0<b"
  shows num_farey: "num_farey (farey a b) = a div (gcd a b)" 
    and denom_farey: "denom_farey (farey a b) = b div (gcd a b)"
proof -
  have "0  Fract a b" "Fract a b  1"
    using assms by (auto simp: Fract_le_one_iff zero_le_Fract_iff)
  with assms show "num_farey (farey a b) = a div (gcd a b)" "denom_farey (farey a b) = b div (gcd a b)"
    by (auto simp: num_farey_def denom_farey_def farey_def quotient_of_Fract Rat.normalize_def Let_def)
qed

lemma
  assumes "coprime a b" "0<b"
  shows num_farey_Fract [simp]: "num_farey (Fract a b) = a"
  and denom_farey_Fract [simp]: "denom_farey (Fract a b) = b"
  using Fract_of_int_quotient assms num_denom_farey_eqI by force+

lemma num_farey_0 [simp]: "num_farey 0 = 0"
  and denom_farey_0 [simp]: "denom_farey 0 = 1"
  and num_farey_1 [simp]: "num_farey 1 = 1"
  and denom_farey_1 [simp]: "denom_farey 1 = 1"
  by (auto simp: num_farey_def denom_farey_def)

lemma num_farey_neq_denom: "denom_farey x  1  num_farey x  denom_farey x"
  by (metis denom_farey_0 div_0 div_self num_farey_1 rat_of_farey_conv_num_denom)

lemma num_farey_0_iff [simp]: "num_farey x = 0  x = 0"
  unfolding num_farey_def
  by (metis div_0 eq_fst_iff of_int_0 quotient_of_div rat_zero_code)

lemma denom_farey_le1_cases:
  assumes "denom_farey x  1" "x  {0..1}"
  shows "x = 0  x = 1"
proof -
  consider "num_farey x = 0" | "num_farey x = 1" "denom_farey x = 1"
    using assms num_farey_le_denom [of x] num_farey_nonneg [of x] by linarith
  then show ?thesis
    by (metis denom_farey_1 farey_eqI num_farey_0_iff num_farey_1)
qed

definition mediant :: "farey  farey  farey" where 
  "mediant  λx y. Fract (fst (quotient_of x) + fst (quotient_of y)) 
                         (snd (quotient_of x) + snd (quotient_of y))"

lemma mediant_eq_Fract:
  "mediant x y = Fract (num_farey x + num_farey y) (denom_farey x + denom_farey y)"
  by (simp add: denom_farey_def num_farey_def mediant_def)

lemma mediant_eq_farey:
  assumes "x  {0..1}" "y  {0..1}"
  shows "mediant x y = farey (num_farey x + num_farey y) (denom_farey x + denom_farey y)"
proof -
  have "0  num_farey x + num_farey y"
    using assms num_farey_nonneg by auto
  moreover have "num_farey x + num_farey y  denom_farey x + denom_farey y"
    by (meson add_mono assms num_farey_le_denom)
  ultimately show ?thesis
    by (simp add: add_pos_pos denom_farey_pos Fract_of_int_quotient rat_of_farey mediant_eq_Fract)
qed


definition farey_unimodular :: "farey  farey  bool" where
  "farey_unimodular x y 
     denom_farey x * num_farey y - num_farey x * denom_farey y = 1"

lemma farey_unimodular_imp_less:
  assumes "farey_unimodular x y"
  shows   "x < y"
  using assms
  by (auto simp: farey_unimodular_def rat_less_code denom_farey_def num_farey_def)

lemma denom_mediant: "denom_farey (mediant x y)  denom_farey x + denom_farey y"
  using quotient_of_denom_pos' [of x] quotient_of_denom_pos' [of y]
  by (simp add: mediant_eq_Fract denom_farey_def num_farey_def quotient_of_Fract normalize_def Let_def int_div_le_self)

lemma unimodular_imp_both_coprime:
  fixes a:: "'a::{algebraic_semidom,comm_ring_1}"
  assumes "b*c - a*d = 1" 
  shows   "coprime a b" "coprime c d"
  using mult.commute by (metis assms coprimeI dvd_diff dvd_mult2)+

lemma unimodular_imp_coprime:
  fixes a:: "'a::{algebraic_semidom,comm_ring_1}"
  assumes "b*c - a*d = 1" 
  shows   "coprime (a+c) (b+d)"
proof (rule coprimeI)
  fix k 
  assume k: "k dvd (a+c)" "k dvd (b+d)"
  moreover have "(b+d)*c = (a+c)*d + 1"
    using assms by (simp add: algebra_simps)
  ultimately show "is_unit k"
    by (metis add_diff_cancel_left' dvd_diff dvd_mult2)
qed

definition fareys :: "int  rat list"
  where "fareys n  sorted_list_of_set {x  {0..1}. denom_farey x  n}"

lemma strict_sorted_fareys: "sorted_wrt (<) (fareys n)"
  by (simp add: fareys_def)

lemma farey_set_UN_farey: "{x  {0..1}. denom_farey x  n} = (b  {1..n}. a  {0..b}. {farey a b})"
proof -
  have "b  {1..n}. a  {0..b}. x = farey a b"
    if "denom_farey x  n" "x  {0..1}" for x :: farey
    unfolding Bex_def
    using that denom_farey_pos int_one_le_iff_zero_less num_farey_le_denom num_farey_nonneg
    by fastforce
  moreover have "b a::int. 0  b  b div gcd a b  b"
    by (metis div_0 int_div_le_self nless_le)
  ultimately show ?thesis
    by (auto simp: denom_farey) (meson order_trans)
qed

lemma farey_set_UN_farey': "{x  {0..1}. denom_farey x  n} = (b  {1..n}. a  {0..b}. if coprime a b then {farey a b} else {})"
proof -
  have "aa ba. farey a b = farey aa ba  0  aa  aa  ba  1  ba  ba  n  coprime aa ba"
    if "1  b" and "b  n" and "0  a" and "a  b" for a b
  proof -
    let ?a = "a div gcd a b"
    let ?b = "b div gcd a b"
    have "coprime ?a ?b"
      by (metis div_gcd_coprime not_one_le_zero b1)
    moreover have "farey a b = farey ?a ?b"
      using Fract_coprime farey_def by presburger
    moreover have "?a  ?b  ?b  n"
      by (smt (verit, best) gcd_pos_int int_div_le_self that zdiv_mono1)
    ultimately show ?thesis
      using that by (metis denom_farey denom_farey_pos div_int_pos_iff gcd_ge_0_int int_one_le_iff_zero_less)
  qed
  then show ?thesis
    unfolding farey_set_UN_farey
    by (fastforce split: if_splits)
qed

lemma farey_set_UN_Fract: "{x  {0..1}. denom_farey x  n} = (b  {1..n}. a  {0..b}. {Fract a b})"
  unfolding farey_set_UN_farey
  by (simp add: Fract_of_int_quotient farey_def)

lemma farey_set_UN_Fract': "{x  {0..1}. denom_farey x  n} = (b  {1..n}. a  {0..b}. if coprime a b then {Fract a b} else {})"
  unfolding farey_set_UN_farey'
  by (simp add: Fract_of_int_quotient farey_def)

lemma finite_farey_set: "finite {x  {0..1}. denom_farey x  n}"
  unfolding farey_set_UN_farey by blast

lemma denom_in_fareys_iff: "x  set (fareys n)  denom_farey x  int n  x  {0..1}"
  using finite_farey_set by (auto simp: fareys_def)

lemma denom_fareys_leI: "x  set (fareys n)  denom_farey x  n"
  using finite_farey_set by (auto simp: fareys_def)

lemma denom_fareys_leD: "denom_farey x  int n; x  {0..1}  x  set (fareys n)"
  using denom_in_fareys_iff by blast

lemma fareys_increasing_1: "set (fareys n)  set (fareys (1 + n))"
  using farey_set_UN_farey by (force simp: fareys_def)

definition fareys_new :: "int  rat set" where
  "fareys_new n  {Fract a n| a. coprime a n  a  {0..n}}"

lemma fareys_new_0 [simp]: "fareys_new 0 = {}"
  by (auto simp: fareys_new_def)

lemma fareys_new_1 [simp]: "fareys_new 1 = {0,1}"
proof -
  have "Fract a 1 = 1"
    if a: "Fract a 1  0" "0  a" "a  1" for a
    by (metis One_rat_def int_one_le_iff_zero_less nless_le order_antisym_conv
        rat_number_collapse(1) that) 
  moreover have "a. Fract a 1 = 0  0  a  a  1"
    using rat_number_expand(1) by auto
  moreover have "a. Fract a 1 = 1  0  a  a  1"
    using One_rat_def by fastforce
  ultimately show ?thesis
    by (auto simp: fareys_new_def)
qed

lemma fareys_new_not01:
  assumes "n>1"
  shows "0  (fareys_new n)" "1  (fareys_new n)"
  using assms by (simp_all add: Fract_of_int_quotient fareys_new_def)

lemma inj_num_farey: "inj_on num_farey (fareys_new n)"
proof (cases "n=1")
  case True
  then show ?thesis
    using fareys_new_1 by auto
next
  case False
  then show ?thesis
  proof -
    have "Fract a n = Fract a' n"
      if "coprime a n" "0  a" "a  n"
        and "coprime a' n" "0  a'" "a'  n"
        and "num_farey (Fract a n) = num_farey (Fract a' n)"
      for a a'
    proof -
      from that
      obtain "a < n" "a' < n"
        using False by force+
      with that show ?thesis
        by auto
    qed
    with False show ?thesis
      by (auto simp: inj_on_def fareys_new_def)
  qed
qed
 
lemma finite_fareys_new [simp]: "finite (fareys_new n)"
  by (auto simp: fareys_new_def)

lemma card_fareys_new:
  assumes "n > 1"
  shows "card (fareys_new (int n)) = totient n"
proof -
  have "bij_betw num_farey (fareys_new n) (int ` totatives n)"
  proof -
    have "a>0. a  n  coprime a n  num_farey x = int a"
      if x: "x  fareys_new n" for x
    proof -
      obtain a where a: "x = Fract a (int n)" "coprime a (int n)" "0  a" "a  int n"
        using x by (auto simp: fareys_new_def)
      then have "a > 0"
        using assms less_le_not_le by fastforce
      moreover have "coprime (nat a) n"
        by (metis a(2,3) coprime_int_iff nat_0_le)
      ultimately have "num_farey x = int (nat a)"
        using a num_farey by auto
      then show ?thesis
        using 0 < a coprime (nat a) n a(4) nat_le_iff zero_less_nat_eq by blast
    qed
    moreover have "xfareys_new n. int a = num_farey x"
      if "0 < a" "a  n" "coprime a n" for a 
    proof -
      have §: "coprime (int a) (int n)" "0  (int a)" "(int a)  int n"
        using that by auto
      then have "Fract (int a) (int n) = Fract (int a) (int n)"
        using Fract_of_int_quotient assms rat_of_farey by auto
      with § have "Fract (int a) (int n)  fareys_new n"
        by (auto simp: fareys_new_def)
      then have "int a = num_farey (Fract (int a) n)"
        using coprime (int a) (int n) assms by auto
      then show ?thesis
        using Fract (int a) (int n)  fareys_new (int n) by blast
    qed
    ultimately show ?thesis
      by (auto simp add: totatives_def bij_betw_def inj_num_farey comp_inj_on image_iff)
  qed
  then show ?thesis
    unfolding totient_def by (metis bij_betw_same_card bij_betw_of_nat)
qed

lemma disjoint_fareys_plus1: 
  assumes "n > 0"
  shows "disjnt (set (fareys n)) (fareys_new (1 + n))"
proof -
  have False
    if §: "0  a" "a  1 + n" "coprime a (1 + n)"
      "1  d" "d  n" "Fract a (1 + n) = Fract c d" "0  c" "c  d" "coprime c d"
    for a c d
  proof (cases "c<d")
    case True
    have alen: "a  n"
      using nle_le that by fastforce
    have "d = denom_farey (Fract c d)"
      using that by force
    also have " = 1 + n"
      using denom_farey_Fract that by fastforce
    finally show False
      using that(5) by fastforce
  next
    case False
    with c  d have "c=d" by auto
    with that have "d=1" by force
    with that have "1 + n = 1"
      by (metis One_rat_def c = d assms denom_farey_1 denom_farey_Fract le_imp_0_less
          order_less_le)
    then show ?thesis
      using c = d § assms by auto
  qed
  then show ?thesis
  unfolding fareys_def farey_set_UN_Fract' fareys_new_def disjnt_iff
    by auto
qed


lemma set_fareys_plus1: "set (fareys (1 + n)) = set (fareys n)  fareys_new (1 + n)"
proof -
  have "b1. b  n  (a0. a  b  coprime a b  Fract c d = Fract a b)"
    if "Fract c d  fareys_new (1 + n)"
      and "coprime c d" "1  d" "d  1 + n" "0  c" "c  d"
    for c d
  proof (cases "d = 1 + n")
    case True
    with that show ?thesis
      by (auto simp: fareys_new_def)
  qed (use that in auto)
  moreover have "d1. d  1 + n  (c0. c  d  coprime c d  x = Fract c d)"
    if "x  fareys_new (1 + n)" for x
    using that nle_le by (fastforce simp add: fareys_new_def)
  ultimately show ?thesis
    unfolding fareys_def farey_set_UN_Fract' by fastforce
qed

lemma length_fareys_Suc: 
  assumes "n>0"
  shows "length (fareys (1 + int n)) = length (fareys n) + totient (Suc n)"
proof -
  have "length (fareys (1 + int n)) = card (set (fareys (1 + int n)))"
    by (metis fareys_def finite_farey_set sorted_list_of_set.sorted_key_list_of_set_unique)
  also have " = card (set (fareys n)) + card (fareys_new(1 + int n))"
    using disjoint_fareys_plus1 assms by (simp add: set_fareys_plus1 card_Un_disjnt)
  also have " = card (set (fareys n)) + totient (Suc n)"
    using assms card_fareys_new by force
  also have " = length (fareys n) + totient (Suc n)"
    using fareys_def finite_farey_set by auto
  finally show ?thesis .
qed


lemma fareys_0 [simp]: "fareys 0 = []"
  unfolding fareys_def farey_set_UN_farey
  by simp

lemma fareys_1 [simp]: "fareys 1 = [0, 1]"
proof -
  have "{x  {0..1}. denom_farey x  1} = {0,1}"
    using denom_farey_le1_cases by auto
  then show ?thesis
    by (simp add: fareys_def)
qed

lemma fareys_2 [simp]: "fareys 2 = [0, farey 1 2, 1]"
proof -
  have §: "denom_farey x  2  denom_farey x = 1  denom_farey x = 2" for x
    using denom_farey_pos [of x] by auto
  have "{x  {0..1}. denom_farey x  2} = {farey 0 1, farey 1 2, farey 1 1}"
  proof -
    have "x = farey 1 1"
      if "x  farey 0 1" "x  {0..1}" "denom_farey x = 1" for x
      using that denom_farey_le1_cases order.eq_iff rat_of_farey by auto
    moreover have False
      if "x  farey 0 1" "x  farey 1 2" "denom_farey x = 2" "x  {0..1}" for x
      using that num_farey_neq_denom
      by (smt (verit) farey_0 farey_num_denom_eq num_farey_le_denom num_farey_nonneg)
    moreover have "denom_farey (farey 1 1) = 1"
       by (simp add: Fract_of_int_quotient farey_def)
    ultimately show ?thesis
      by (auto simp: farey_set_UN_farey §)
  qed
  also have " = {0, 1/2, 1::rat}"
    by (simp add: farey_def Fract_of_int_quotient)
  finally show ?thesis
    by (simp add: fareys_def farey_def Fract_of_int_quotient)
qed

lemma length_fareys_1: 
  shows "length (fareys 1) = 1 + totient 1"
  by simp

lemma length_fareys: "n>0  length (fareys n) = 1 + (k=1..n. totient k)"
proof (induction n)
  case (Suc n)
  then show ?case
    by (cases "n=0") (auto simp add: length_fareys_Suc)
qed auto

lemma subseq_fareys_1: "subseq (fareys n) (fareys (1 + n))"
  by (metis fareys_increasing_1 strict_sorted_fareys sorted_subset_imp_subseq strict_sorted_imp_sorted)

lemma monotone_fareys: "monotone (≤) subseq fareys"
  using refl_transp_add1_int [OF _ _ subseq_order.reflp_on_le subseq_order.transp_on_le]
  by (metis monotoneI subseq_fareys_1)

lemma farey_unimodular_0_1 [simp, intro]: "farey_unimodular 0 1"
  by (auto simp: farey_unimodular_def)

text ‹Apostol's Theorem 5.2 for integers›
lemma mediant_lies_betw_int:
  fixes a b c d::int
  assumes "rat_of_int a / of_int b < of_int c / of_int d" "b>0" "d>0"
  shows "rat_of_int a / of_int b < (of_int a + of_int c) / (of_int b + of_int d)" 
        "(rat_of_int a + of_int c) / (of_int b + of_int d) < of_int c / of_int d"
    using assms by (simp_all add: field_split_simps)

text ‹Apostol's Theorem 5.2›
theorem mediant_inbetween:
  fixes x y::farey
  assumes "x < y"
  shows "x < mediant x y" "mediant x y < y"
  using assms mediant_lies_betw_int Fract_of_int_quotient
  by (metis denom_farey_pos mediant_eq_Fract of_int_add rat_of_farey_conv_num_denom)+

lemma sublist_fareysD:
  assumes "sublist [x,y] (fareys n)"
  obtains "x  set (fareys n)" "y  set (fareys n)"
  by (meson assms list.set_intros set_mono_sublist subsetD)

text ‹Adding the denominators of two consecutive Farey fractions›
lemma sublist_fareys_add_denoms:
  fixes a b c d::int
  defines "x  Fract a b"
  defines "y  Fract c d"
  assumes sub: "sublist [x,y] (fareys (int n))" and "b>0" "d>0" "coprime a b" "coprime c d"
  shows "b + d > n"
proof (rule ccontr)
  have §: "x < y" "z  set (fareys n). z  x  z  y"
    using sorted_two_sublist strict_sorted_fareys sub by blast+
  assume "¬ int n < b + d"
  with assms have "denom_farey (mediant x y)  int n"
    by (metis denom_farey_Fract denom_mediant dual_order.trans leI)
  then have "mediant x y  set (fareys n)"
    by (metis sub atLeastAtMost_iff denom_in_fareys_iff farey01 mediant_eq_farey
        sublist_fareysD)
  moreover have "x < mediant x y" "mediant x y < y"
    by (simp_all add: mediant_inbetween x < y)
  ultimately show False
    using § x_def y_def by fastforce
qed

subsection ‹Apostol's Theorems 5.3--5.5›

theorem consec_subset_fareys:
  fixes a b c d::int
  assumes abcd: "0  Fract a b" "Fract a b < Fract c d" "Fract c d  1"
    and consec: "b*c - a*d = 1"
    and max: "max b d  n" "n < b+d" 
    and "b>0" 
  shows "sublist [Fract a b, Fract c d] (fareys n)"
proof (rule ccontr)
  assume con: "¬ ?thesis"
  have "d > 0"
    using max by force
  have "coprime a b" "coprime c d"
    using consec unimodular_imp_both_coprime by blast+
  with b > 0 d > 0 have "denom_farey (Fract a b) = b" "denom_farey (Fract c d) = d"
    by auto
  moreover have "bn" "dn"
    using max by auto
  ultimately have ab: "Fract a b  set (fareys n)" and cd: "Fract c d  set (fareys n)"
    using abcd finite_farey_set by (auto simp: fareys_def)
  then obtain xs us where us: "fareys n = xs @ [Fract a b] @ us"
    using abcd by (metis append_Cons append_Nil split_list)
  have "Fract c d  set us"
    using abcd cd strict_sorted_fareys [of n]
    by (fastforce simp add: us sorted_wrt_append)
  then obtain ys zs where yz: "fareys n = xs @ [Fract a b] @ ys @ [Fract c d] @ zs"
    by (metis split_list us append_Cons append_Nil)
  with con have "ys  []"
    by (metis Cons_eq_append_conv sublist_appendI)
  then obtain h k where hk: "coprime h k" "Fract h k  set ys"  "k>0"
    by (metis Rat_cases list.set_sel(1))
  then have hk_fareys: "Fract h k  set (fareys n)" 
    by (auto simp: yz)
  have less: "Fract a b < Fract h k" "Fract h k < Fract c d" 
    using hk strict_sorted_fareys [of n] by (auto simp add: yz sorted_wrt_append)
  with b > 0 d > 0 hk have *: "k*a < h*b" "d*h < c*k"
    by (simp_all add: Fract_of_int_quotient mult.commute divide_simps flip: of_int_mult)
  have "k  n"
    using hk by (metis hk_fareys denom_fareys_leI denom_farey_Fract)
  have "k = (b*c - a*d)*k"
    by (simp add: consec)
  also have " = b*(c*k - d*h) + d*(b*h - a*k)"
    by (simp add: algebra_simps)
  finally have k: "k = b * (c * k - d * h) + d * (b * h - a * k)" .  
  moreover have "c*k - d*h  1" "b*h - a*k  1"
    using b > 0 d > 0 * by (auto simp: mult.commute)
  ultimately have "b * (c * k - d * h) + d * (b * h - a * k)  b+d"
    by (metis b > 0 d > 0 add_mono mult.right_neutral mult_left_mono
        order_le_less)
  then show False
    using k  n max k by force
qed

lemma farey_unimodular_mediant:
  assumes "farey_unimodular x y"
  shows "farey_unimodular x (mediant x y)" "farey_unimodular (mediant x y) y"
  using assms quotient_of_denom_pos' [of x] quotient_of_denom_pos' [of y]
  unfolding farey_unimodular_def
  by (auto simp: mediant_eq_Fract denom_farey_def num_farey_def quotient_of_Fract unimodular_imp_coprime algebra_simps)

text ‹Apostol's Theorem 5.4›
theorem mediant_unimodular:
  fixes a b c d::int
  assumes abcd: "0  Fract a b" "Fract a b < Fract c d" "Fract c d  1"
    and consec: "b*c - a*d = 1"
    and 0: "b>0" "d>0" 
  defines "h  a+c"
  defines "k  b+d"
  obtains "Fract a b < Fract h k" "Fract h k < Fract c d" "coprime h k"
          "b*h - a*k = 1"  "c*k - d*h = 1"
proof
  show "Fract a b < Fract h k" "Fract h k < Fract c d"
    using abcd 0
    by (simp_all add: Fract_of_int_quotient h_def k_def distrib_left distrib_right divide_simps)
  show "coprime h k"
    by (simp add: consec unimodular_imp_coprime h_def k_def)
  show "b * h - a * k = 1"
    by (simp add: consec distrib_left h_def k_def)
  show "c * k - d * h = 1"
    by (simp add: consec h_def distrib_left k_def mult.commute)
qed

text ‹Apostol's Theorem 5.5, first part: "Each fraction in @{term"F(n+1)"} which is not in @{term"F(n)"}
      is the mediant of a pair of consecutive fractions in @{term"F(n)"}

lemma get_consecutive_parents:
  fixes m n::int
  assumes "coprime m n" "0<m" "m<n"
  obtains a b c d where "m = a+c" "n = b+d" "b*c - a*d = 1" "a0" "b>0" "c>0" "d>0" "a<b" "cd"
proof (cases "m=1")
  case True
  show ?thesis
  proof
    show "m = 0 + 1" "n = 1 + (n-1)"
      by (auto simp: True)
  qed (use True m<n in auto)
next
  case False
  then obtain d c where *: "n*c - m*d = 1" "0 < d" "d < n" "0 < c" "c < m"
    using coprime_unimodular_int
 [of n m] coprime_commute assms by (smt (verit) coprime_commute)
  then have **: "n * (c - d) + (n - m) * d = 1"
    by (metis mult_diff_mult)
  show ?thesis
  proof
    show "cd"
      using * ** m<n by (smt (verit) mult_le_0_iff)
    show "(n-d) * c - (m-c) * d = 1"
      using * by (simp add: algebra_simps)
    with * m<n show "m-c < n-d"
      by (smt (verit, best) mult_mono)
  qed (use * in auto)
qed

theorem fareys_new_eq_mediant:
  assumes "x  fareys_new n" "n>1"
  obtains a b c d where 
    "sublist [Fract a b, Fract c d] (fareys (n-1))" 
    "x = mediant (Fract a b) (Fract c d)" "coprime a b" "coprime c d" "a0" "b>0" "c>0" "d>0" 
proof -
  obtain m where m: "coprime m n" "0  m" "m  n" "x = Fract m n"
    using assms nless_le zero_less_imp_eq_int by (force simp: fareys_new_def)
  moreover
  have "x  0" "x  1"
    using assms fareys_new_not01 by auto
  with m have "0<m" "m<n"
    using n>1 of_nat_le_0_iff by fastforce+
  ultimately
  obtain a b c d where 
    abcd: "m = a+c" "n = b+d" "b*c - a*d = 1" "a0" "b>0" "c>0" "d>0" "a<b" "cd"
    by (metis get_consecutive_parents)
  show thesis
  proof
    have "Fract a b < Fract c d"
      using abcd mult.commute[of b c] by force
    with consec_subset_fareys
    show "sublist [Fract a b, Fract c d] (fareys (n-1))"
      using Fract_le_one_iff abcd zero_le_Fract_iff by auto
    show "x = mediant (Fract a b) (Fract c d)"
      using abcd x = Fract m n mediant_eq_Fract unimodular_imp_both_coprime by fastforce
    show "coprime a b" "coprime c d"
      using b * c - a * d = 1 unimodular_imp_both_coprime by blast+
  qed (use abcd in auto)
qed

text ‹Apostol's Theorem 5.5, second part: "Moreover, if @{term"a/b<c/d"} are consecutive in any @{term"F(n)"},
then they satisfy the unimodular relation @{term"bc - ad = 1"}.›
theorem consec_imp_unimodular:
  assumes "sublist [Fract a b, Fract c d] (fareys (int n))" "b>0" "d>0" "coprime a b" "coprime c d"
  shows  "b*c - a*d = 1"
  using assms 
proof (induction n arbitrary: a b c d)
  case 0
  then show ?case
    by auto
next
  case (Suc n)
  show ?case
  proof (cases "n=0")
    case True
    with Suc.prems have "Fract a b = 0" "Fract c d = 1"
      by (auto simp add: sublist_Cons_right)
    with Suc.prems obtain "a=0" "b=1" "c=1" "d=1"
      by (auto simp: Fract_of_int_quotient)
    then show ?thesis
      by auto
  next
    case False
    have "Fract a b < Fract c d" 
      and ab: "Fract a b  set (fareys (1 + int n))" and cd: "Fract c d  set (fareys (1 + int n))"
      using strict_sorted_fareys [of "Suc n"] Suc.prems
      by (auto simp add: sublist_def sorted_wrt_append)
    have con: "z  Fract a b  Fract c d  z" if "z  set (fareys (int n))" for z
    proof -
      have "z  set (fareys (1 + int n))"
        using fareys_increasing_1 that by blast
      with Suc.prems strict_sorted_fareys [of "1 + int n"] show ?thesis
        by (fastforce simp add: sublist_def sorted_wrt_append)
    qed
    show ?thesis
    proof (cases "Fract a b  set (fareys n)  Fract c d  set (fareys n)")
      case True
      then have "sublist [Fract a b, Fract c d] (fareys n)"
        using con Fract a b < Fract c d sorted_two_sublist strict_sorted_fareys by blast
      then show ?thesis
        by (simp add: Suc) 
    next
      case False
      have notboth: False if §: "Fract a b  fareys_new (1+n)" "Fract c d  fareys_new (1+n)"
      proof -
        obtain a' b' c' d' where eq':
          "sublist [Fract a' b', Fract c' d'] (fareys n)" "Fract a b = mediant (Fract a' b') (Fract c' d')"
          by (smt (verit) n0 § fareys_new_eq_mediant of_nat_Suc of_nat_le_0_iff plus_1_eq_Suc)
        then have abcd': "Fract a' b'  set (fareys n)" "Fract c' d'  set (fareys n)"
          by (auto simp: sublist_def)
        have con': "z  Fract a' b'  Fract c' d'  z" if "z  set (fareys n)" for z
          by (metis eq'(1) sorted_two_sublist strict_sorted_fareys that)
        have "Fract a' b' < Fract c' d'"
          using eq'(1) sorted_two_sublist [OF strict_sorted_fareys] by blast         
        then obtain A: "Fract a' b' < Fract a b" "Fract a b < Fract c' d'"
          using eq'(2) mediant_inbetween by presburger
          obtain a'' b'' c'' d'' where eq'':
          "sublist [Fract a'' b'', Fract c'' d''] (fareys n)" "Fract c d = mediant (Fract a'' b'') (Fract c'' d'')"
          by (smt (verit) § n0 fareys_new_eq_mediant of_nat_Suc of_nat_le_0_iff plus_1_eq_Suc)
        then have abcd'': "Fract a'' b''  set (fareys n)" "Fract c'' d''  set (fareys n)"
          by (auto simp: sublist_def)
        then have "Fract c'' d''  set (fareys (1 + int n))"
          using fareys_increasing_1 by blast
        have con'': "z  Fract a'' b''  Fract c'' d''  z" if "z  set (fareys n)" for z
          using sorted_two_sublist [OF strict_sorted_fareys] eq''(1) that by blast
        then obtain "Fract a'' b'' < Fract c'' d''" "Fract a'' b'' < Fract c d" "Fract c d < Fract c'' d''"
          by (metis eq'' sorted_two_sublist [OF strict_sorted_fareys] mediant_inbetween)
        with A show False
          using con' con'' abcd' abcd'' con Fract a b < Fract c d
          by (metis eq'(2) eq''(2) dual_order.strict_trans1 not_less_iff_gr_or_eq)
      qed
      consider "Fract a b  fareys_new (1+n)" | "Fract c d  fareys_new (1+n)"
        using False set_fareys_plus1 [of n]
        by (metis (mono_tags, opaque_lifting) Suc.prems(1) Suc_eq_plus1_left UnE list.set_intros(1) of_nat_Suc set_mono_sublist
            set_subset_Cons subset_iff)
      then show ?thesis
      proof cases
        case 1
        then obtain a' b' c' d' where eq:
          "sublist [Fract a' b', Fract c' d'] (fareys n)" 
          "Fract a b = mediant (Fract a' b') (Fract c' d')" "coprime a' b'" "coprime c' d'" "b'>0" "d'>0" 
          by (smt (verit) n0 fareys_new_eq_mediant of_nat_Suc of_nat_le_0_iff plus_1_eq_Suc)
        then have abcd': "Fract a' b'  set (fareys n)" "Fract c' d'  set (fareys n)"
          by (auto simp: sublist_def)
        have con': "z  Fract a' b'  Fract c' d'  z" if "z  set (fareys n)" for z
          using eq(1) sorted_two_sublist strict_sorted_fareys that by blast
        obtain "Fract a' b' < Fract c' d'" "Fract a b < Fract c' d'"
          using eq by (simp add: mediant_inbetween(2) sorted_two_sublist strict_sorted_fareys)
        then have "Fract c d  Fract c' d'"
          using con abcd' linorder_not_less by blast
        moreover have "Fract c' d'  Fract c d" 
          if "Fract c d  set (fareys n)"
          by (metis con' Fract a b < Fract c d Fract a' b' < Fract c' d' eq(2) order.trans linorder_not_less mediant_inbetween(1)
              nless_le that)
        ultimately have "Fract c' d' = Fract c d"
          using notboth "1" cd set_fareys_plus1 by auto
        with Suc.prems obtain "c' = c" "d' = d"
          by (metis 0 < d' coprime c' d' denom_farey_Fract num_farey_Fract)
        then have uni: "b'*c - a'*d = 1"
          using Suc eq by blast
        then obtain "a = a' + c" "b = b' + d"
          using eq Suc.prems apply (simp add: mediant_eq_Fract)
          by (metis c' = c d' = d denom_farey_Fract num_farey_Fract pos_add_strict
              unimodular_imp_coprime)
        with uni show ?thesis
          by (auto simp: algebra_simps)
      next
        case 2
        then obtain a' b' c' d' where eq:
          "sublist [Fract a' b', Fract c' d'] (fareys n)" 
          "Fract c d = mediant (Fract a' b') (Fract c' d')" "coprime a' b'" "coprime c' d'" "b'>0" "d'>0" 
          by (smt (verit) n0 fareys_new_eq_mediant of_nat_Suc of_nat_le_0_iff plus_1_eq_Suc)
        then have abcd': "Fract a' b'  set (fareys n)" "Fract c' d'  set (fareys n)"
          by (auto simp: sublist_def)
        have con': "z  Fract a' b'  Fract c' d'  z" if "z  set (fareys n)" for z
          using eq(1) sorted_two_sublist strict_sorted_fareys that by blast
        obtain "Fract a' b' < Fract c' d'" "Fract a' b' < Fract c d"
          using eq mediant_inbetween
          by (metis sorted_two_sublist strict_sorted_fareys)
        then have "Fract a' b'  Fract a b"
          using con abcd' linorder_not_less by blast
        moreover have "Fract a b  Fract a' b'" 
          if "Fract a b  set (fareys n)"
          by (metis Fract a b < Fract c d Fract a' b' < Fract c' d' con' order.strict_trans2 eq(2) mediant_inbetween(2)
              not_less_iff_gr_or_eq that)
        ultimately have "Fract a' b' = Fract a b"
          using notboth 2 ab set_fareys_plus1 by auto
        with Suc.prems obtain "a' = a" "b' = b"
          by (metis 0 < b' coprime a' b' denom_farey_Fract num_farey_Fract)
        then have uni: "b*c' - a*d' = 1"
          using Suc.IH Suc.prems eq by blast
        then obtain "c = a + c'" "d = b + d'"
          using eq Suc.prems apply (simp add: mediant_eq_Fract)
          by (metis a' = a b' = b denom_farey_Fract num_farey_Fract pos_add_strict
              unimodular_imp_coprime)
        with uni show ?thesis
          by (auto simp: algebra_simps)
      qed
    qed
  qed
qed

subsection ‹Ford circles›

definition Ford_center :: "rat  complex" where
  "Ford_center r  (λ(h,k). Complex (h/k) (1/(2 * k^2))) (quotient_of r)"

definition Ford_radius :: "rat  real" where
  "Ford_radius r  (λ(h,k). 1/(2 * k^2)) (quotient_of r)"

definition Ford_tan :: "[rat,rat]  bool" where
  "Ford_tan r s  dist (Ford_center r) (Ford_center s) = Ford_radius r + Ford_radius s"

definition Ford_circle :: "rat  complex set" where
  "Ford_circle r  sphere (Ford_center r) (Ford_radius r)"

lemma Im_Ford_center [simp]: "Im (Ford_center r) = Ford_radius r"
  by (auto simp: Ford_center_def Ford_radius_def split: prod.splits)

lemma Ford_radius_nonneg: "Ford_radius r  0"
  by (simp add: Ford_radius_def split: prod.splits)

lemma two_Ford_tangent:
  assumes r: "(a,b) = quotient_of r" and s: "(c,d) = quotient_of s"
  shows "(dist (Ford_center r) (Ford_center s))^2 - (Ford_radius r + Ford_radius s)^2 
       = ((a*d - b*c)^2 - 1) / (b*d)^2"
proof -
  obtain 0: "b > 0" "d > 0"
    by (metis assms quotient_of_denom_pos)
  have 1: "dist (Ford_center r) (Ford_center s) ^ 2 = (a/b - c/d)^2 + (1/(2*b^2) - 1/(2*d^2)) ^ 2"
    using assms by (force simp: Ford_center_def dist_norm complex_norm complex_diff split: prod.splits)
  have 2: "(Ford_radius r + Ford_radius s) ^ 2 = (1/(2*b^2) + 1/(2*d^2)) ^ 2"
    using assms by (force simp: Ford_radius_def split: prod.splits)
  show ?thesis
    using 0 unfolding 1 2 by (simp add: field_simps eval_nat_numeral)
qed

text ‹Apostol's Theorem 5.6›
lemma two_Ford_tangent_iff:
  assumes r: "(a,b) = quotient_of r" and s: "(c,d) = quotient_of s"
  shows "Ford_tan r s  ¦b * c - a * d¦ = 1"
proof -
  obtain 0: "b > 0" "d > 0"
    by (metis assms quotient_of_denom_pos)
  have "Ford_tan r s  dist (Ford_center r) (Ford_center s) ^ 2 = (Ford_radius r + Ford_radius s) ^ 2"
    using Ford_radius_nonneg by (simp add: Ford_tan_def)
  also have "  ((a*d - b*c)^2 - 1) / (b*d)^2 = 0"
    using two_Ford_tangent [OF assms] by (simp add: diff_eq_eq)
  also have "  ¦b * c - a * d¦ = 1"
    using 0 by (simp add: abs_square_eq_1 abs_minus_commute flip: of_int_mult of_int_diff)
  finally show ?thesis .
qed

text ‹Also Apostol's Theorem 5.6: Distinct Ford circles do not overlap›
lemma Ford_no_overlap:
  assumes "r  s"
  shows "dist (Ford_center r) (Ford_center s)  Ford_radius r + Ford_radius s"
proof -
  obtain a b c d where r: "(a,b) = quotient_of r" and s: "(c,d) = quotient_of s"
                 and "b>0" "d>0"
    by (metis quotient_of_denom_pos surj_pair)
  moreover have "a  c  b  d"
    using assms r s quotient_of_inject by force
  ultimately have "a * d  c * b"
    by (metis Fract_of_int_quotient assms eq_rat(1) less_irrefl quotient_of_div)
  then have "(a*d - b*c)^2  (1::int)"
    by (simp add: mult.commute int_one_le_iff_zero_less)
  then have "((a*d - b*c)^2 - 1) / (b*d)^2  (0::real)"
    by (simp add: divide_simps mult_less_0_iff flip: of_int_mult of_int_power)
  then show ?thesis
    using two_Ford_tangent [OF r s]
    by (metis (no_types, lifting) ge_iff_diff_ge_0 of_int_1 of_int_diff of_int_mult 
              of_int_power power2_le_imp_le zero_le_dist)
qed

lemma Ford_aux1:
  assumes "a0"
  shows "cmod (Complex (b / (a * (a2 + b2))) (1 / (2 * a2) - inverse (a2 + b2))) = 1 / (2 * a2)"
       (is "cmod ?z = ?r")
proof -
  have "(2 * a2) * cmod ?z = cmod ((2 * a2) * ?z)"
    by (simp add: norm_mult power2_eq_square)
  also have " = cmod (Complex (2*a*b / (a2 + b2)) (1 - (2 * a2) / (a2 + b2)))"
    unfolding complex_of_real_mult_Complex inverse_eq_divide
    using a0 by (simp add: power2_eq_square mult.assoc right_diff_distrib)
  also have " = cmod (Complex (2*a*b) ((a2 + b2) - (2 * a2)) / (a2 + b2))"
    unfolding Complex_divide_complex_of_real diff_divide_distrib
    using assms by force
  also have " = cmod (Complex (2*a*b) ((a2 + b2) - (2 * a2))) / (a2 + b2)"
    by (smt (verit) norm_divide norm_of_real not_sum_power2_lt_zero)
  also have " = sqrt ((a2 + b2)^2) / (a2 + b2)"
    unfolding power2_eq_square complex_norm
    by (simp add: algebra_simps)
  also have " = 1"
    using assms by auto
  finally show ?thesis
    by (metis inverse_eq_divide inverse_unique)
qed

lemma Ford_aux2:
  assumes "a0"
  shows "cmod (Complex (a / (b * (b2 + a2)) - 1 / (a * b)) (1 / (2 * a2) - inverse (b2 + a2))) = 1 / (2 * a2)"
       (is "cmod ?z = ?r")
proof -
  have "a / (b * (b2 + a2)) - 1 / (a * b) = -b / (a * (b2 + a2))"
    by (simp add: divide_simps power2_eq_square)
  then have "cmod ?z = cmod (Complex (b / (a * (a2 + b2))) (1 / (2 * a2) - inverse (a2 + b2)))"
    by (simp add: cmod_neg_real add.commute)
  also have " = 1 / (2 * a2)"
    using Ford_aux1 assms by simp
  finally show ?thesis .
qed

definition Radem_trans :: "rat  complex  complex" where
  "Radem_trans  λr τ. let (h,k) = quotient_of r in - 𝗂 * of_int k ^ 2 * (τ - of_rat r)"

text ‹Theorem 5.8 first part›
lemma Radem_trans_image: "Radem_trans r ` Ford_circle r = sphere (of_rat (1/2)) (1/2)"
proof -
  obtain h k where r: "quotient_of r = (h,k)" and "k>0" and req: "r = of_int h / of_int k"
    using quotient_of_denom_pos quotient_of_div by fastforce
  have "Radem_trans r ` Ford_circle r = ((*) (-𝗂 * of_int k ^ 2)) ` (λτ. τ - of_rat r) ` Ford_circle r"
    by (simp add: Radem_trans_def r image_comp)
  also have " = ((*) (-𝗂 * of_int k ^ 2)) ` sphere (Ford_center r - of_rat r) (Ford_radius r)"
    by (simp add: Ford_circle_def flip: sphere_translation_subtract)
  also have " = sphere (- 𝗂 * (of_int k)2 * (Ford_center r - r))
                          (cmod (- 𝗂 * (of_int k)2) * Ford_radius r)"
    using k>0 by (intro sphere_cscale) auto
  also have " = sphere (of_rat (1/2)) (1/2)"
  proof -
    have "(- 𝗂 * (of_int k)2 * (Ford_center r - r)) = 1/2"
      using k>0
      apply (simp add: Ford_center_def r algebra_simps Complex_eq)
      by (simp add: of_rat_divide req)
    moreover 
    have "(cmod (- 𝗂 * (of_int k)2) * Ford_radius r) = 1/2"
      using k>0
      by (simp add: norm_mult norm_power Ford_radius_def r)
    ultimately show ?thesis
      by (metis of_rat_1 of_rat_divide of_rat_numeral_eq)
  qed
  finally show ?thesis .
qed

locale three_Ford =
  fixes N::nat
  fixes h1 k1 h k h2 k2::int
  assumes sub1: "sublist [Fract h1 k1, Fract h k] (fareys (int N))"
  assumes sub2: "sublist [Fract h k, Fract h2 k2] (fareys (int N))"
  assumes coprime: "coprime h1 k1" "coprime h k" "coprime h2 k2"
  assumes k_pos: "k1 > 0" "k > 0" "k2 > 0"

begin

definition "r1  Fract h1 k1"
definition "r  Fract h k"
definition "r2  Fract h2 k2"

lemma N_pos: "N>0"
  using sub1 gr0I by force

lemma r_eq_quotient:
  "(h1,k1) = quotient_of r1" "(h,k) = quotient_of r" "(h2,k2) = quotient_of r2"
  by (simp_all add: coprime k_pos quotient_of_Fract r1_def r_def r2_def)

lemma r_eq_divide:
  "r1 = of_int h1 / of_int k1" "r = of_int h / of_int k" "r2 = of_int h2/ of_int k2"
  by (simp_all add: Fract_of_int_quotient of_rat_divide r1_def r2_def r_def)

lemma collapse_r:
  "real_of_int h1 / of_int k1 = of_rat r1" 
  "real_of_int h / of_int k = of_rat r" "real_of_int h2/ of_int k2 = of_rat r2"
  by (simp_all add: of_rat_divide r_eq_divide)

lemma unimod1: "k1*h - h1*k = 1"
  and unimod2: "k*h2 - h*k2 = 1"
  using consec_imp_unimodular coprime k_pos sub1 sub2 by blast+

lemma r_less: "r1 < r" "r < r2"
  using r1_def r_def r2_def sub1 sub2 sorted_two_sublist [OF strict_sorted_fareys] by auto

lemma r01:
  obtains "r1  {0..1}" "r  {0..1}" "r2  {0..1}"
  by (metis denom_in_fareys_iff r1_def r2_def r_def sub1 sub2 sublist_fareysD)

lemma atMost_N:
  obtains "k1  N" "k  N" "k2  N"
  by (metis denom_farey_def denom_in_fareys_iff prod.sel(2) r1_def r2_def r_def r_eq_quotient sub1 sub2
      sublist_fareysD)

lemma greaterN1: "k1 + k > N"
  using sublist_fareys_add_denoms coprime k_pos sub1 by blast

lemma greaterN2: "k + k2 > N"
  using sublist_fareys_add_denoms coprime k_pos sub2 by blast

definition "alpha1  Complex (h/k - k1 / of_int(k * (k2 + k12))) (inverse (of_int (k2 + k12)))"
definition "alpha2  Complex (h/k + k2 / of_int(k * (k2 + k22))) (inverse (of_int (k2 + k22)))"

definition "zed1  Complex (k2) (k*k1) / ((k2 + k12))"
definition "zed2  Complex (k2) (- k*k2) / ((k2 + k22))"

text ‹Apostol's Theorem 5.7›
lemma three_Ford_tangent:
  obtains "alpha1  Ford_circle r" "alpha1  Ford_circle r1"
          "alpha2  Ford_circle r" "alpha2  Ford_circle r2"
proof
  show "alpha1  Ford_circle r"
    using k_pos Ford_aux1 r_eq_quotient
    by (force simp: alpha1_def Ford_circle_def Ford_center_def dist_norm complex_diff 
        Ford_radius_def split: prod.splits)
  have 1: "real_of_int h1 / real_of_int k1 = real_of_int h / real_of_int k - 1 / (k1*k)"
    using unimod1 k_pos
    by (simp add: divide_simps) (simp add: algebra_simps flip: of_int_mult of_int_diff)
  show "alpha1  Ford_circle r1"
    using k_pos Ford_aux2 r_eq_quotient
    by (force simp: alpha1_def Ford_circle_def Ford_center_def dist_norm complex_diff 1 Ford_radius_def split: prod.splits)
  show "alpha2  Ford_circle r"
    using k_pos Ford_aux1 [of k k2] cmod_neg_real r_eq_quotient
    by (force simp add: alpha2_def Ford_circle_def Ford_center_def dist_norm complex_diff Ford_radius_def split: prod.splits)
  have 2: "real_of_int h / real_of_int k = real_of_int h2 / real_of_int k2 - 1 / (k*k2)"
    using unimod2 k_pos
    by (simp add: divide_simps) (simp add: algebra_simps flip: of_int_mult of_int_diff)
  show "alpha2  Ford_circle r2"
    using k_pos Ford_aux2 [of k2 k] cmod_neg_real r_eq_quotient
    apply (simp add: alpha2_def Ford_circle_def Ford_center_def dist_norm complex_diff 2 Ford_radius_def split: prod.splits)
    by (smt (verit) mult.commute prod.sel)
qed

text ‹Theorem 5.8 second part, for alpha1›
lemma Radem_trans_alpha1: "Radem_trans r alpha1 = zed1"
proof -
  have "Radem_trans r alpha1 = ((*) (-𝗂 * of_int k ^ 2)) ((λτ. τ - of_rat r) alpha1)"
    by (metis Radem_trans_def prod.simps(2) r_eq_quotient(2))
  also have " = ((*) (-𝗂 * of_int k ^ 2)) (Complex (- k1 / of_int(k * (k2 + k12))) (inverse (of_int (k2 + k12))))"
    using k_pos by (simp add: alpha1_def r_def of_rat_rat Complex_eq)
  also have " = zed1"
    unfolding complex_eq_iff by (simp add: zed1_def inverse_eq_divide power2_eq_square)
  finally show ?thesis .
qed

text ‹Theorem 5.8 second part, for alpha2›
lemma Radem_trans_alpha2: "Radem_trans r alpha2 = zed2"
proof -
  have "Radem_trans r alpha2 = ((*) (-𝗂 * of_int k ^ 2)) ((λτ. τ - of_rat r) alpha2)"
    by (metis Radem_trans_def prod.simps(2) r_eq_quotient(2))
  also have " = ((*) (-𝗂 * of_int k ^ 2)) (Complex (k2 / of_int(k * (k2 + k22))) (inverse (of_int (k2 + k22))))"
    using k_pos by (simp add: alpha2_def r_def of_rat_rat Complex_eq)
  also have " = zed2"
    unfolding complex_eq_iff by (simp add: zed2_def inverse_eq_divide power2_eq_square)
  finally show ?thesis .
qed

text ‹Theorem 5.9, for zed1›
lemma cmod_zed1: "cmod zed1 = k / sqrt (k2 + k12)"
proof -
  have "cmod zed1 ^ 2 = (k^4 + k2 * k12) / (k2 + k12)^2"
    by (simp add: zed1_def cmod_def divide_simps)
  also have " = (of_int k) ^ 2 / (k2 + k12)"
    by (simp add: eval_nat_numeral divide_simps) argo
  finally have "cmod zed1 ^ 2 = (of_int k) ^ 2 / (k2 + k12)" .
  with k_pos real_sqrt_divide show ?thesis
    unfolding cmod_def by force
qed

text ‹Theorem 5.9, for zed2›
lemma cmod_zed2: "cmod zed2 = k / sqrt (k2 + k22)"
proof -
  have "cmod zed2 ^ 2 = (k^4 + k2 * k22) / (k2 + k22)^2"
    by (simp add: zed2_def cmod_def divide_simps)
  also have " = (of_int k) ^ 2 / (k2 + k22)"
    by (simp add: eval_nat_numeral divide_simps) argo
  finally have "cmod zed2 ^ 2 = (of_int k) ^ 2 / (k2 + k22)" .
  with k_pos real_sqrt_divide show ?thesis
    unfolding cmod_def by force
qed


text ‹The last part of theorem 5.9›

lemma RMS_calc:
  assumes "k' > 0" "k' + k > int N"
  shows "k / sqrt (k2 + k'2) < sqrt 2 * k/N"
proof -
  have §: "(k + k')/2  sqrt ((k2 + k'2) / 2)"
    using sum_squared_le_sum_of_squares_2 by simp
  have "N / sqrt 2 < (N+1) / sqrt 2"
    by (simp add: divide_strict_right_mono)
  also have "  (k + k') / sqrt 2"
    using assms by (simp add: divide_simps)
  also have " = (k + k')/2 * sqrt 2"
    by (metis nonzero_divide_eq_eq real_div_sqrt times_divide_eq_right zero_le_numeral
        zero_neq_numeral)
  also have "  sqrt (k2 + k'2)"
    using § by (simp add: le_divide_eq real_sqrt_divide)
  finally have 1: "real N / sqrt 2 < sqrt (real_of_int (k2 + k'2))" .
  with N_pos k_pos not_sum_power2_lt_zero show ?thesis
    by (force simp add: cmod_zed1 mult.commute divide_simps)
qed

lemma on_chord_bounded_cmod:
  assumes "z  closed_segment zed1 zed2"
  shows "cmod z < sqrt 2 * k/N"
proof -
  have "cmod z  max (cmod zed1) (cmod zed2)"
    using segment_furthest_le [OF assms, of 0] by auto
  moreover
  obtain "cmod zed1 < sqrt 2 * k/N" "cmod zed2 < sqrt 2 * k/N"
    using RMS_calc cmod_zed1 cmod_zed2 greaterN1 greaterN2 k_pos by force
  ultimately show ?thesis
    using assms cmod_zed1 cmod_zed2 by linarith
qed

end

end