Theory Continued_Fractions.Quadratic_Irrationals

(*
  File:     Continued_Fractions/Quadratic Irrationals.thy
  Author:   Manuel Eberl, University of Innsbruck
*)
section ‹Quadratic Irrationals›
theory Quadratic_Irrationals
imports
  Continued_Fractions
  "HOL-Computational_Algebra.Computational_Algebra"
  "HOL-Library.Discrete_Functions"
  "Coinductive.Coinductive_Stream"
begin

lemma snth_cycle:
  assumes "xs ≠ []"
  shows   "snth (cycle xs) n = xs ! (n mod length xs)"
proof (induction n rule: less_induct)
  case (less n)
  have "snth (shift xs (cycle xs)) n = xs ! (n mod length xs)"
  proof (cases "n < length xs")
    case True
    thus ?thesis
      by (subst shift_snth_less) auto
  next
    case False
    have "0 < length xs"
      using assms by simp
    also have "… ≤ n"
      using False by simp
    finally have "n > 0" .

    from False have "snth (shift xs (cycle xs)) n = snth (cycle xs) (n - length xs)"
      by (subst shift_snth_ge) auto
    also have "… = xs ! ((n - length xs) mod length xs)"
      using assms ‹n > 0› by (intro less) auto
    also have "(n - length xs) mod length xs = n mod length xs"
      using False by (simp add: mod_if)
    finally show ?thesis .
  qed
  also have "shift xs (cycle xs) = cycle xs"
    by (rule cycle_decomp [symmetric]) fact
  finally show ?case .
qed

subsection ‹Basic results on rationality of square roots›

lemma inverse_in_Rats_iff [simp]: "inverse (x :: real) ∈ ℚ ⟷ x ∈ ℚ"
  by (auto simp: inverse_eq_divide divide_in_Rats_iff1)

lemma nonneg_sqrt_nat_or_irrat:
  assumes "x ^ 2 = real a" and "x ≥ 0"
  shows   "x ∈ ℕ ∨ x ∉ ℚ"
proof safe
  assume "x ∉ ℕ" and "x ∈ ℚ"
  from Rats_abs_nat_div_natE[OF this(2)]
    obtain p q :: nat where q_nz [simp]: "q ≠ 0" and "abs x = p / q" and coprime: "coprime p q" .
  with ‹x ≥ 0› have x: "x = p / q"
      by simp
  with assms have "real (q ^ 2) * real a = real (p ^ 2)"
    by (simp add: field_simps)
  also have "real (q ^ 2) * real a = real (q ^ 2 * a)"
    by simp
  finally have "p ^ 2 = q ^ 2 * a"
    by (subst (asm) of_nat_eq_iff) auto
  hence "q ^ 2 dvd p ^ 2"
    by simp
  hence "q dvd p"
    by simp
  with coprime have "q = 1"
    by auto
  with x and ‹x ∉ ℕ› show False
    by simp
qed

text ‹
  A square root of a natural number is either an integer or irrational.
›
corollary sqrt_nat_or_irrat:
  assumes "x ^ 2 = real a"
  shows   "x ∈ ℤ ∨ x ∉ ℚ"
proof (cases "x ≥ 0")
  case True
  with nonneg_sqrt_nat_or_irrat[OF assms this]
    show ?thesis by (auto simp: Nats_altdef2)
next
  case False
  from assms have "(-x) ^ 2 = real a"
    by simp
  moreover from False have "-x ≥ 0"
    by simp
  ultimately have "-x ∈ ℕ ∨ -x ∉ ℚ"
    by (rule nonneg_sqrt_nat_or_irrat)
  thus ?thesis
    by (auto simp: Nats_altdef2 minus_in_Ints_iff)
qed

corollary sqrt_nat_or_irrat':
  "sqrt (real a) ∈ ℕ ∨ sqrt (real a) ∉ ℚ"
  using nonneg_sqrt_nat_or_irrat[of "sqrt a" a] by auto

text ‹
  The square root of a natural number ‹n› is again a natural number iff ‹n is a perfect square.›
›
corollary sqrt_nat_iff_is_square:
  "sqrt (real n) ∈ ℕ ⟷ is_square n"
proof
  assume "sqrt (real n) ∈ ℕ"
  then obtain k where "sqrt (real n) = real k" by (auto elim!: Nats_cases)
  hence "sqrt (real n) ^ 2 = real (k ^ 2)" by (simp only: of_nat_power)
  also have "sqrt (real n) ^ 2 = real n" by simp
  finally have "n = k ^ 2" by (simp only: of_nat_eq_iff)
  thus "is_square n" by blast
qed (auto elim!: is_nth_powerE)

corollary irrat_sqrt_nonsquare: "¬is_square n ⟹ sqrt (real n) ∉ ℚ"
  using sqrt_nat_or_irrat'[of n] by (auto simp: sqrt_nat_iff_is_square)

lemma sqrt_of_nat_in_Rats_iff: "sqrt (real n) ∈ ℚ ⟷ is_square n"
  using irrat_sqrt_nonsquare[of n] sqrt_nat_iff_is_square[of n] Nats_subset_Rats by blast

lemma Discrete_sqrt_altdef: "floor_sqrt n = nat ⌊sqrt n⌋"
proof -
  have "real (floor_sqrt n ^ 2) ≤ sqrt n ^ 2"
    by simp
  hence "floor_sqrt n ≤ sqrt n"
    unfolding of_nat_power by (rule power2_le_imp_le) auto
  moreover have "real (Suc (floor_sqrt n) ^ 2) > real n"
    unfolding of_nat_less_iff by (rule Suc_floor_sqrt_power2_gt)
  hence "real (floor_sqrt n + 1) ^ 2 > sqrt n ^ 2"
    unfolding of_nat_power by simp
  hence "real (floor_sqrt n + 1) > sqrt n"
    by (rule power2_less_imp_less) auto
  hence "floor_sqrt n + 1 > sqrt n" by simp
  ultimately show ?thesis by linarith
qed


subsection ‹Definition of quadratic irrationals›

text ‹
  Irrational real numbers $x$ that satisfy a quadratic equation $a x^2 + b x + c = 0$
  with ‹a›, ‹b›, ‹c› not all equal to 0 are called ∗‹quadratic irrationals›. These are of the form
  $p + q \sqrt{d}$ for rational numbers ‹p›, ‹q› and a positive integer ‹d›.
›
inductive quadratic_irrational :: "real ⇒ bool" where
  "x ∉ ℚ ⟹ real_of_int a * x ^ 2 + real_of_int b * x + real_of_int c = 0 ⟹
     a ≠ 0 ∨ b ≠ 0 ∨ c ≠ 0 ⟹ quadratic_irrational x"

lemma quadratic_irrational_sqrt [intro]:
  assumes "¬is_square n"
  shows   "quadratic_irrational (sqrt (real n))"
  using irrat_sqrt_nonsquare[OF assms]
  by (intro quadratic_irrational.intros[of "sqrt n" 1 0 "-int n"]) auto

lemma quadratic_irrational_uminus [intro]:
  assumes "quadratic_irrational x"
  shows   "quadratic_irrational (-x)"
  using assms
proof induction
  case (1 x a b c)
  thus ?case by (intro quadratic_irrational.intros[of "-x" a "-b" c]) auto
qed

lemma quadratic_irrational_uminus_iff [simp]:
  "quadratic_irrational (-x) ⟷ quadratic_irrational x"
  using quadratic_irrational_uminus[of x] quadratic_irrational_uminus[of "-x"] by auto

lemma quadratic_irrational_plus_int [intro]:
  assumes "quadratic_irrational x"
  shows   "quadratic_irrational (x + of_int n)"
  using assms
proof induction
  case (1 x a b c)
  define x' where "x' = x + of_int n"
  define a' b' c' where
     "a' = a" and "b' = b - 2 * of_int n * a" and
     "c' = a * of_int n ^ 2 - b * of_int n + c"
  from 1 have "0 = a * (x' - of_int n) ^ 2 + b * (x' - of_int n) + c"
    by (simp add: x'_def)
  also have "… = a' * x' ^ 2 + b' * x' + c'"
    by (simp add: algebra_simps a'_def b'_def c'_def power2_eq_square)
  finally have "… = 0" ..
  moreover have "x' ∉ ℚ"
    using 1 by (auto simp: x'_def add_in_Rats_iff2)
  moreover have "a' ≠ 0 ∨ b' ≠ 0 ∨ c' ≠ 0"
    using 1 by (auto simp: a'_def b'_def c'_def)
  ultimately show ?case
    by (intro quadratic_irrational.intros[of "x + of_int n" a' b' c']) (auto simp: x'_def)
qed

lemma quadratic_irrational_plus_int_iff [simp]:
  "quadratic_irrational (x + of_int n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int[of x n]
        quadratic_irrational_plus_int[of "x + of_int n" "-n"] by auto

lemma quadratic_irrational_minus_int_iff [simp]:
  "quadratic_irrational (x - of_int n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "-n"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_plus_nat_iff [simp]:
  "quadratic_irrational (x + of_nat n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "int n"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_minus_nat_iff [simp]:
  "quadratic_irrational (x - of_nat n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "-int n"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_plus_1_iff [simp]:
  "quadratic_irrational (x + 1) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x 1]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_minus_1_iff [simp]:
  "quadratic_irrational (x - 1) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "-1"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_plus_numeral_iff [simp]:
  "quadratic_irrational (x + numeral n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "numeral n"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_minus_numeral_iff [simp]:
  "quadratic_irrational (x - numeral n) ⟷ quadratic_irrational x"
  using quadratic_irrational_plus_int_iff[of x "-numeral n"]
  by (simp del: quadratic_irrational_plus_int_iff)

lemma quadratic_irrational_inverse:
  assumes "quadratic_irrational x"
  shows   "quadratic_irrational (inverse x)"
  using assms
proof induction
  case (1 x a b c)
  from 1 have "x ≠ 0" by auto
  have "0 = (real_of_int a * x2 + real_of_int b * x + real_of_int c) / x ^ 2"
    by (subst 1) simp
  also have "… = real_of_int c * (inverse x) ^ 2 + real_of_int b * inverse x + real_of_int a"
    using ‹x ≠ 0› by (simp add: field_simps power2_eq_square)
  finally have "… = 0" ..
  thus ?case using 1
    by (intro quadratic_irrational.intros[of "inverse x" c b a]) auto
qed

lemma quadratic_irrational_inverse_iff [simp]:
  "quadratic_irrational (inverse x) ⟷ quadratic_irrational x"
  using quadratic_irrational_inverse[of x] quadratic_irrational_inverse[of "inverse x"]
  by (cases "x = 0") auto

lemma quadratic_irrational_cfrac_remainder_iff:
  "quadratic_irrational (cfrac_remainder c n) ⟷ quadratic_irrational (cfrac_lim c)"
proof (cases "cfrac_length c = ∞")
  case False
  thus ?thesis
    by (auto simp: quadratic_irrational.simps)
next
  case [simp]: True
  show ?thesis
  proof (induction n)
    case (Suc n)
    from Suc.prems have "cfrac_remainder c (Suc n) =
                           inverse (cfrac_remainder c n - of_int (cfrac_nth c n))"
      by (subst cfrac_remainder_Suc) (auto simp: field_simps)
    also have "quadratic_irrational … ⟷ quadratic_irrational (cfrac_remainder c n)"
      by simp
    also have "… ⟷ quadratic_irrational (cfrac_lim c)"
      by (rule Suc.IH)
    finally show ?case .
  qed auto
qed
                
subsection ‹Real solutions of quadratic equations›

text ‹
  For the next result, we need some basic properties of real solutions to quadratic equations.
›
lemma quadratic_equation_reals:
  fixes a b c :: real
  defines "f ≡ (λx. a * x ^ 2 + b * x + c)"
  defines "discr ≡ (b^2 - 4 * a * c)"
  shows   "{x. f x = 0} =
             (if a = 0 then
                (if b = 0 then if c = 0 then UNIV else {} else {-c/b})
              else if discr ≥ 0 then {(-b + sqrt discr) / (2 * a), (-b - sqrt discr) / (2 * a)}
                                else {})" (is ?th1)                
proof (cases "a = 0")
  case [simp]: True
  show ?th1
  proof (cases "b = 0")
    case [simp]: True
    hence "{x. f x = 0} = (if c = 0 then UNIV else {})"
      by (auto simp: f_def)
    thus ?th1 by simp
  next
    case False
    hence "{x. f x = 0} = {-c / b}" by (auto simp: f_def field_simps)
    thus ?th1 using False by simp
  qed
next
  case [simp]: False
  show ?th1
  proof (cases "discr ≥ 0")
    case True
    {
      fix x :: real
      have "f x = a * (x - (-b + sqrt discr) / (2 * a)) * (x - (-b - sqrt discr) / (2 * a))"
        using True by (simp add: f_def field_simps discr_def power2_eq_square)
      also have "… = 0 ⟷ x ∈ {(-b + sqrt discr) / (2 * a), (-b - sqrt discr) / (2 * a)}"
        by simp
      finally have "f x = 0 ⟷ …" .
    }
    hence "{x. f x = 0} = {(-b + sqrt discr) / (2 * a), (-b - sqrt discr) / (2 * a)}"
      by blast
    thus ?th1 using True by simp
  next
    case False
    {
      fix x :: real
      assume x: "f x = 0"
      have "0 ≤ (x + b / (2 * a)) ^ 2" by simp
      also have "f x = a * ((x + b / (2 * a)) ^ 2 - b ^ 2 / (4 * a ^ 2) + c / a)"
        by (simp add: field_simps power2_eq_square f_def)
      with x have "(x + b / (2 * a)) ^ 2 - b ^ 2 / (4 * a ^ 2) + c / a = 0"
        by simp
      hence "(x + b / (2 * a)) ^ 2 = b ^ 2 / (4 * a ^ 2) - c / a"
        by (simp add: algebra_simps)
      finally have "0 ≤ (b2 / (4 * a2) - c / a) * (4 * a2)"
        by (intro mult_nonneg_nonneg) auto
      also have "… = b2 - 4 * a * c" by (simp add: field_simps power2_eq_square)
      also have "… < 0" using False by (simp add: discr_def)
      finally have False by simp
    }
    hence "{x. f x = 0} = {}" by auto
    thus ?th1 using False by simp
  qed
qed

lemma finite_quadratic_equation_solutions_reals:
  fixes a b c :: real
  defines "discr ≡ (b^2 - 4 * a * c)"
  shows   "finite {x. a * x ^ 2 + b * x + c = 0} ⟷ a ≠ 0 ∨ b ≠ 0 ∨ c ≠ 0"
  by (subst quadratic_equation_reals)
     (auto simp: discr_def card_eq_0_iff infinite_UNIV_char_0 split: if_split)

lemma card_quadratic_equation_solutions_reals:
  fixes a b c :: real
  defines "discr ≡ (b^2 - 4 * a * c)"
  shows   "card {x. a * x ^ 2 + b * x + c = 0} =
             (if a = 0 then
                (if b = 0 then 0 else 1)
              else if discr ≥ 0 then if discr = 0 then 1 else 2 else 0)" (is ?th1)                
  by (subst quadratic_equation_reals)
     (auto simp: discr_def card_eq_0_iff infinite_UNIV_char_0 split: if_split)

lemma card_quadratic_equation_solutions_reals_le_2:
  "card {x :: real. a * x ^ 2 + b * x + c = 0} ≤ 2"
  by (subst card_quadratic_equation_solutions_reals) auto

lemma quadratic_equation_solution_rat_iff:
  fixes a b c :: int and x y :: real
  defines "f ≡ (λx::real. a * x ^ 2 + b * x + c)"
  defines "discr ≡ nat (b ^ 2 - 4 * a * c)"
  assumes "a ≠ 0" "f x = 0"
  shows   "x ∈ ℚ ⟷ is_square discr"
proof -
  define discr' where "discr' ≡ real_of_int (b ^ 2 - 4 * a * c)"
  from assms have "x ∈ {x. f x = 0}" by simp
  with ‹a ≠ 0› have "discr' ≥ 0" unfolding discr'_def f_def of_nat_diff
    by (subst (asm) quadratic_equation_reals) (auto simp: discr_def split: if_splits)
  hence *: "sqrt (discr') = sqrt (real discr)" unfolding of_int_0_le_iff discr_def discr'_def
    by (simp add: algebra_simps nat_diff_distrib)
  from ‹x ∈ {x. f x = 0}› have "x = (-b + sqrt discr) / (2 * a) ∨ x = (-b - sqrt discr) / (2 * a)"
    using ‹a ≠ 0› * unfolding discr'_def f_def
    by (subst (asm) quadratic_equation_reals) (auto split: if_splits)
  thus ?thesis using ‹a ≠ 0›
    by (auto simp: sqrt_of_nat_in_Rats_iff divide_in_Rats_iff2 diff_in_Rats_iff2 diff_in_Rats_iff1)
qed


subsection ‹Periodic continued fractions and quadratic irrationals›

text ‹
  We now show the main result: A positive irrational number has a periodic continued 
  fraction expansion iff it is a quadratic irrational.

  In principle, this statement naturally also holds for negative numbers, but the current 
  formalisation of continued fractions only supports non-negative numbers. It also holds for
  rational numbers in some sense, since their continued fraction expansion is finite to begin with.
›
theorem periodic_cfrac_imp_quadratic_irrational:
  assumes [simp]: "cfrac_length c = ∞"
      and period: "l > 0" "⋀k. k ≥ N ⟹ cfrac_nth c (k + l) = cfrac_nth c k"
  shows   "quadratic_irrational (cfrac_lim c)"
proof -
  define h' and k' where "h' = conv_num_int (cfrac_drop N c)" 
                     and "k' = conv_denom_int (cfrac_drop N c)"
  define x' where "x' = cfrac_remainder c N"

  have c_pos: "cfrac_nth c n > 0" if "n ≥ N" for n
  proof -
    from assms(1,2) have "cfrac_nth c (n + l) > 0" by auto
    with assms(3)[OF that] show ?thesis by simp
  qed
  have k'_pos: "k' n > 0" if "n ≠ -1" "n ≥ -2" for n
    using that by (auto simp: k'_def conv_denom_int_def intro!: conv_denom_pos)
  have k'_nonneg: "k' n ≥ 0" if "n ≥ -2" for n
    using that by (auto simp: k'_def conv_denom_int_def intro!: conv_denom_pos)
  have "cfrac_nth c (n + (N + l)) = cfrac_nth c (n + N)" for n
    using period(2)[of "n + N"] by (simp add: add_ac)
  have "cfrac_drop (N + l) c = cfrac_drop N c"
    by (rule cfrac_eqI) (use period(2)[of "n + N" for n] in ‹auto simp: algebra_simps›)
  hence x'_altdef: "x' = cfrac_remainder c (N + l)"
    by (simp add: x'_def cfrac_remainder_def)
  have x'_pos: "x' > 0" unfolding x'_def
    using c_pos by (intro cfrac_remainder_pos) auto

  define A where "A = (k' (int l - 1))"
  define B where "B = k' (int l - 2) - h' (int l - 1)"
  define C where "C = -(h' (int l - 2))"

  have pos: "(k' (int l - 1) * x' + k' (int l - 2)) > 0"
    using x'_pos ‹l > 0›
    by (intro add_pos_nonneg mult_pos_pos) (auto intro!: k'_pos k'_nonneg)
  have "cfrac_remainder c N = conv' (cfrac_drop N c) l (cfrac_remainder c (l + N))"
    unfolding cfrac_remainder_def cfrac_drop_add
    by (subst (2) cfrac_remainder_def [symmetric]) (auto simp: conv'_cfrac_remainder)
  hence "x' = conv' (cfrac_drop N c) l x'"
    by (subst (asm) add.commute) (simp only: x'_def [symmetric] x'_altdef [symmetric])
  also have "… = (h' (int l - 1) * x' + h' (int l - 2)) / (k' (int l - 1) * x' + k' (int l - 2))"
    using conv'_num_denom_int[OF x'_pos, of _ l] unfolding h'_def k'_def
    by (simp add: mult_ac)
  finally have "x' * (k' (int l - 1) * x' + k' (int l - 2)) = (h' (int l - 1) * x' + h' (int l - 2))"
    using pos by (simp add: divide_simps)
  hence quadratic: "A * x' ^ 2 + B * x' + C = 0"
    by (simp add: algebra_simps power2_eq_square A_def B_def C_def)
  moreover have "x' ∉ ℚ" unfolding x'_def
    by auto
  moreover have "A > 0" using ‹l > 0› by (auto simp: A_def intro!: k'_pos)
  ultimately have "quadratic_irrational x'" using ‹x' ∉ ℚ›
    by (intro quadratic_irrational.intros[of x' A B C]) simp_all
  thus ?thesis
    using assms by (simp add: x'_def quadratic_irrational_cfrac_remainder_iff)
qed


lift_definition pperiodic_cfrac :: "nat list ⇒ cfrac" is
  "λxs. if xs = [] then (0, LNil) else
        (int (hd xs), llist_of_stream (cycle (map (λn. n- 1) (tl xs @ [hd xs]))))" .

definition periodic_cfrac :: "int list ⇒ int list ⇒ cfrac" where
  "periodic_cfrac xs ys = cfrac_of_stream (Stream.shift xs (Stream.cycle ys))"

lemma periodic_cfrac_Nil [simp]: "pperiodic_cfrac [] = 0"
  unfolding zero_cfrac_def by transfer auto

lemma cfrac_length_pperiodic_cfrac [simp]:
  "xs ≠ [] ⟹ cfrac_length (pperiodic_cfrac xs) = ∞"
  by transfer auto

lemma cfrac_nth_pperiodic_cfrac:
  assumes "xs ≠ []" and "0 ∉ set xs"
  shows   "cfrac_nth (pperiodic_cfrac xs) n = xs ! (n mod length xs)"
  using assms
proof (transfer, goal_cases)
  case (1 xs n)
  show ?case
  proof (cases n)
    case (Suc n')
    have "int (cycle (tl (map (λn. n - 1) xs) @ [hd (map (λn. n - 1) xs)]) !! n') + 1 = 
          int (stl (cycle (map (λn. n - 1) xs)) !! n') + 1"
      by (subst cycle.sel(2) [symmetric]) (rule refl)
    also have "… = int (cycle (map (λn. n - 1) xs) !! n) + 1"
      by (simp add: Suc del: cycle.sel)
    also have "… = int (xs ! (n mod length xs) - 1) + 1"
      by (simp add: snth_cycle ‹xs ≠ []›)
    also have "xs ! (n mod length xs) ∈ set xs"
      using ‹xs ≠ []› by (auto simp: set_conv_nth)
    with 1 have "xs ! (n mod length xs) > 0"
      by (intro Nat.gr0I) auto
    hence "int (xs ! (n mod length xs) - 1) + 1 = int (xs ! (n mod length xs))"
      by simp
    finally show ?thesis
      using Suc 1 by (simp add: hd_conv_nth map_tl)
  qed (use 1 in ‹auto simp: hd_conv_nth›)
qed

definition pperiodic_cfrac_info :: "nat list ⇒ int × int × int"where
  "pperiodic_cfrac_info xs =
     (let l = length xs;
          h = conv_num_fun (λn. xs ! n);
          k = conv_denom_fun (λn. xs ! n);
          A = k (l - 1);
          B = h (l - 1) - (if l = 1 then 0 else k (l - 2));
          C = (if l = 1 then -1 else -h (l - 2))
      in  (B^2-4*A*C, B, 2 * A))"

lemma conv_gen_cong:
  assumes "∀k∈{n..N}. f k = f' k"
  shows   "conv_gen f (a,b,n) N = conv_gen f' (a,b,n) N"
  using assms 
proof (induction "N - n" arbitrary: a b n N)
  case (Suc d n N a b)
  have "conv_gen f (b, b * f n + a, Suc n) N = conv_gen f' (b, b * f n + a, Suc n) N"
    using Suc(2,3) by (intro Suc) auto
  moreover have "f n = f' n"
    using bspec[OF Suc.prems, of n] Suc(2) by auto
  ultimately show ?case
    by (subst (1 2) conv_gen.simps) auto
qed (auto simp: conv_gen.simps)

lemma
  assumes "∀k≤n. c k = cfrac_nth c' k"
  shows   conv_num_fun_eq': "conv_num_fun c n = conv_num c' n"
    and   conv_denom_fun_eq': "conv_denom_fun c n = conv_denom c' n"
proof -
  have "conv_num c' n = conv_gen (cfrac_nth c') (0, 1, 0) n"
    unfolding conv_num_code ..
  also have "… = conv_gen c (0, 1, 0) n"
    unfolding conv_num_fun_def using assms by (intro conv_gen_cong) auto
  finally show "conv_num_fun c n = conv_num c' n"
    by (simp add: conv_num_fun_def)
next
  have "conv_denom c' n = conv_gen (cfrac_nth c') (1, 0, 0) n"
    unfolding conv_denom_code ..
  also have "… = conv_gen c (1, 0, 0) n"
    unfolding conv_denom_fun_def using assms by (intro conv_gen_cong) auto
  finally show "conv_denom_fun c n = conv_denom c' n"
    by (simp add: conv_denom_fun_def)
qed

lemma gcd_minus_commute_left: "gcd (a - b :: 'a :: ring_gcd) c = gcd (b - a) c"
  by (metis gcd.commute gcd_neg2 minus_diff_eq)

lemma gcd_minus_commute_right: "gcd c (a - b :: 'a :: ring_gcd) = gcd c (b - a)"
  by (metis gcd_neg2 minus_diff_eq)

lemma periodic_cfrac_info_aux:
  fixes D E F :: int
  assumes "pperiodic_cfrac_info xs = (D, E, F)"
  assumes "xs ≠ []" "0 ∉ set xs"
  shows   "cfrac_lim (pperiodic_cfrac xs) = (sqrt D + E) / F"
    and   "D > 0" and "F > 0"
proof -
  define c where "c = pperiodic_cfrac xs"
  have [simp]: "cfrac_length c = ∞"
    using assms by (simp add: c_def)
  define h and k where "h = conv_num_int c" and "k = conv_denom_int c"
  define x where "x = cfrac_lim c"
  define l where "l = length xs"

  define A where "A = (k (int l - 1))"
  define B where "B = k (int l - 2) - h (int l - 1)"
  define C where "C = -(h (int l - 2))"
  define discr where "discr = B ^ 2 - 4 * A * C"

  have "l > 0"
    using assms by (simp add: l_def)
  have c_pos: "cfrac_nth c n > 0" for n
    using assms by (auto simp: c_def cfrac_nth_pperiodic_cfrac set_conv_nth)
  have x_pos: "x > 0"
    unfolding x_def by (intro cfrac_lim_pos c_pos)
  have h_pos: "h n > 0" if "n > -2" for n
    using that unfolding h_def by (auto simp: conv_num_int_def intro: conv_num_pos' c_pos)
  have k_pos: "k n > 0" if "n > -1" for n
    using that unfolding k_def by (auto simp: conv_denom_int_def)
  have k_nonneg: "k n ≥ 0" for n
    unfolding k_def by (auto simp: conv_denom_int_def)

  have pos: "(k (int l - 1) * x + k (int l - 2)) > 0"
    using x_pos ‹l > 0›
    by (intro add_pos_nonneg mult_pos_pos) (auto intro!: k_pos k_nonneg)
  have "cfrac_drop l c = c"
    using assms by (intro cfrac_eqI) (auto simp: c_def cfrac_nth_pperiodic_cfrac l_def)

  have "x = conv' c l (cfrac_remainder c l)"
    unfolding x_def by (rule conv'_cfrac_remainder[symmetric]) auto
  also have "… = conv' c l x"
    unfolding cfrac_remainder_def ‹cfrac_drop l c = c› x_def ..
  finally have "x = conv' c l x" .
  also have "… = (h (int l - 1) * x + h (int l - 2)) / (k (int l - 1) * x + k (int l - 2))"
    using conv'_num_denom_int[OF x_pos, of _ l] unfolding h_def k_def
    by (simp add: mult_ac)
  finally have "x * (k (int l - 1) * x + k (int l - 2)) = (h (int l - 1) * x + h (int l - 2))"
    using pos by (simp add: divide_simps)
  hence quadratic: "A * x ^ 2 + B * x + C = 0"
    by (simp add: algebra_simps power2_eq_square A_def B_def C_def)

  have "A > 0" using ‹l > 0› by (auto simp: A_def intro!: k_pos)
  have discr_altdef: "discr = (k (int l-2) - h (int l-1)) ^ 2 + 4 * k (int l-1) * h (int l-2)"
    by (simp add: discr_def A_def B_def C_def)

  have "0 < 0 + 4 * A * 1"
    using ‹A > 0› by simp
  also have "0 + 4 * A * 1 ≤ discr"
    unfolding discr_altdef A_def using h_pos[of "int l - 2"] ‹l > 0›
    by (intro add_mono mult_mono order.refl k_nonneg mult_nonneg_nonneg) auto
  finally have "discr > 0" .

  have "x ∈ {x. A * x ^ 2 + B * x + C = 0}"
    using quadratic by simp
  hence x_cases: "x = (-B - sqrt discr) / (2 * A) ∨ x = (-B + sqrt discr) / (2 * A)"
    unfolding quadratic_equation_reals of_int_diff using ‹A > 0›
    by (auto split: if_splits simp: discr_def)

  have "B ^ 2 < discr"
    unfolding discr_def by (auto intro!: mult_pos_pos k_pos h_pos ‹l > 0› simp: A_def C_def)
  hence "¦B¦ < sqrt discr"
    using ‹discr > 0› by (simp add: real_less_rsqrt)

  have "x = (if x ≥ 0 then (sqrt discr - B) / (2 * A) else -(sqrt discr + B) / (2 * A))"
    using x_cases
  proof
    assume x: "x = (-B - sqrt discr) / (2 * A)"
    have "(-B - sqrt discr) / (2 * A) < 0"
      using ‹¦B¦ < sqrt discr› ‹A > 0› by (intro divide_neg_pos) auto
    also note x[symmetric]
    finally show ?thesis using x by simp
  next
    assume x: "x = (-B + sqrt discr) / (2 * A)"
    have "(-B + sqrt discr) / (2 * A) > 0"
      using ‹¦B¦ < sqrt discr› ‹A > 0› by (intro divide_pos_pos) auto
    also note x[symmetric]
    finally show ?thesis using x by simp
  qed
  also have "x ≥ 0 ⟷ floor x ≥ 0"
    by auto
  also have "floor x = floor (cfrac_lim c)"
    by (simp add: x_def)
  also have "… = cfrac_nth c 0"
    by (subst cfrac_nth_0_conv_floor) auto
  also have "… = int (hd xs)"
    using assms unfolding c_def by (subst cfrac_nth_pperiodic_cfrac) (auto simp: hd_conv_nth)
  finally have x_eq: "x = (sqrt discr - B) / (2 * A)"
    by simp


  define h' where "h' = conv_num_fun (λn. int (xs ! n))"
  define k' where "k' = conv_denom_fun (λn. int (xs ! n))"
  have num_eq: "h' i = h i"
    if "i < l" for i using that assms unfolding h'_def h_def
    by (subst conv_num_fun_eq'[where c' = c]) (auto simp: c_def l_def cfrac_nth_pperiodic_cfrac)
  have denom_eq: "k' i = k i"
    if "i < l" for i using that assms unfolding k'_def k_def
    by (subst conv_denom_fun_eq'[where c' = c]) (auto simp: c_def l_def cfrac_nth_pperiodic_cfrac)

  have 1: "h (int l - 1) = h' (l - 1)"
    by (subst num_eq) (use ‹l > 0› in ‹auto simp: of_nat_diff›)
  have 2: "k (int l - 1) = k' (l - 1)"
    by (subst denom_eq) (use ‹l > 0› in ‹auto simp: of_nat_diff›)
  have 3: "h (int l - 2) = (if l = 1 then 1 else h' (l - 2))"
    using ‹l > 0› num_eq[of "l - 2"] by (auto simp: h_def nat_diff_distrib)
  have 4: "k (int l - 2) = (if l = 1 then 0 else k' (l - 2))"
    using ‹l > 0› denom_eq[of "l - 2"] by (auto simp: k_def nat_diff_distrib)
  
  have "pperiodic_cfrac_info xs =
          (let A = k (int l - 1);
               B = h (int l - 1) - (if l = 1 then 0 else k (int l - 2));
               C = (if l = 1 then -1 else - h (int l - 2))
           in (B2 - 4 * A * C, B, 2 * A))"
    unfolding pperiodic_cfrac_info_def Let_def using 1 2 3 4 ‹l > 0›
    by (auto simp: num_eq denom_eq h'_def k'_def l_def of_nat_diff)
  also have "… = (B2 - 4 * A * C, -B, 2 * A)"
    by (simp add: Let_def A_def B_def C_def h_def k_def algebra_simps power2_commute)
  finally have per_eq: "pperiodic_cfrac_info xs = (discr, -B, 2 * A)"
    by (simp add: discr_def)

  show "x = (sqrt (real_of_int D) + real_of_int E) / real_of_int F"
    using per_eq assms by (simp add: x_eq)
  show "D > 0" "F > 0"
    using assms per_eq ‹discr > 0› ‹A > 0› by auto
qed

text ‹
  We can now compute surd representations for (purely) periodic continued fractions, e.g.
  $[1, 1, 1, \ldots] = \frac{\sqrt{5} + 1}{2}$:
›
value "pperiodic_cfrac_info [1]"

text ‹
  We can now compute surd representations for periodic continued fractions, e.g.
  $[\overline{1,1,1,1,6}] = \frac{\sqrt{13} + 3}{4}$:
›
value "pperiodic_cfrac_info [1,1,1,1,6]"



text ‹
  With a little bit of work, one could also easily derive from this a version for
  non-purely periodic continued fraction.
›


text ‹
  Next, we show that any quadratic irrational has a periodic continued fraction expansion.
›
theorem quadratic_irrational_imp_periodic_cfrac:
  assumes "quadratic_irrational (cfrac_lim e)"
  obtains N l where "l > 0" and "⋀n m. n ≥ N ⟹ cfrac_nth e (n + m * l) = cfrac_nth e n"
                and "cfrac_remainder e (N + l) = cfrac_remainder e N"
                and "cfrac_length e = ∞"
proof -
  have [simp]: "cfrac_length e = ∞"
    using assms by (auto simp: quadratic_irrational.simps)
  note [intro] = assms(1)
  define x where "x = cfrac_lim e"
  from assms obtain a b c :: int where
    nontrivial: "a ≠ 0 ∨ b ≠ 0 ∨ c ≠ 0" and
          root: "a * x^2 + b * x + c = 0" (is "?f x = 0")
    by (auto simp: quadratic_irrational.simps x_def)

  define f where "f = ?f"
  define h and k where "h = conv_num e" and "k = conv_denom e"
  define X where "X = cfrac_remainder e"
  have [simp]: "k i > 0" "k i ≠ 0" for i
    using conv_denom_pos[of e i] by (auto simp: k_def)
  have k_leI: "k i ≤ k j" if "i ≤ j" for i j
    by (auto simp: k_def intro!: conv_denom_leI that)
  have k_nonneg: "k n ≥ 0" for n
    by (auto simp: k_def)
  have k_ge_1: "k n ≥ 1" for n
    using k_leI[of 0 n] by (simp add: k_def)
    
  define R where "R = conv e"
  define A where "A = (λn. a * h (n - 1) ^ 2 + b * h (n - 1) * k (n - 1) + c * k (n - 1) ^ 2)"
  define B where "B = (λn. 2 * a * h (n - 1) * h (n - 2) + b * (h (n - 1) * k (n - 2) + h (n - 2) * k (n - 1)) + 2 * c * k (n - 1) * k (n - 2))"
  define C where "C = (λn. a * h (n - 2) ^ 2 + b * h (n - 2) * k (n - 2) + c * k (n - 2) ^ 2)"

  define A' where "A' = nat ⌊2 * ¦a¦ * ¦x¦ + ¦a¦ + ¦b¦⌋"
  define B' where "B' = nat ⌊(3 / 2) * (2 * ¦a¦ * ¦x¦ + ¦b¦) + 9 / 4 * ¦a¦⌋"

  have [simp]: "X n ∉ ℚ" for n unfolding X_def
    by simp
  from this[of 0] have [simp]: "x ∉ ℚ"
    unfolding X_def by (simp add: x_def)

  have "a ≠ 0"
  proof
    assume "a = 0"
    with root and nontrivial have "x = 0 ∨ x = -c / b"
      by (auto simp: divide_simps add_eq_0_iff)
    hence "x ∈ ℚ" by (auto simp del: ‹x ∉ ℚ›)
    thus False by simp
  qed

  have bounds: "(A n, B n, C n) ∈ {-A'..A'} × {-B'..B'} × {-A'..A'}"
   and X_root: "A n * X n ^ 2 + B n * X n + C n = 0" if n: "n ≥ 2" for n
  proof -
    define n' where "n' = n - 2"
    have n': "n = Suc (Suc n')" using ‹n ≥ 2› unfolding n'_def by simp
    have *: "of_int (k (n - Suc 0)) * X n + of_int (k (n - 2)) ≠ 0"
    proof
      assume "of_int (k (n - Suc 0)) * X n + of_int (k (n - 2)) = 0"
      hence "X n = -k (n - 2) / k (n - 1)" by (auto simp: divide_simps mult_ac)
      also have "… ∈ ℚ" by auto
      finally show False by simp
    qed
  
    let ?denom = "(k (n - 1) * X n + k (n - 2))"
    have "0 = 0 * ?denom ^ 2" by simp
    also have "0 * ?denom ^ 2 = (a * x ^ 2 + b * x + c) * ?denom ^ 2" using root by simp
    also have "… = a * (x * ?denom) ^ 2 + b * ?denom * (x * ?denom) + c * ?denom * ?denom"
      by (simp add: algebra_simps power2_eq_square)
    also have "x * ?denom = h (n - 1) * X n + h (n - 2)"
      using cfrac_lim_eq_num_denom_remainder_aux[of "n - 2" e] ‹n ≥ 2›
      by (simp add: numeral_2_eq_2 Suc_diff_Suc x_def k_def h_def X_def)
    also have "a * … ^ 2 + b * ?denom * … + c * ?denom * ?denom = A n * X n ^ 2 + B n * X n + C n"
      by (simp add: A_def B_def C_def power2_eq_square algebra_simps)
    finally show "A n * X n ^ 2 + B n * X n + C n = 0" ..
  
    have f_abs_bound: "¦f (R n)¦ ≤ (2 * ¦a¦ * ¦x¦ + ¦b¦) * (1 / (k n * k (Suc n))) +
                                      ¦a¦ * (1 / (k n * k (Suc n))) ^ 2" for n
    proof -
      have "¦f (R n)¦ = ¦?f (R n) - ?f x¦" by (simp add: root f_def)
      also have "?f (R n) - ?f x = (R n - x) * (2 * a * x + b) + (R n - x) ^ 2 * a"
        by (simp add: power2_eq_square algebra_simps)
      also have "¦…¦ ≤ ¦(R n - x) * (2 * a * x + b)¦ + ¦(R n - x) ^ 2 * a¦"
        by (rule abs_triangle_ineq)
      also have "… = ¦2 * a * x + b¦ * ¦R n - x¦ + ¦a¦ * ¦R n - x¦ ^ 2"
        by (simp add: abs_mult)
      also have "… ≤ ¦2 * a * x + b¦ * (1 / (k n * k (Suc n))) + ¦a¦ * (1 / (k n * k (Suc n))) ^ 2"
        unfolding x_def R_def using cfrac_lim_minus_conv_bounds[of n e]
        by (intro add_mono mult_left_mono power_mono) (auto simp: k_def)
      also have "¦2 * a * x + b¦ ≤ 2 * ¦a¦ * ¦x¦ + ¦b¦"
        by (rule order.trans[OF abs_triangle_ineq]) (auto simp: abs_mult)
      hence "¦2 * a * x + b¦ * (1 / (k n * k (Suc n))) + ¦a¦ * (1 / (k n * k (Suc n))) ^ 2 ≤
               … * (1 / (k n * k (Suc n))) + ¦a¦ * (1 / (k n * k (Suc n))) ^ 2"
        by (intro add_mono mult_right_mono) (auto intro!: mult_nonneg_nonneg k_nonneg)
      finally show "¦f (R n)¦ ≤ …"
        by (simp add: mult_right_mono add_mono divide_left_mono)
    qed
  
    have h_eq_conv_k: "h i = R i * k i" for i
      using conv_denom_pos[of e i] unfolding R_def
      by (subst conv_num_denom) (auto simp: h_def k_def)
  
    have "A n = k (n - 1) ^ 2 * f (R (n - 1))" for n
      by (simp add: algebra_simps A_def n' k_def power2_eq_square h_eq_conv_k f_def)
    have A_bound: "¦A i¦ ≤ A'" if "i > 0" for i
    proof -
      have "k i > 0"
        by simp
      hence "k i ≥ 1"
        by linarith
      have "A i = k (i - 1) ^ 2 * f (R (i - 1))"
        by (simp add: algebra_simps A_def k_def power2_eq_square h_eq_conv_k f_def)
      also have "¦…¦ = k (i - 1) ^ 2 * ¦f (R (i - 1))¦"
        by (simp add: abs_mult f_def)
      also have "… ≤ k (i - 1) ^ 2 * ((2 * ¦a¦ * ¦x¦ + ¦b¦) * (1 / (k (i - 1) * k (Suc (i - 1)))) +
                        ¦a¦ * (1 / (k (i - 1) * k (Suc (i - 1)))) ^ 2)"
        by (intro mult_left_mono f_abs_bound) auto
      also have "… = k (i - 1) / k i * (2 * ¦a¦ * ¦x¦ + ¦b¦) + ¦a¦ / k i ^ 2" using ‹i > 0›
        by (simp add: power2_eq_square field_simps)
      also have "… ≤ 1 * (2 * ¦a¦ * ¦x¦ + ¦b¦) + ¦a¦ / 1" using ‹i > 0› ‹k i ≥ 1›
        by (intro add_mono divide_left_mono mult_right_mono)
           (auto intro!: k_leI one_le_power simp: of_nat_ge_1_iff)
      also have "… = 2 * ¦a¦ * ¦x¦ + ¦a¦ + ¦b¦" by simp
      finally show ?thesis unfolding A'_def by linarith
    qed
  
    have "C n = A (n - 1)" by (simp add: A_def C_def n')
    hence C_bound: "¦C n¦ ≤ A'" using A_bound[of "n - 1"] n by simp
  
    have "B n = k (n - 1) * k (n - 2) *
                  (f (R (n - 1)) + f (R (n - 2)) - a * (R (n - 1) - R (n - 2)) ^ 2)"
      by (simp add: B_def h_eq_conv_k algebra_simps power2_eq_square f_def)
    also have "¦…¦ = k (n - 1) * k (n - 2) * 
                       ¦f (R (n - 1)) + f (R (n - 2)) - a * (R (n - 1) - R (n - 2)) ^ 2¦"
      by (simp add: abs_mult k_nonneg)
    also have "… ≤ k (n - 1) * k (n - 2) * 
                      (((2 * ¦a¦ * ¦x¦ + ¦b¦) * (1 / (k (n - 1) * k (Suc (n - 1)))) +
                          ¦a¦ * (1 / (k (n - 1) * k (Suc (n - 1)))) ^ 2) +                      
                       ((2 * ¦a¦ * ¦x¦ + ¦b¦) * (1 / (k (n - 2) * k (Suc (n - 2)))) +
                          ¦a¦ * (1 / (k (n - 2) * k (Suc (n - 2)))) ^ 2) +
                        ¦a¦ * ¦R (Suc (n - 2)) - R (n - 2)¦ ^ 2)" (is "_ ≤ _ * (?S1 + ?S2 + ?S3)")
      by (intro mult_left_mono order.trans[OF abs_triangle_ineq4] order.trans[OF abs_triangle_ineq] 
            add_mono f_abs_bound order.refl)
         (insert n, auto simp: abs_mult Suc_diff_Suc numeral_2_eq_2 k_nonneg)
    also have "¦R (Suc (n - 2)) - R (n - 2)¦ = 1 / (k (n - 2) * k (Suc (n - 2)))"
      unfolding R_def k_def by (rule abs_diff_successive_convs)
    also have "of_int (k (n - 1) * k (n - 2)) * (?S1 + ?S2 + ¦a¦ * … ^ 2) = 
                 (k (n - 2) / k n + 1) * (2 * ¦a¦ * ¦x¦ + ¦b¦) + 
                 ¦a¦ * (k (n - 2) / (k (n - 1) * k n ^ 2) + 2 / (k (n - 1) * k (n - 2)))"
      (is "_ = ?S") using n by (simp add: field_simps power2_eq_square numeral_2_eq_2 Suc_diff_Suc)
    also {
      have A: "2 * real_of_int (k (n - 2)) ≤ of_int (k n)"
        using conv_denom_plus2_ratio_ge[of e "n - 2"] n
        by (simp add: numeral_2_eq_2 Suc_diff_Suc k_def)
      have "fib (Suc 2) ≤ k 2" unfolding k_def by (intro conv_denom_lower_bound)
      also have "… ≤ k n" by (intro k_leI n)
      finally have "k n ≥ 2" by (simp add: numeral_3_eq_3)
      hence B: "of_int (k (n - 2)) * 2 ^ 2 ≤ (of_int (k (n - 1)) * (of_int (k n))2 :: real)"
        by (intro mult_mono power_mono) (auto intro: k_leI k_nonneg)
      have C: "1 * 1 ≤ real_of_int (k (n - 1)) * of_int (k (n - 2))" using k_ge_1
        by (intro mult_mono) (auto simp: Suc_le_eq of_nat_ge_1_iff k_nonneg)
      note A B C
    }
    hence "?S ≤ (1 / 2 + 1) * (2 * ¦a¦ * ¦x¦ + ¦b¦) + ¦a¦ * (1 / 4 + 2)"
      by (intro add_mono mult_right_mono mult_left_mono) (auto simp: field_simps)
    also have "… = (3 / 2) * (2 * ¦a¦ * ¦x¦ + ¦b¦) + 9 / 4 * ¦a¦" by simp
    finally have B_bound: "¦B n¦ ≤ B'" unfolding B'_def by linarith
    from A_bound[of n] B_bound C_bound n
    show "(A n, B n, C n) ∈ {-A'..A'} × {-B'..B'} × {-A'..A'}" by auto
  qed

  have A_nz: "A n ≠ 0" if "n ≥ 1" for n
    using that
  proof (induction n rule: dec_induct)
    case base
    show ?case
    proof
      assume "A 1 = 0"
      hence "real_of_int (A 1) = 0" by simp
      also have "real_of_int (A 1) =
                   real_of_int a * of_int (cfrac_nth e 0) ^ 2 +
                   real_of_int b * cfrac_nth e 0 + real_of_int c"
        by (simp add: A_def h_def k_def)
      finally have root': "… = 0" .

      have "cfrac_nth e 0 ∈ ℚ" by auto
      also from root' and ‹a ≠ 0› have "?this ⟷ is_square (nat (b2 - 4 * a * c))"
        by (intro quadratic_equation_solution_rat_iff) auto
      also from root and ‹a ≠ 0› have "… ⟷ x ∈ ℚ"
        by (intro quadratic_equation_solution_rat_iff [symmetric]) auto
      finally show False using ‹x ∉ ℚ› by contradiction
    qed
  next
    case (step m)
    hence nz: "C (Suc m) ≠ 0" by (simp add: C_def A_def)
    show "A (Suc m) ≠ 0"
    proof
      assume [simp]: "A (Suc m) = 0"
      have "X (Suc m) > 0" unfolding X_def
        by (intro cfrac_remainder_pos) auto
      with X_root[of "Suc m"] step.hyps nz have "X (Suc m) = -C (Suc m) / B (Suc m)"
        by (auto simp: divide_simps mult_ac)
      also have "… ∈ ℚ" by auto
      finally show False by simp
    qed
  qed 

  have "finite ({-A'..A'} × {-B'..B'} × {-A'..A'})" by auto
  from this and bounds have "finite ((λn. (A n, B n, C n)) ` {2..})"
    by (blast intro: finite_subset)
  moreover have "infinite ({2..} :: nat set)" by (simp add: infinite_Ici)
  ultimately have "∃k1∈{2..}. infinite {n ∈ {2..}. (A n, B n, C n) = (A k1, B k1, C k1)}"
    by (intro pigeonhole_infinite)
  then obtain k0 where k0: "k0 ≥ 2" "infinite {n ∈ {2..}. (A n, B n, C n) = (A k0, B k0, C k0)}"
    by auto
  from infinite_countable_subset[OF this(2)] obtain g :: "nat ⇒ _"
    where g: "inj g" "range g ⊆ {n∈{2..}. (A n, B n, C n) = (A k0, B k0, C k0)}" by blast
  hence g_ge_2: "g k ≥ 2" for k by auto
  from g have [simp]: "A (g k) = A k0" "B (g k) = B k0" "C (g k) = C k0" for k
    by auto

  from g(1) have [simp]: "g k1 = g k2 ⟷ k1 = k2" for k1 k2 by (auto simp: inj_def)
  define z where "z = (A k0, B k0, C k0)"
  let ?h = "λk. (A (g k), B (g k), C (g k))"
  from g have g': "distinct [g 1, g 2, g 3]" "?h 0 = z" "?h 1 = z" "?h 2 = z"
    by (auto simp: z_def)
  have fin: "finite {x :: real. A k0 * x ^ 2 + B k0 * x + C k0 = 0}" using A_nz[of k0] k0(1)
    by (subst finite_quadratic_equation_solutions_reals) auto
  from X_root[of "g 0"] X_root[of "g 1"] X_root[of "g 2"] g_ge_2 g
    have "(X ∘ g) ` {0, 1, 2} ⊆ {x. A k0 * x ^ 2 + B k0 * x + C k0 = 0}"
    by auto
  hence "card ((X ∘ g) ` {0, 1, 2}) ≤ card …"
    by (intro card_mono fin) auto
  also have "… ≤ 2"
    by (rule card_quadratic_equation_solutions_reals_le_2)
  also have "… < card {0, 1, 2 :: nat}" by simp
  finally have "¬inj_on (X ∘ g) {0, 1, 2}"
    by (rule pigeonhole)
  then obtain m1 m2 where
    m12: "m1 ∈ {0, 1, 2}" "m2 ∈ {0, 1, 2}" "X (g m1) = X (g m2)" "m1 ≠ m2"
    unfolding inj_on_def o_def by blast
  define n and l where "n = min (g m1) (g m2)" and "l = nat ¦int (g m1) - g m2¦"
  with m12 g' have l: "l > 0" "X (n + l) = X n"
    by (auto simp: min_def nat_diff_distrib split: if_splits)

  from l have "cfrac_lim (cfrac_drop (n + l) e) = cfrac_lim (cfrac_drop n e)"
    by (simp add: X_def cfrac_remainder_def)
  hence "cfrac_drop (n + l) e = cfrac_drop n e"
    by (simp add: cfrac_lim_eq_iff)
  hence "cfrac_nth (cfrac_drop (n + l) e) = cfrac_nth (cfrac_drop n e)"
    by (simp only:)
  hence period: "cfrac_nth e (n + l + k) = cfrac_nth e (n + k)" for k
    by (simp add: fun_eq_iff add_ac)
  have period: "cfrac_nth e (k + l) = cfrac_nth e k" if "k ≥ n" for k
    using period[of "k - n"] that by (simp add: add_ac)
  have period: "cfrac_nth e (k + m * l) = cfrac_nth e k" if "k ≥ n" for k m
    using that
  proof (induction m)
    case (Suc m)
    have "cfrac_nth e (k + Suc m * l) = cfrac_nth e (k + m * l + l)"
      by (simp add: algebra_simps)
    also have "… = cfrac_nth e (k + m * l)"
      using Suc.prems by (intro period) auto
    also have "… = cfrac_nth e k"
      using Suc.prems by (intro Suc.IH) auto
    finally show ?case .
  qed simp_all

  from this and l and that[of l n] show ?thesis by (simp add: X_def)
qed

theorem periodic_cfrac_iff_quadratic_irrational:
  assumes "x ∉ ℚ" "x ≥ 0"
  shows   "quadratic_irrational x ⟷ 
           (∃N l. l > 0 ∧ (∀n≥N. cfrac_nth (cfrac_of_real x) (n + l) =
                                  cfrac_nth (cfrac_of_real x) n))"
proof safe
  assume *: "quadratic_irrational x"
  with assms have **: "quadratic_irrational (cfrac_lim (cfrac_of_real x))" by auto
  obtain N l where Nl: "l > 0"
    "⋀n m. N ≤ n ⟹ cfrac_nth (cfrac_of_real x) (n + m * l) = cfrac_nth (cfrac_of_real x) n"
    "cfrac_remainder (cfrac_of_real x) (N + l) = cfrac_remainder (cfrac_of_real x) N"
    "cfrac_length (cfrac_of_real x) = ∞"
    using quadratic_irrational_imp_periodic_cfrac [OF **] by metis
  show "∃N l. l > 0 ∧ (∀n≥N. cfrac_nth (cfrac_of_real x) (n + l) = cfrac_nth (cfrac_of_real x) n)"
    by (rule exI[of _ N], rule exI[of _ l]) (insert Nl(1) Nl(2)[of _ 1], auto)
next
  fix N l assume "l > 0" "∀n≥N. cfrac_nth (cfrac_of_real x) (n + l) = cfrac_nth (cfrac_of_real x) n"
  hence "quadratic_irrational (cfrac_lim (cfrac_of_real x))" using assms
    by (intro periodic_cfrac_imp_quadratic_irrational[of _ l N]) auto
  with assms show "quadratic_irrational x"
    by simp
qed

text ‹
  The following result can e.g. be used to show that a number is ∗‹not› a quadratic
  irrational.
›
lemma quadratic_irrational_cfrac_nth_range_finite:
  assumes "quadratic_irrational (cfrac_lim e)"
  shows   "finite (range (cfrac_nth e))" 
proof -
  from quadratic_irrational_imp_periodic_cfrac[OF assms] obtain l N
    where period: "l > 0" "⋀m n. n ≥ N ⟹ cfrac_nth e (n + m * l) = cfrac_nth e n"
    by metis
  have "cfrac_nth e k ∈ cfrac_nth e ` {..<N+l}" for k
  proof (cases "k < N + l")
    case False
    define n m where "n = N + (k - N) mod l" and "m = (k - N) div l"
    have "cfrac_nth e n ∈ cfrac_nth e ` {..<N+l}"
      using ‹l > 0› by (intro imageI) (auto simp: n_def)
    also have "cfrac_nth e n = cfrac_nth e (n + m * l)"
      by (subst period) (auto simp: n_def)
    also have "n + m * l = k"
      using False by (simp add: n_def m_def)
    finally show ?thesis .
  qed auto
  hence "range (cfrac_nth e) ⊆ cfrac_nth e ` {..<N+l}"
    by blast
  thus ?thesis by (rule finite_subset) auto
qed

end