Theory 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"
  "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: "Discrete.sqrt n = nat sqrt n"
proof -
  have "real (Discrete.sqrt n ^ 2)  sqrt n ^ 2"
    by simp
  hence "Discrete.sqrt n  sqrt n"
    unfolding of_nat_power by (rule power2_le_imp_le) auto
  moreover have "real (Suc (Discrete.sqrt n) ^ 2) > real n"
    unfolding of_nat_less_iff by (rule Suc_sqrt_power2_gt)
  hence "real (Discrete.sqrt n + 1) ^ 2 > sqrt n ^ 2"
    unfolding of_nat_power by simp
  hence "real (Discrete.sqrt n + 1) > sqrt n"
    by (rule power2_less_imp_less) auto
  hence "Discrete.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 "kn. 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  (nN. 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  (nN. 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" "nN. 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