Theory Sqrt_Nat_Cfrac
section ‹Continued fraction expansions for square roots of naturals›
theory Sqrt_Nat_Cfrac
imports
Quadratic_Irrationals
"HOL-Library.While_Combinator"
"HOL-Library.IArray"
begin
text ‹
In this section, we shall explore the continued fraction expansion of $\sqrt{D}$, where $D$
is a natural number.
›
lemma butlast_nth [simp]: "n < length xs - 1 ⟹ butlast xs ! n = xs ! n"
by (induction xs arbitrary: n) (auto simp: nth_Cons split: nat.splits)
text ‹
The following is the length of the period in the continued fraction expansion of
$\sqrt{D}$ for a natural number $D$.
›
definition sqrt_nat_period_length :: "nat ⇒ nat" where
"sqrt_nat_period_length D =
(if is_square D then 0
else (LEAST l. l > 0 ∧ (∀n. cfrac_nth (cfrac_of_real (sqrt D)) (Suc n + l) =
cfrac_nth (cfrac_of_real (sqrt D)) (Suc n))))"
text ‹
Next, we define a more workable representation for the continued fraction expansion of
$\sqrt{D}$ consisting of the period length, the natural number $\lfloor\sqrt{D}\rfloor$, and
the content of the period.
›
definition sqrt_cfrac_info :: "nat ⇒ nat × nat × nat list" where
"sqrt_cfrac_info D =
(sqrt_nat_period_length D, Discrete.sqrt D,
map (λn. nat (cfrac_nth (cfrac_of_real (sqrt D)) (Suc n))) [0..<sqrt_nat_period_length D])"
lemma sqrt_nat_period_length_square [simp]: "is_square D ⟹ sqrt_nat_period_length D = 0"
by (auto simp: sqrt_nat_period_length_def)
definition sqrt_cfrac :: "nat ⇒ cfrac"
where "sqrt_cfrac D = cfrac_of_real (sqrt (real D))"
context
fixes D D' :: nat
defines "D' ≡ nat ⌊sqrt D⌋"
begin
text ‹
A number $\alpha = \frac{\sqrt D + p}{q}$ for ‹p, q ∈ ℕ› is called a ∗‹reduced quadratic surd›
if $\alpha > 1$ and $bar\alpha \in (-1;0)$, where $\bar\alpha$ denotes the conjugate
$\frac{-\sqrt D + p}{q}$.
It is furthermore called ∗‹associated› to $D$ if ‹q› divides ‹D - p⇧2›.
›
definition red_assoc :: "nat × nat ⇒ bool" where
"red_assoc = (λ(p, q).
q > 0 ∧ q dvd (D - p⇧2) ∧ (sqrt D + p) / q > 1 ∧ (-sqrt D + p) / q ∈ {-1<..<0})"
text ‹
The following two functions convert between a surd represented as a pair of natural numbers
and the actual real number and its conjugate:
›
definition surd_to_real :: "nat × nat ⇒ real"
where "surd_to_real = (λ(p, q). (sqrt D + p) / q)"
definition surd_to_real_cnj :: "nat × nat ⇒ real"
where "surd_to_real_cnj = (λ(p, q). (-sqrt D + p) / q)"
text ‹
The next function performs a single step in the continued fraction expansion of $\sqrt{D}$.
›
definition sqrt_remainder_step :: "nat × nat ⇒ nat × nat" where
"sqrt_remainder_step = (λ(p, q). let X = (p + D') div q; p' = X * q - p in (p', (D - p'⇧2) div q))"
text ‹
If we iterate this step function starting from the surd
$\frac{1}{\sqrt{D} - \lfloor\sqrt{D}\rfloor}$, we get the entire expansion.
›
definition sqrt_remainder_surd :: "nat ⇒ nat × nat"
where "sqrt_remainder_surd = (λn. (sqrt_remainder_step ^^ n) (D', D - D'⇧2))"
context
fixes sqrt_cfrac_nth :: "nat ⇒ nat" and l
assumes nonsquare: "¬is_square D"
defines "sqrt_cfrac_nth ≡ (λn. case sqrt_remainder_surd n of (p, q) ⇒ (D' + p) div q)"
defines "l ≡ sqrt_nat_period_length D"
begin
lemma D'_pos: "D' > 0"
using nonsquare by (auto simp: D'_def of_nat_ge_1_iff intro: Nat.gr0I)
lemma D'_sqr_less_D: "D'⇧2 < D"
proof -
have "D' ≤ sqrt D" by (auto simp: D'_def)
hence "real D' ^ 2 ≤ sqrt D ^ 2" by (intro power_mono) auto
also have "… = D" by simp
finally have "D'⇧2 ≤ D" by simp
moreover from nonsquare have "D ≠ D'⇧2" by auto
ultimately show ?thesis by simp
qed
lemma red_assoc_imp_irrat:
assumes "red_assoc pq"
shows "surd_to_real pq ∉ ℚ"
proof
assume rat: "surd_to_real pq ∈ ℚ"
with assms rat show False using irrat_sqrt_nonsquare[OF nonsquare]
by (auto simp: field_simps red_assoc_def surd_to_real_def divide_in_Rats_iff2 add_in_Rats_iff1)
qed
lemma surd_to_real_cnj_irrat:
assumes "red_assoc pq"
shows "surd_to_real_cnj pq ∉ ℚ"
proof
assume rat: "surd_to_real_cnj pq ∈ ℚ"
with assms rat show False using irrat_sqrt_nonsquare[OF nonsquare]
by (auto simp: field_simps red_assoc_def surd_to_real_cnj_def divide_in_Rats_iff2 diff_in_Rats_iff1)
qed
lemma surd_to_real_nonneg [intro]: "surd_to_real pq ≥ 0"
by (auto simp: surd_to_real_def case_prod_unfold divide_simps intro!: divide_nonneg_nonneg)
lemma surd_to_real_pos [intro]: "red_assoc pq ⟹ surd_to_real pq > 0"
by (auto simp: surd_to_real_def case_prod_unfold divide_simps red_assoc_def
intro!: divide_nonneg_nonneg)
lemma surd_to_real_nz [simp]: "red_assoc pq ⟹ surd_to_real pq ≠ 0"
by (auto simp: surd_to_real_def case_prod_unfold divide_simps red_assoc_def
intro!: divide_nonneg_nonneg)
lemma surd_to_real_cnj_nz [simp]: "red_assoc pq ⟹ surd_to_real_cnj pq ≠ 0"
using surd_to_real_cnj_irrat[of pq] by auto
lemma red_assoc_step:
assumes "red_assoc pq"
defines "X ≡ (D' + fst pq) div snd pq"
defines "pq' ≡ sqrt_remainder_step pq"
shows "red_assoc pq'"
"surd_to_real pq' = 1 / frac (surd_to_real pq)"
"surd_to_real_cnj pq' = 1 / (surd_to_real_cnj pq - X)"
"X > 0" "X * snd pq ≤ 2 * D'" "X = nat ⌊surd_to_real pq⌋"
"X = nat ⌊-1 / surd_to_real_cnj pq'⌋"
proof -
obtain p q where [simp]: "pq = (p, q)" by (cases pq)
obtain p' q' where [simp]: "pq' = (p', q')" by (cases pq')
define α where "α = (sqrt D + p) / q"
define α' where "α' = 1 / frac α"
define cnj_α' where "cnj_α' = (-sqrt D + (X * q - int p)) / ((D - (X * q - int p)⇧2) div q)"
from assms(1) have "α > 0" "q > 0"
by (auto simp: α_def red_assoc_def)
from assms(1) nonsquare have "α ∉ ℚ"
by (auto simp: α_def red_assoc_def divide_in_Rats_iff2 add_in_Rats_iff2 irrat_sqrt_nonsquare)
hence α'_pos: "frac α > 0" using Ints_subset_Rats by auto
from ‹pq' = (p', q')› have p'_def: "p' = X * q - p" and q'_def: "q' = (D - p'⇧2) div q"
unfolding pq'_def sqrt_remainder_step_def X_def by (auto simp: Let_def add_ac)
have "D' + p = ⌊sqrt D + p⌋"
by (auto simp: D'_def)
also have "… div int q = ⌊(sqrt D + p) / q⌋"
by (subst floor_divide_real_eq_div [symmetric]) auto
finally have X_altdef: "X = nat ⌊(sqrt D + p) / q⌋"
unfolding X_def zdiv_int [symmetric] by auto
have nz: "sqrt (real D) + (X * q - real p) ≠ 0"
proof
assume "sqrt (real D) + (X * q - real p) = 0"
hence "sqrt (real D) = real p - X * q" by (simp add: algebra_simps)
also have "… ∈ ℚ" by auto
finally show False using irrat_sqrt_nonsquare nonsquare by blast
qed
from assms(1) have "real (p ^ 2) ≤ sqrt D ^ 2"
unfolding of_nat_power by (intro power_mono) (auto simp: red_assoc_def field_simps)
also have "sqrt D ^ 2 = D" by simp
finally have "p⇧2 ≤ D" by (subst (asm) of_nat_le_iff)
have "frac α = α - X"
by (simp add: X_altdef frac_def α_def)
also have "… = (sqrt D - (X * q - int p)) / q"
using ‹q > 0› by (simp add: field_simps α_def)
finally have "1 / frac α = q / (sqrt D - (X * q - int p))"
by simp
also have "… = q * (sqrt D + (X * q - int p)) /
((sqrt D - (X * q - int p)) * (sqrt D + (X * q - int p)))" (is "_ = ?A / ?B")
using nz by (subst mult_divide_mult_cancel_right) auto
also have "?B = real_of_int (D - int p ^ 2 + 2 * X * p * q - int X ^ 2 * q ^ 2)"
by (auto simp: algebra_simps power2_eq_square)
also have "q dvd (D - p ^ 2)" using assms(1) by (auto simp: red_assoc_def)
with ‹p⇧2 ≤ D› have "int q dvd (int D - int p ^ 2)"
unfolding of_nat_power [symmetric] by (subst of_nat_diff [symmetric]) auto
hence "D - int p ^ 2 + 2 * X * p * q - int X ^ 2 * q ^ 2 = q * ((D - (X * q - int p)⇧2) div q)"
by (auto simp: power2_eq_square algebra_simps)
also have "?A / … = (sqrt D + (X * q - int p)) / ((D - (X * q - int p)⇧2) div q)"
unfolding of_int_mult of_int_of_nat_eq
by (rule mult_divide_mult_cancel_left) (insert ‹q > 0›, auto)
finally have α': "α' = …" by (simp add: α'_def)
have dvd: "q dvd (D - (X * q - int p)⇧2)"
using assms(1) ‹int q dvd (int D - int p ^ 2)›
by (auto simp: power2_eq_square algebra_simps)
have "X ≤ (sqrt D + p) / q" unfolding X_altdef by simp
moreover have "X ≠ (sqrt D + p) / q"
proof
assume "X = (sqrt D + p) / q"
hence "sqrt D = q * X - real p" using ‹q > 0› by (auto simp: field_simps)
also have "… ∈ ℚ" by auto
finally show False using irrat_sqrt_nonsquare[OF nonsquare] by simp
qed
ultimately have "X < (sqrt D + p) / q" by simp
hence *: "(X * q - int p) < sqrt D"
using ‹q > 0› by (simp add: field_simps)
moreover
have pos: "real_of_int (int D - (int X * int q - int p)⇧2) > 0"
proof (cases "X * q ≥ p")
case True
hence "real p ≤ real X * real q" unfolding of_nat_mult [symmetric] of_nat_le_iff .
hence "real_of_int ((X * q - int p) ^ 2) < sqrt D ^ 2" using *
unfolding of_int_power by (intro power_strict_mono) auto
also have "… = D" by simp
finally show ?thesis by simp
next
case False
hence less: "real X * real q < real p"
unfolding of_nat_mult [symmetric] of_nat_less_iff by auto
have "(real X * real q - real p)⇧2 = (real p - real X * real q)⇧2"
by (simp add: power2_eq_square algebra_simps)
also have "… ≤ real p ^ 2" using less by (intro power_mono) auto
also have "… < sqrt D ^ 2"
using ‹q > 0› assms(1) unfolding of_int_power
by (intro power_strict_mono) (auto simp: red_assoc_def field_simps)
also have "… = D" by simp
finally show ?thesis by simp
qed
hence pos': "int D - (int X * int q - int p)⇧2 > 0"
by (subst (asm) of_int_0_less_iff)
from pos have "real_of_int ((int D - (int X * int q - int p)⇧2) div q) > 0"
using ‹q > 0› dvd by (subst real_of_int_div) (auto intro!: divide_pos_pos)
ultimately have cnj_neg: "cnj_α' < 0" unfolding cnj_α'_def using dvd
unfolding of_int_0_less_iff by (intro divide_neg_pos) auto
have "(p - sqrt D) / q < 0"
using assms(1) by (auto simp: red_assoc_def X_altdef le_nat_iff)
also have "X ≥ 1"
using assms(1) by (auto simp: red_assoc_def X_altdef le_nat_iff)
hence "0 ≤ real X - 1" by simp
finally have "q < sqrt D + int q * X - p"
using ‹q > 0› by (simp add: field_simps)
hence "q * (sqrt D - (int q * X - p)) < (sqrt D + (int q * X - p)) * (sqrt D - (int q * X - p))"
using * by (intro mult_strict_right_mono) (auto simp: red_assoc_def X_altdef field_simps)
also have "… = D - (int q * X - p) ^ 2"
by (simp add: power2_eq_square algebra_simps)
finally have "cnj_α' > -1"
using dvd pos ‹q > 0› by (simp add: real_of_int_div field_simps cnj_α'_def)
from cnj_neg and this have "cnj_α' ∈ {-1<..<0}" by auto
have "α' > 1" using ‹frac α > 0›
by (auto simp: α'_def field_simps frac_lt_1)
have "0 = 1 + (-1 :: real)"
by simp
also have "1 + -1 < α' + cnj_α'"
using ‹cnj_α' > -1› and ‹α' > 1› by (intro add_strict_mono)
also have "α' + cnj_α' = 2 * (real X * q - real p) / ((int D - (int X * q - int p)⇧2) div int q)"
by (simp add: α' cnj_α'_def add_divide_distrib [symmetric])
finally have "real X * q - real p > 0" using pos dvd ‹q > 0›
by (subst (asm) zero_less_divide_iff, subst (asm) (1 2 3) real_of_int_div)
(auto simp: field_simps)
hence "real (X * q) > real p" unfolding of_nat_mult by simp
hence p_less_Xq: "p < X * q" by (simp only: of_nat_less_iff)
from pos' and p_less_Xq have "int D > int ((X * q - p)⇧2)"
by (subst of_nat_power) (auto simp: of_nat_diff)
hence pos'': "D > (X * q - p)⇧2" unfolding of_nat_less_iff .
from dvd have "int q dvd int (D - (X * q - p)⇧2)"
using p_less_Xq pos'' by (subst of_nat_diff) (auto simp: of_nat_diff)
with dvd have dvd': "q dvd (D - (X * q - p)⇧2)"
by simp
have α'_altdef: "α' = (sqrt D + p') / q'"
using dvd dvd' pos'' p_less_Xq α'
by (simp add: real_of_int_div p'_def q'_def real_of_nat_div mult_ac of_nat_diff)
have cnj_α'_altdef: "cnj_α' = (-sqrt D + p') / q'"
using dvd dvd' pos'' p_less_Xq unfolding cnj_α'_def
by (simp add: real_of_int_div p'_def q'_def real_of_nat_div mult_ac of_nat_diff)
from dvd' have dvd'': "q' dvd (D - p'⇧2)"
by (auto simp: mult_ac p'_def q'_def)
have "real ((D - p'⇧2) div q) > 0" unfolding p'_def
by (subst real_of_nat_div[OF dvd'], rule divide_pos_pos) (insert ‹q > 0› pos'', auto)
hence "q' > 0" unfolding q'_def of_nat_0_less_iff .
show "red_assoc pq'" using ‹α' > 1› and ‹cnj_α' ∈ _› and dvd'' and ‹q' > 0›
by (auto simp: red_assoc_def α'_altdef cnj_α'_altdef)
from assms(1) have "real p < sqrt D"
by (auto simp add: field_simps red_assoc_def)
hence "p ≤ D'" unfolding D'_def by linarith
with * have "real (X * q) < sqrt (real D) + D'"
by simp
thus "X * snd pq ≤ 2 * D'" unfolding D'_def ‹pq = (p, q)› snd_conv by linarith
have "(sqrt D + p') / q' = α'"
by (rule α'_altdef [symmetric])
also have "α' = 1 / frac ((sqrt D + p) / q)"
by (simp add: α'_def α_def)
finally show "surd_to_real pq' = 1 / frac (surd_to_real pq)" by (simp add: surd_to_real_def)
from ‹X ≥ 1› show "X > 0" by simp
from X_altdef show "X = nat ⌊surd_to_real pq⌋" by (simp add: surd_to_real_def)
have "sqrt (real D) < real p + 1 * real q"
using assms(1) by (auto simp: red_assoc_def field_simps)
also have "… ≤ real p + real X * real q"
using ‹X > 0› by (intro add_left_mono mult_right_mono) (auto simp: of_nat_ge_1_iff)
finally have "sqrt (real D) < …" .
have "real p < sqrt D"
using assms(1) by (auto simp add: field_simps red_assoc_def)
also have "… ≤ sqrt D + q * X"
by linarith
finally have less: "real p < sqrt D + X * q" by (simp add: algebra_simps)
moreover have "D + p * p' + X * q * sqrt D = q * q' + p * sqrt D + p' * sqrt D + X * p' * q"
using dvd' pos'' p_less_Xq ‹q > 0› unfolding p'_def q'_def of_nat_mult of_nat_add
by (simp add: power2_eq_square field_simps of_nat_diff real_of_nat_div)
ultimately show *: "surd_to_real_cnj pq' = 1 / (surd_to_real_cnj pq - X)"
using ‹q > 0› ‹q' > 0› by (auto simp: surd_to_real_cnj_def field_simps)
have **: "a = nat ⌊y⌋" if "x ≥ 0" "x < 1" "real a + x = y" for a :: nat and x y :: real
using that by linarith
from assms(1) have surd_to_real_cnj: "surd_to_real_cnj (p, q) ∈ {-1<..<0}"
by (auto simp: surd_to_real_cnj_def red_assoc_def)
have "surd_to_real_cnj (p, q) < X"
using assms(1) less by (auto simp: surd_to_real_cnj_def field_simps red_assoc_def)
hence "real X = surd_to_real_cnj (p, q) - 1 / surd_to_real_cnj (p', q')" using *
using surd_to_real_cnj_irrat assms(1) ‹red_assoc pq'› by (auto simp: field_simps)
thus "X = nat ⌊-1 / surd_to_real_cnj pq'⌋" using surd_to_real_cnj
by (intro **[of "-surd_to_real_cnj (p, q)"]) auto
qed
lemma red_assoc_denom_2D:
assumes "red_assoc (p, q)"
defines "X ≡ (D' + p) div q"
assumes "X > D'"
shows "q = 1"
proof -
have "X * q ≤ 2 * D'" "X > 0"
using red_assoc_step(4,5)[OF assms(1)] by (simp_all add: X_def)
note this(1)
also have "2 * D' < 2 * X"
by (intro mult_strict_left_mono assms) auto
finally have "q < 2" using ‹X > 0› by simp
moreover from assms(1) have "q > 0" by (auto simp: red_assoc_def)
ultimately show ?thesis by simp
qed
lemma red_assoc_denom_1:
assumes "red_assoc (p, 1)"
shows "p = D'"
proof -
from assms have "sqrt D > p" "sqrt D < real p + 1"
by (auto simp: red_assoc_def)
thus "p = D'" unfolding D'_def
by linarith
qed
lemma red_assoc_begin:
"red_assoc (D', D - D'⇧2)"
"surd_to_real (D', D - D'⇧2) = 1 / frac (sqrt D)"
"surd_to_real_cnj (D', D - D'⇧2) = -1 / (sqrt D + D')"
proof -
have pos: "D > 0" "D' > 0"
using nonsquare by (auto simp: D'_def of_nat_ge_1_iff intro!: Nat.gr0I)
have "sqrt D ≠ D'"
using irrat_sqrt_nonsquare[OF nonsquare] by auto
moreover have "sqrt D ≥ 0" by simp
hence "D' ≤ sqrt D" unfolding D'_def by linarith
ultimately have less: "D' < sqrt D" by simp
have "sqrt D ≠ D' + 1"
using irrat_sqrt_nonsquare[OF nonsquare] by auto
moreover have "sqrt D ≥ 0" by simp
hence "D' ≥ sqrt D - 1" unfolding D'_def by linarith
ultimately have gt: "D' > sqrt D - 1" by simp
from less have "real D' ^ 2 < sqrt D ^ 2" by (intro power_strict_mono) auto
also have "… = D" by simp
finally have less': "D'⇧2 < D" unfolding of_nat_power [symmetric] of_nat_less_iff .
moreover have "real D' * (real D' - 1) < sqrt D * (sqrt D - 1)"
using less pos
by (intro mult_strict_mono diff_strict_right_mono) (auto simp: of_nat_ge_1_iff)
hence "D'⇧2 + sqrt D < D' + D"
by (simp add: field_simps power2_eq_square)
moreover have "(sqrt D - 1) * sqrt D < real D' * (real D' + 1)"
using pos gt by (intro mult_strict_mono) auto
hence "D < sqrt D + D'⇧2 + D'" by (simp add: power2_eq_square field_simps)
ultimately show "red_assoc (D', D - D'⇧2)"
by (auto simp: red_assoc_def field_simps of_nat_diff less)
have frac: "frac (sqrt D) = sqrt D - D'" unfolding frac_def D'_def
by auto
show "surd_to_real (D', D - D'⇧2) = 1 / frac (sqrt D)" unfolding surd_to_real_def
using less less' pos by (subst frac) (auto simp: of_nat_diff power2_eq_square field_simps)
have "surd_to_real_cnj (D', D - D'⇧2) = -((sqrt D - D') / (D - D'⇧2))"
using less less' pos by (auto simp: surd_to_real_cnj_def field_simps)
also have "real (D - D'⇧2) = (sqrt D - D') * (sqrt D + D')"
using less' by (simp add: power2_eq_square algebra_simps of_nat_diff)
also have "(sqrt D - D') / … = 1 / (sqrt D + D')"
using less by (subst nonzero_divide_mult_cancel_left) auto
finally show "surd_to_real_cnj (D', D - D'⇧2) = -1 / (sqrt D + D')" by simp
qed
lemma cfrac_remainder_surd_to_real:
assumes "red_assoc pq"
shows "cfrac_remainder (cfrac_of_real (surd_to_real pq)) n =
surd_to_real ((sqrt_remainder_step ^^ n) pq)"
using assms(1)
proof (induction n arbitrary: pq)
case 0
hence "cfrac_lim (cfrac_of_real (surd_to_real pq)) = surd_to_real pq"
by (intro cfrac_lim_of_real red_assoc_imp_irrat 0)
thus ?case using 0
by auto
next
case (Suc n)
obtain p q where [simp]: "pq = (p, q)" by (cases pq)
have "surd_to_real ((sqrt_remainder_step ^^ Suc n) pq) =
surd_to_real ((sqrt_remainder_step ^^ n) (sqrt_remainder_step (p, q)))"
by (subst funpow_Suc_right) auto
also have "… = cfrac_remainder (cfrac_of_real (surd_to_real (sqrt_remainder_step (p, q)))) n"
using red_assoc_step(1)[of "(p, q)"] Suc.prems
by (intro Suc.IH [symmetric]) (auto simp: sqrt_remainder_step_def Let_def add_ac)
also have "surd_to_real (sqrt_remainder_step (p, q)) = 1 / frac (surd_to_real (p, q))"
using red_assoc_step(2)[of "(p, q)"] Suc.prems
by (auto simp: sqrt_remainder_step_def Let_def add_ac surd_to_real_def)
also have "cfrac_of_real … = cfrac_tl (cfrac_of_real (surd_to_real (p, q)))"
using Suc.prems Ints_subset_Rats red_assoc_imp_irrat by (subst cfrac_tl_of_real) auto
also have "cfrac_remainder … n = cfrac_remainder (cfrac_of_real (surd_to_real (p, q))) (Suc n)"
by (simp add: cfrac_drop_Suc_right cfrac_remainder_def)
finally show ?case by simp
qed
lemma red_assoc_step' [intro]: "red_assoc pq ⟹ red_assoc (sqrt_remainder_step pq)"
using red_assoc_step(1)[of pq]
by (simp add: sqrt_remainder_step_def case_prod_unfold add_ac Let_def)
lemma red_assoc_steps [intro]: "red_assoc pq ⟹ red_assoc ((sqrt_remainder_step ^^ n) pq)"
by (induction n) auto
lemma floor_sqrt_less_sqrt: "D' < sqrt D"
proof -
have "D' ≤ sqrt D" unfolding D'_def by auto
moreover have "sqrt D ≠ D'"
using irrat_sqrt_nonsquare[OF nonsquare] by auto
ultimately show ?thesis by auto
qed
lemma red_assoc_bounds:
assumes "red_assoc pq"
shows "pq ∈ (SIGMA p:{0<..D'}. {Suc D' - p..D' + p})"
proof -
obtain p q where [simp]: "pq = (p, q)" by (cases pq)
from assms have *: "p < sqrt D"
by (auto simp: red_assoc_def field_simps)
hence p: "p ≤ D'" unfolding D'_def by linarith
from assms have "p > 0" by (auto intro!: Nat.gr0I simp: red_assoc_def)
have "q > sqrt D - p" "q < sqrt D + p"
using assms by (auto simp: red_assoc_def field_simps)
hence "q ≥ D' + 1 - p" "q ≤ D' + p"
unfolding D'_def by linarith+
with p ‹p > 0› show ?thesis by simp
qed
lemma surd_to_real_cnj_eq_iff:
assumes "red_assoc pq" "red_assoc pq'"
shows "surd_to_real_cnj pq = surd_to_real_cnj pq' ⟷ pq = pq'"
proof
assume eq: "surd_to_real_cnj pq = surd_to_real_cnj pq'"
from assms have pos: "snd pq > 0" "snd pq' > 0" by (auto simp: red_assoc_def)
have "snd pq = snd pq'"
proof (rule ccontr)
assume "snd pq ≠ snd pq'"
with eq have "sqrt D = (real (fst pq' * snd pq) - fst pq * snd pq') / (real (snd pq) - snd pq')"
using pos by (auto simp: field_simps surd_to_real_cnj_def case_prod_unfold)
also have "… ∈ ℚ" by auto
finally show False using irrat_sqrt_nonsquare[OF nonsquare] by auto
qed
moreover from this eq pos have "fst pq = fst pq'"
by (auto simp: surd_to_real_cnj_def case_prod_unfold)
ultimately show "pq = pq'" by (simp add: prod_eq_iff)
qed auto
lemma red_assoc_sqrt_remainder_surd [intro]: "red_assoc (sqrt_remainder_surd n)"
by (auto simp: sqrt_remainder_surd_def intro!: red_assoc_begin)
lemma surd_to_real_sqrt_remainder_surd:
"surd_to_real (sqrt_remainder_surd n) = cfrac_remainder (cfrac_of_real (sqrt D)) (Suc n)"
proof (induction n)
case 0
from nonsquare have "D > 0" by (auto intro!: Nat.gr0I)
with red_assoc_begin show ?case using nonsquare irrat_sqrt_nonsquare[OF nonsquare]
using Ints_subset_Rats cfrac_drop_Suc_right cfrac_remainder_def cfrac_tl_of_real
sqrt_remainder_surd_def by fastforce
next
case (Suc n)
have "surd_to_real (sqrt_remainder_surd (Suc n)) =
surd_to_real (sqrt_remainder_step (sqrt_remainder_surd n))"
by (simp add: sqrt_remainder_surd_def)
also have "… = 1 / frac (surd_to_real (sqrt_remainder_surd n))"
using red_assoc_step[OF red_assoc_sqrt_remainder_surd[of n]] by simp
also have "surd_to_real (sqrt_remainder_surd n) =
cfrac_remainder (cfrac_of_real (sqrt D)) (Suc n)" (is "_ = ?X")
by (rule Suc.IH)
also have "⌊cfrac_remainder (cfrac_of_real (sqrt (real D))) (Suc n)⌋ =
cfrac_nth (cfrac_of_real (sqrt (real D))) (Suc n)"
using irrat_sqrt_nonsquare[OF nonsquare] by (intro floor_cfrac_remainder) auto
hence "1 / frac ?X = cfrac_remainder (cfrac_of_real (sqrt D)) (Suc (Suc n))"
using irrat_sqrt_nonsquare[OF nonsquare]
by (subst cfrac_remainder_Suc[of "Suc n"])
(simp_all add: frac_def cfrac_length_of_real_irrational)
finally show ?case .
qed
lemma sqrt_cfrac: "sqrt_cfrac_nth n = cfrac_nth (cfrac_of_real (sqrt D)) (Suc n)"
proof -
have "cfrac_nth (cfrac_of_real (sqrt D)) (Suc n) =
⌊cfrac_remainder (cfrac_of_real (sqrt D)) (Suc n)⌋"
using irrat_sqrt_nonsquare[OF nonsquare] by (subst floor_cfrac_remainder) auto
also have "cfrac_remainder (cfrac_of_real (sqrt D)) (Suc n) = surd_to_real (sqrt_remainder_surd n)"
by (rule surd_to_real_sqrt_remainder_surd [symmetric])
also have "nat ⌊surd_to_real (sqrt_remainder_surd n)⌋ = sqrt_cfrac_nth n"
unfolding sqrt_cfrac_nth_def using red_assoc_step(6)[OF red_assoc_sqrt_remainder_surd[of n]]
by (simp add: case_prod_unfold)
finally show ?thesis
by (simp add: nat_eq_iff)
qed
lemma sqrt_cfrac_pos: "sqrt_cfrac_nth k > 0"
using red_assoc_step(4)[OF red_assoc_sqrt_remainder_surd[of k]]
by (simp add: sqrt_cfrac_nth_def case_prod_unfold)
lemma snd_sqrt_remainder_surd_pos: "snd (sqrt_remainder_surd n) > 0"
using red_assoc_sqrt_remainder_surd[of n] by (auto simp: red_assoc_def)
lemma
shows period_nonempty: "l > 0"
and period_length_le_aux: "l ≤ D' * (D' + 1)"
and sqrt_remainder_surd_periodic: "⋀n. sqrt_remainder_surd n = sqrt_remainder_surd (n mod l)"
and sqrt_cfrac_periodic: "⋀n. sqrt_cfrac_nth n = sqrt_cfrac_nth (n mod l)"
and sqrt_remainder_surd_smallest_period:
"⋀n. n ∈ {0<..<l} ⟹ sqrt_remainder_surd n ≠ sqrt_remainder_surd 0"
and snd_sqrt_remainder_surd_gt_1: "⋀n. n < l - 1 ⟹ snd (sqrt_remainder_surd n) > 1"
and sqrt_cfrac_le: "⋀n. n < l - 1 ⟹ sqrt_cfrac_nth n ≤ D'"
and sqrt_remainder_surd_last: "sqrt_remainder_surd (l - 1) = (D', 1)"
and sqrt_cfrac_last: "sqrt_cfrac_nth (l - 1) = 2 * D'"
and sqrt_cfrac_palindrome: "⋀n. n < l - 1 ⟹ sqrt_cfrac_nth (l - n - 2) = sqrt_cfrac_nth n"
and sqrt_cfrac_smallest_period:
"⋀l'. l' > 0 ⟹ (⋀k. sqrt_cfrac_nth (k + l') = sqrt_cfrac_nth k) ⟹ l' ≥ l"
proof -
note [simp] = sqrt_remainder_surd_def
define f where "f = sqrt_remainder_surd"
have *[intro]: "red_assoc (f n)" for n
unfolding f_def by (rule red_assoc_sqrt_remainder_surd)
define S where "S = (SIGMA p:{0<..D'}. {Suc D' - p..D' + p})"
have [intro]: "finite S" by (simp add: S_def)
have "card S = (∑p=1..D'. 2 * p)" unfolding S_def
by (subst card_SigmaI) (auto intro!: sum.cong)
also have "… = D' * (D' + 1)"
by (induction D') (auto simp: power2_eq_square)
finally have [simp]: "card S = D' * (D' + 1)" .
have "D' * (D' + 1) + 1 = card {..D' * (D' + 1)}" by simp
define k1 where
"k1 = (LEAST k1. k1 ≤ D' * (D' + 1) ∧ (∃k2. k2 ≤ D' * (D' + 1) ∧ k1 ≠ k2 ∧ f k1 = f k2))"
define k2 where
"k2 = (LEAST k2. k2 ≤ D' * (D' + 1) ∧ k1 ≠ k2 ∧ f k1 = f k2)"
have "f ` {..D' * (D' + 1)} ⊆ S" unfolding S_def
using red_assoc_bounds[OF *] by blast
hence "card (f ` {..D' * (D' + 1)}) ≤ card S"
by (intro card_mono) auto
also have "card S = D' * (D' + 1)" by simp
also have "… < card {..D' * (D' + 1)}" by simp
finally have "¬inj_on f {..D' * (D' + 1)}"
by (rule pigeonhole)
hence "∃k1. k1 ≤ D' * (D' + 1) ∧ (∃k2. k2 ≤ D' * (D' + 1) ∧ k1 ≠ k2 ∧ f k1 = f k2)"
by (auto simp: inj_on_def)
from LeastI_ex[OF this, folded k1_def]
have "k1 ≤ D' * (D' + 1)" "∃k2≤D' * (D' + 1). k1 ≠ k2 ∧ f k1 = f k2" by auto
moreover from LeastI_ex[OF this(2), folded k2_def]
have "k2 ≤ D' * (D' + 1)" "k1 ≠ k2" "f k1 = f k2" by auto
moreover have "k1 ≤ k2"
proof (rule ccontr)
assume "¬(k1 ≤ k2)"
hence "k2 ≤ D' * (D' + 1) ∧ (∃k2'. k2' ≤ D' * (D' + 1) ∧ k2 ≠ k2' ∧ f k2 = f k2')"
using ‹k1 ≤ D' * (D' + 1)› and ‹k1 ≠ k2› and ‹f k1 = f k2› by auto
hence "k1 ≤ k2" unfolding k1_def by (rule Least_le)
with ‹¬(k1 ≤ k2)› show False by simp
qed
ultimately have k12: "k1 < k2" "k2 ≤ D' * (D' + 1)" "f k1 = f k2" by auto
have [simp]: "k1 = 0"
proof (cases k1)
case (Suc k1')
define k2' where "k2' = k2 - 1"
have Suc': "k2 = Suc k2'" using k12 by (simp add: k2'_def)
have nz: "surd_to_real_cnj (sqrt_remainder_step (f k1')) ≠ 0"
"surd_to_real_cnj (sqrt_remainder_step (f k2')) ≠ 0"
using surd_to_real_cnj_nz[OF *[of k2]] surd_to_real_cnj_nz[OF *[of k1]]
by (simp_all add: f_def Suc Suc')
define a where "a = (D' + fst (f k1)) div snd (f k1)"
define a' where "a' = (D' + fst (f k1')) div snd (f k1')"
define a'' where "a'' = (D' + fst (f k2')) div snd (f k2')"
have "a' = nat ⌊- 1 / surd_to_real_cnj (sqrt_remainder_step (f k1'))⌋"
using red_assoc_step[OF *[of k1']] by (simp add: a'_def)
also have "sqrt_remainder_step (f k1') = f k1"
by (simp add: Suc f_def)
also have "f k1 = f k2" by fact
also have "f k2 = sqrt_remainder_step (f k2')" by (simp add: Suc' f_def)
also have "nat ⌊- 1 / surd_to_real_cnj (sqrt_remainder_step (f k2'))⌋ = a''"
using red_assoc_step[OF *[of k2']] by (simp add: a''_def)
finally have a'_a'': "a' = a''" .
have "surd_to_real_cnj (f k2') ≠ a''"
using surd_to_real_cnj_irrat[OF *[of k2']] by auto
hence "surd_to_real_cnj (f k2') = 1 / surd_to_real_cnj (sqrt_remainder_step (f k2')) + a''"
using red_assoc_step(3)[OF *[of k2'], folded a''_def] nz
by (simp add: field_simps)
also have "… = 1 / surd_to_real_cnj (sqrt_remainder_step (f k1')) + a'"
using k12 by (simp add: a'_a'' k12 Suc Suc' f_def)
also have nz': "surd_to_real_cnj (f k1') ≠ a'"
using surd_to_real_cnj_irrat[OF *[of k1']] by auto
hence "1 / surd_to_real_cnj (sqrt_remainder_step (f k1')) + a' = surd_to_real_cnj (f k1')"
using red_assoc_step(3)[OF *[of k1'], folded a'_def] nz nz'
by (simp add: field_simps)
finally have "f k1' = f k2'"
by (subst (asm) surd_to_real_cnj_eq_iff) auto
with k12 have "k1' ≤ D' * (D' + 1) ∧ (∃k2≤D' * (D' + 1). k1' ≠ k2 ∧ f k1' = f k2)"
by (auto simp: Suc Suc' intro!: exI[of _ k2'])
hence "k1 ≤ k1'" unfolding k1_def by (rule Least_le)
thus "k1 = 0" by (simp add: Suc)
qed auto
have smallest_period: "f k ≠ f 0" if "k ∈ {0<..<k2}" for k
proof
assume "f k = f 0"
hence "k ≤ D' * (D' + 1) ∧ k1 ≠ k ∧ f k1 = f k"
using k12 that by auto
hence "k2 ≤ k" unfolding k2_def by (rule Least_le)
with that show False by auto
qed
have snd_f_gt_1: "snd (f k) > 1" if "k < k2 - 1" for k
proof -
have "snd (f k) ≠ 1"
proof
assume "snd (f k) = 1"
hence "f k = (D', 1)" using red_assoc_denom_1[of "fst (f k)"] *[of k]
by (cases "f k") auto
hence "sqrt_remainder_step (f k) = (D', D - D'⇧2)" by (auto simp: sqrt_remainder_step_def)
hence "f (Suc k) = f 0" by (simp add: f_def)
moreover have "f (Suc k) ≠ f 0"
using that by (intro smallest_period) auto
ultimately show False by contradiction
qed
moreover have "snd (f k) > 0" using *[of k] by (auto simp: red_assoc_def)
ultimately show ?thesis by simp
qed
have sqrt_cfrac_le: "sqrt_cfrac_nth k ≤ D'" if "k < k2 - 1" for k
proof -
define p and q where "p = fst (f k)" and "q = snd (f k)"
have "q ≥ 2" using snd_f_gt_1[of k] that by (auto simp: q_def)
also have "sqrt_cfrac_nth k * q ≤ D' * 2"
using red_assoc_step(5)[OF *[of k]]
by (simp add: sqrt_cfrac_nth_def p_def q_def case_prod_unfold f_def)
finally show ?thesis by simp
qed
have last: "f (k2 - 1) = (D', 1)"
proof -
define p and q where "p = fst (f (k2 - 1))" and "q = snd (f (k2 - 1))"
have pq: "f (k2 - 1) = (p, q)" by (simp add: p_def q_def)
have "sqrt_remainder_step (f (k2 - 1)) = f (Suc (k2 - 1))"
by (simp add: f_def)
also from k12 have "Suc (k2 - 1) = k2" by simp
also have "f k2 = f 0"
using k12 by simp
also have "f 0 = (D', D - D'⇧2)" by (simp add: f_def)
finally have eq: "sqrt_remainder_step (f (k2 - 1)) = (D', D - D'⇧2)" .
hence "(D - D'⇧2) div q = D - D'⇧2" unfolding sqrt_remainder_step_def Let_def pq
by auto
moreover have "q > 0" using *[of "k2 - 1"]
by (auto simp: red_assoc_def q_def)
ultimately have "q = 1" using D'_sqr_less_D
by (subst (asm) div_eq_dividend_iff) auto
hence "p = D'"
using red_assoc_denom_1[of p] *[of "k2 - 1"] unfolding pq by auto
with ‹q = 1› show "f (k2 - 1) = (D', 1)" unfolding pq by simp
qed
have period: "sqrt_remainder_surd n = sqrt_remainder_surd (n mod k2)" for n
unfolding sqrt_remainder_surd_def using k12
by (metis ‹k1 = 0› f_def funpow_mod_eq funpow_0 sqrt_remainder_surd_def)
have period': "sqrt_cfrac_nth k = sqrt_cfrac_nth (k mod k2)" for k
using period[of k] by (simp add: sqrt_cfrac_nth_def)
have k2_le: "l ≥ k2" if "l > 0" "⋀k. sqrt_cfrac_nth (k + l) = sqrt_cfrac_nth k" for l
proof (rule ccontr)
assume *: "¬(l ≥ k2)"
hence "sqrt_cfrac_nth (k2 - Suc l) = sqrt_cfrac_nth (k2 - 1)"
using that(2)[of "k2 - Suc l"] by simp
also have "… = 2 * D'"
using last by (simp add: sqrt_cfrac_nth_def f_def)
finally have "2 * D' = sqrt_cfrac_nth (k2 - Suc l)" ..
also have "… ≤ D'" using k12 that *
by (intro sqrt_cfrac_le diff_less_mono2) auto
finally show False using D'_pos by simp
qed
have "l = (LEAST l. 0 < l ∧ (∀n. int (sqrt_cfrac_nth (n + l)) = int (sqrt_cfrac_nth n)))"
using nonsquare unfolding sqrt_cfrac_def
by (simp add: l_def sqrt_nat_period_length_def sqrt_cfrac)
hence l_altdef: "l = (LEAST l. 0 < l ∧ (∀n. sqrt_cfrac_nth (n + l) = sqrt_cfrac_nth n))"
by simp
have [simp]: "D ≠ 0" using nonsquare by (auto intro!: Nat.gr0I)
have "∃l. l > 0 ∧ (∀k. sqrt_cfrac_nth (k + l) = sqrt_cfrac_nth k)"
proof (rule exI, safe)
fix k show "sqrt_cfrac_nth (k + k2) = sqrt_cfrac_nth k"
using period'[of k] period'[of "k + k2"] k12 by simp
qed (insert k12, auto)
from LeastI_ex[OF this, folded l_altdef]
have l: "l > 0" "⋀k. sqrt_cfrac_nth (k + l) = sqrt_cfrac_nth k"
by (simp_all add: sqrt_cfrac)
have "l ≤ k2" unfolding l_altdef
by (rule Least_le) (subst (1 2) period', insert k12, auto)
moreover have "k2 ≤ l" using k2_le l by blast
ultimately have [simp]: "l = k2" by auto
define x' where "x' = (λk. -1 / surd_to_real_cnj (f k))"
{
fix k :: nat
have nz: "surd_to_real_cnj (f k) ≠ 0" "surd_to_real_cnj (f (Suc k)) ≠ 0"
using surd_to_real_cnj_nz[OF *, of k] surd_to_real_cnj_nz[OF *, of "Suc k"]
by (simp_all add: f_def)
have "surd_to_real_cnj (f k) ≠ sqrt_cfrac_nth k"
using surd_to_real_cnj_irrat[OF *[of k]] by auto
hence "x' (Suc k) = sqrt_cfrac_nth k + 1 / x' k"
using red_assoc_step(3)[OF *[of k]] nz
by (simp add: field_simps sqrt_cfrac_nth_def case_prod_unfold f_def x'_def)
} note x'_Suc = this
have x'_nz: "x' k ≠ 0" for k
using surd_to_real_cnj_nz[OF *[of k]] by (auto simp: x'_def)
have x'_0: "x' 0 = real D' + sqrt D"
using red_assoc_begin by (simp add: x'_def f_def)
define c' where "c' = cfrac (λn. sqrt_cfrac_nth (l - Suc n))"
define c'' where "c'' = cfrac (λn. if n = 0 then 2 * D' else sqrt_cfrac_nth (n - 1))"
have nth_c' [simp]: "cfrac_nth c' n = sqrt_cfrac_nth (l - Suc n)" for n
unfolding c'_def by (subst cfrac_nth_cfrac) (auto simp: is_cfrac_def intro!: sqrt_cfrac_pos)
have nth_c'' [simp]: "cfrac_nth c'' n = (if n = 0 then 2 * D' else sqrt_cfrac_nth (n - 1))" for n
unfolding c''_def by (subst cfrac_nth_cfrac) (auto simp: is_cfrac_def intro!: sqrt_cfrac_pos)
have "conv' c' n (x' (l - n)) = x' l" if "n ≤ l" for n
using that
proof (induction n)
case (Suc n)
have "x' l = conv' c' n (x' (l - n))"
using Suc.prems by (intro Suc.IH [symmetric]) auto
also have "l - n = Suc (l - Suc n)"
using Suc.prems by simp
also have "x' … = cfrac_nth c' n + 1 / x' (l - Suc n)"
by (subst x'_Suc) simp
also have "conv' c' n … = conv' c' (Suc n) (x' (l - Suc n))"
by (simp add: conv'_Suc_right)
finally show ?case ..
qed simp_all
from this[of l] have conv'_x'_0: "conv' c' l (x' 0) = x' 0"
using k12 by (simp add: x'_def)
have "cfrac_nth (cfrac_of_real (x' 0)) n = cfrac_nth c'' n" for n
proof (cases n)
case 0
thus ?thesis by (simp add: x'_0 D'_def)
next
case (Suc n')
have "sqrt D ∉ ℤ"
using red_assoc_begin(1) red_assoc_begin(2) by auto
hence "cfrac_nth (cfrac_of_real (real D' + sqrt (real D))) (Suc n') =
cfrac_nth (cfrac_of_real (sqrt (real D))) (Suc n')"
by (simp add: cfrac_tl_of_real frac_add_of_nat Ints_add_left_cancel flip: cfrac_nth_tl)
thus ?thesis using x'_nz[of 0]
by (simp add: x'_0 sqrt_cfrac Suc)
qed
show "sqrt_cfrac_nth (l - n - 2) = sqrt_cfrac_nth n" if "n < l - 1" for n
proof -
have "D > 1" using nonsquare by (cases D) (auto intro!: Nat.gr0I)
hence "D' + sqrt D > 0 + 1" using D'_pos by (intro add_strict_mono) auto
hence "x' 0 > 1" by (auto simp: x'_0)
hence "cfrac_nth c' (Suc n) = cfrac_nth (cfrac_of_real (conv' c' l (x' 0))) (Suc n)"
using ‹n < l - 1› using cfrac_of_real_conv' by auto
also have "… = cfrac_nth (cfrac_of_real (x' 0)) (Suc n)"
by (subst conv'_x'_0) auto
also have "… = cfrac_nth c'' (Suc n)" by fact
finally show "sqrt_cfrac_nth (l - n - 2) = sqrt_cfrac_nth n"
by simp
qed
show "l > 0" "l ≤ D' * (D' + 1)" using k12 by simp_all
show "sqrt_remainder_surd n = sqrt_remainder_surd (n mod l)"
"sqrt_cfrac_nth n = sqrt_cfrac_nth (n mod l)" for n
using period[of n] period'[of n] by simp_all
show "sqrt_remainder_surd n ≠ sqrt_remainder_surd 0" if "n ∈ {0<..<l}" for n
using smallest_period[of n] that by (auto simp: f_def)
show "snd (sqrt_remainder_surd n) > 1" if "n < l - 1" for n
using that snd_f_gt_1[of n] by (simp add: f_def)
show "f (l - 1) = (D', 1)" and "sqrt_cfrac_nth (l - 1) = 2 * D'"
using last by (simp_all add: sqrt_cfrac_nth_def f_def)
show "sqrt_cfrac_nth k ≤ D'" if "k < l - 1" for k
using sqrt_cfrac_le[of k] that by simp
show "l' ≥ l" if "l' > 0" "⋀k. sqrt_cfrac_nth (k + l') = sqrt_cfrac_nth k" for l'
using k2_le[of l'] that by auto
qed
theorem cfrac_sqrt_periodic:
"cfrac_nth (cfrac_of_real (sqrt D)) (Suc n) =
cfrac_nth (cfrac_of_real (sqrt D)) (Suc (n mod l))"
using sqrt_cfrac_periodic[of n] by (metis sqrt_cfrac)
theorem cfrac_sqrt_le: "n ∈ {0<..<l} ⟹ cfrac_nth (cfrac_of_real (sqrt D)) n ≤ D'"
using sqrt_cfrac_le[of "n - 1"]
by (metis Suc_less_eq Suc_pred add.right_neutral greaterThanLessThan_iff of_nat_mono
period_nonempty plus_1_eq_Suc sqrt_cfrac)
theorem cfrac_sqrt_last: "cfrac_nth (cfrac_of_real (sqrt D)) l = 2 * D'"
using sqrt_cfrac_last by (metis One_nat_def Suc_pred period_nonempty sqrt_cfrac)
theorem cfrac_sqrt_palindrome:
assumes "n ∈ {0<..<l}"
shows "cfrac_nth (cfrac_of_real (sqrt D)) (l - n) = cfrac_nth (cfrac_of_real (sqrt D)) n"
proof -
have "cfrac_nth (cfrac_of_real (sqrt D)) (l - n) = sqrt_cfrac_nth (l - n - 1)"
using assms by (subst sqrt_cfrac) (auto simp: Suc_diff_Suc)
also have "… = sqrt_cfrac_nth (n - 1)"
using assms by (subst sqrt_cfrac_palindrome [symmetric]) auto
also have "… = cfrac_nth (cfrac_of_real (sqrt D)) n"
using assms by (subst sqrt_cfrac) auto
finally show ?thesis .
qed
lemma sqrt_cfrac_info_palindrome:
assumes "sqrt_cfrac_info D = (a, b, cs)"
shows "rev (butlast cs) = butlast cs"
proof (rule List.nth_equalityI; safe?)
fix i assume "i < length (rev (butlast cs))"
with period_nonempty have "Suc i < length cs" by simp
thus "rev (butlast cs) ! i = butlast cs ! i"
using assms cfrac_sqrt_palindrome[of "Suc i"] period_nonempty unfolding l_def
by (auto simp: sqrt_cfrac_info_def rev_nth algebra_simps Suc_diff_Suc simp del: cfrac.simps)
qed simp_all
lemma sqrt_cfrac_info_last:
assumes "sqrt_cfrac_info D = (a, b, cs)"
shows "last cs = 2 * Discrete.sqrt D"
proof -
from assms show ?thesis using period_nonempty cfrac_sqrt_last
by (auto simp: sqrt_cfrac_info_def last_map l_def D'_def Discrete_sqrt_altdef)
qed
text ‹
The following lemmas allow us to compute the period of the expansion of the square root:
›
lemma while_option_sqrt_cfrac:
defines "step' ≡ (λ(as, pq). ((D' + fst pq) div snd pq # as, sqrt_remainder_step pq))"
defines "b ≡ (λ(_, pq). snd pq ≠ 1)"
defines "initial ≡ ([] :: nat list, (D', D - D'⇧2))"
shows "while_option b step' initial =
Some (rev (map sqrt_cfrac_nth [0..<l -1]), (D', 1))"
proof -
define P where
"P = (λ(as, pq). let n = length as
in n < l ∧ pq = sqrt_remainder_surd n ∧ as = rev (map sqrt_cfrac_nth [0..<n]))"
define μ :: "nat list × (nat × nat) ⇒ nat" where "μ = (λ(as, _). l - length as)"
have [simp]: "P initial" using period_nonempty
by (auto simp: initial_def P_def sqrt_remainder_surd_def)
have step': "P (step' s) ∧ Suc (length (fst s)) < l" if "P s" "b s" for s
proof (cases s)
case (fields as p q)
define n where "n = length as"
from that fields sqrt_remainder_surd_last have "Suc n ≤ l"
by (auto simp: b_def P_def Let_def n_def [symmetric])
moreover from that fields sqrt_remainder_surd_last have "Suc n ≠ l"
by (auto simp: b_def P_def Let_def n_def [symmetric])
ultimately have "Suc n < l" by auto
with that fields sqrt_remainder_surd_last show "P (step' s) ∧ Suc (length (fst s)) < l"
by (simp add: b_def P_def Let_def n_def step'_def sqrt_cfrac_nth_def
sqrt_remainder_surd_def case_prod_unfold)
qed
have [simp]: "length (fst (step' s)) = Suc (length (fst s))" for s
by (simp add: step'_def case_prod_unfold)
have "∃x. while_option b step' initial = Some x"
proof (rule measure_while_option_Some)
fix s assume *: "P s" "b s"
from step'[OF *] show "P (step' s) ∧ μ (step' s) < μ s"
by (auto simp: b_def μ_def case_prod_unfold intro!: diff_less_mono2)
qed auto
then obtain x where x: "while_option b step' initial = Some x" ..
have "P x" by (rule while_option_rule[OF _ x]) (insert step', auto)
have "¬b x" using while_option_stop[OF x] by auto
obtain as p q where [simp]: "x = (as, (p, q))" by (cases x)
define n where "n = length as"
have [simp]: "q = 1" using ‹¬b x› by (auto simp: b_def)
have [simp]: "p = D'" using ‹P x›
using red_assoc_denom_1[of p] by (auto simp: P_def Let_def)
have "n < l" "sqrt_remainder_surd (length as) = (D', Suc 0)"
and as: "as = rev (map sqrt_cfrac_nth [0..<n])" using ‹P x›
by (auto simp: P_def Let_def n_def)
hence "¬(n < l - 1)"
using snd_sqrt_remainder_surd_gt_1[of n] by (intro notI) auto
with ‹n < l› have [simp]: "n = l - 1" by auto
show ?thesis by (simp add: as x)
qed
lemma while_option_sqrt_cfrac_info:
defines "step' ≡ (λ(as, pq). ((D' + fst pq) div snd pq # as, sqrt_remainder_step pq))"
defines "b ≡ (λ(_, pq). snd pq ≠ 1)"
defines "initial ≡ ([], (D', D - D'⇧2))"
shows "sqrt_cfrac_info D =
(case while_option b step' initial of
Some (as, _) ⇒ (Suc (length as), D', rev ((2 * D') # as)))"
proof -
have "nat (cfrac_nth (cfrac_of_real (sqrt (real D))) (Suc k)) = sqrt_cfrac_nth k" for k
by (metis nat_int sqrt_cfrac)
thus ?thesis unfolding assms while_option_sqrt_cfrac
using period_nonempty sqrt_cfrac_last
by (cases l) (auto simp: sqrt_cfrac_info_def D'_def l_def Discrete_sqrt_altdef)
qed
end
end
lemma sqrt_nat_period_length_le: "sqrt_nat_period_length D ≤ nat ⌊sqrt D⌋ * (nat ⌊sqrt D⌋ + 1)"
by (cases "is_square D") (use period_length_le_aux[of D] in auto)
lemma sqrt_nat_period_length_0_iff [simp]:
"sqrt_nat_period_length D = 0 ⟷ is_square D"
using period_nonempty[of D] by (cases "is_square D") auto
lemma sqrt_nat_period_length_pos_iff [simp]:
"sqrt_nat_period_length D > 0 ⟷ ¬is_square D"
using period_nonempty[of D] by (cases "is_square D") auto
lemma sqrt_cfrac_info_code [code]:
"sqrt_cfrac_info D =
(let D' = Discrete.sqrt D
in if D'⇧2 = D then (0, D', [])
else
case while_option
(λ(_, pq). snd pq ≠ 1)
(λ(as, (p, q)). let X = (p + D') div q; p' = X * q - p
in (X # as, p', (D - p'⇧2) div q))
([], D', D - D'⇧2)
of Some (as, _) ⇒ (Suc (length as), D', rev ((2 * D') # as)))"
proof -
define D' where "D' = Discrete.sqrt D"
show ?thesis
proof (cases "is_square D")
case True
hence "D' ^ 2 = D" by (auto simp: D'_def elim!: is_nth_powerE)
thus ?thesis using True
by (simp add: D'_def Let_def sqrt_cfrac_info_def sqrt_nat_period_length_def)
next
case False
hence "D' ^ 2 ≠ D" by (subst eq_commute) auto
thus ?thesis using while_option_sqrt_cfrac_info[OF False]
by (simp add: sqrt_cfrac_info_def D'_def Let_def
case_prod_unfold Discrete_sqrt_altdef add_ac sqrt_remainder_step_def)
qed
qed
lemma sqrt_nat_period_length_code [code]:
"sqrt_nat_period_length D = fst (sqrt_cfrac_info D)"
by (simp add: sqrt_cfrac_info_def)
text ‹
For efficiency reasons, it is often better to use an array instead of a list:
›
definition sqrt_cfrac_info_array where
"sqrt_cfrac_info_array D = (case sqrt_cfrac_info D of (a, b, c) ⇒ (a, b, IArray c))"
lemma fst_sqrt_cfrac_info_array [simp]: "fst (sqrt_cfrac_info_array D) = sqrt_nat_period_length D"
by (simp add: sqrt_cfrac_info_array_def sqrt_cfrac_info_def)
lemma snd_sqrt_cfrac_info_array [simp]: "fst (snd (sqrt_cfrac_info_array D)) = Discrete.sqrt D"
by (simp add: sqrt_cfrac_info_array_def sqrt_cfrac_info_def)
definition cfrac_sqrt_nth :: "nat × nat × nat iarray ⇒ nat ⇒ nat" where
"cfrac_sqrt_nth info n =
(case info of (l, a0, as) ⇒ if n = 0 then a0 else as !! ((n - 1) mod l))"
lemma cfrac_sqrt_nth:
assumes "¬is_square D"
shows "cfrac_nth (cfrac_of_real (sqrt D)) n =
int (cfrac_sqrt_nth (sqrt_cfrac_info_array D) n)" (is "?lhs = ?rhs")
proof (cases n)
case (Suc n')
define l where "l = sqrt_nat_period_length D"
from period_nonempty[OF assms] have "l > 0" by (simp add: l_def)
have "cfrac_nth (cfrac_of_real (sqrt D)) (Suc n') =
cfrac_nth (cfrac_of_real (sqrt D)) (Suc (n' mod l))" unfolding l_def
using cfrac_sqrt_periodic[OF assms, of n'] by simp
also have "… = map (λn. nat (cfrac_nth (cfrac_of_real (sqrt D)) (Suc n))) [0..<l] ! (n' mod l)"
using ‹l > 0› by (subst nth_map) auto
finally show ?thesis using Suc
by (simp add: sqrt_cfrac_info_array_def sqrt_cfrac_info_def l_def cfrac_sqrt_nth_def)
qed (simp_all add: sqrt_cfrac_info_def sqrt_cfrac_info_array_def
Discrete_sqrt_altdef cfrac_sqrt_nth_def)
lemma sqrt_cfrac_code [code]:
"sqrt_cfrac D =
(let info = sqrt_cfrac_info_array D;
(l, a0, _) = info
in if l = 0 then cfrac_of_int (int a0) else cfrac (cfrac_sqrt_nth info))"
proof (cases "is_square D")
case True
hence "sqrt (real D) = of_int (Discrete.sqrt D)"
by (auto elim!: is_nth_powerE)
thus ?thesis using True
by (auto simp: Let_def sqrt_cfrac_info_array_def sqrt_cfrac_info_def sqrt_cfrac_def)
next
case False
have "cfrac_sqrt_nth (sqrt_cfrac_info_array D) n > 0" if "n > 0" for n
proof -
have "int (cfrac_sqrt_nth (sqrt_cfrac_info_array D) n) > 0"
using False that by (subst cfrac_sqrt_nth [symmetric]) auto
thus ?thesis by simp
qed
moreover have "sqrt D ∉ ℚ"
using False irrat_sqrt_nonsquare by blast
ultimately have "sqrt_cfrac D = cfrac (cfrac_sqrt_nth (sqrt_cfrac_info_array D))"
using cfrac_sqrt_nth[OF False]
by (intro cfrac_eqI) (auto simp: sqrt_cfrac_def is_cfrac_def)
thus ?thesis
using False by (simp add: Let_def sqrt_cfrac_info_array_def sqrt_cfrac_info_def)
qed
text ‹
As a test, we determine the continued fraction expansion of $\sqrt{129}$, which is
$[11; \overline{2, 1, 3, 1, 6, 1, 3, 1, 2, 22}]$ (a period length of 10):
›
value "let info = sqrt_cfrac_info_array 129 in info"
value "sqrt_nat_period_length 129"
text ‹
We can also compute convergents of $\sqrt{129}$ and observe that the difference between
the square of the convergents and 129 vanishes quickly::
›
value "map (conv (sqrt_cfrac 129)) [0..<10]"
value "map (λn. ¦conv (sqrt_cfrac 129) n ^ 2 - 129¦) [0..<20]"
end