# Theory Real_Impl_Auxiliary

(* Title: Implementing field extensions of the form Q[sqrt(b)] Author: René Thiemann <rene.thiemann@uibk.ac.at> Maintainer: René Thiemann License: LGPL *) (* Copyright 2009-2014 René Thiemann This file is part of IsaFoR/CeTA. IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>. *) section ‹Auxiliary lemmas which might be moved into the Isabelle distribution.› theory Real_Impl_Auxiliary imports "HOL-Computational_Algebra.Primes" begin lemma multiplicity_prime: assumes p: "prime (i :: nat)" and ji: "j ≠ i" shows "multiplicity j i = 0" using assms by (metis dvd_refl prime_nat_iff multiplicity_eq_zero_iff multiplicity_unit_left multiplicity_zero) end

# Theory Prime_Product

(* Title: Implementing field extensions of the form Q[sqrt(b)] Author: René Thiemann <rene.thiemann@uibk.ac.at> Maintainer: René Thiemann License: LGPL *) (* Copyright 2009-2014 René Thiemann This file is part of IsaFoR/CeTA. IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>. *) section ‹Prime products› theory Prime_Product imports Real_Impl_Auxiliary Sqrt_Babylonian.Sqrt_Babylonian begin text ‹ Prime products are natural numbers where no prime factor occurs more than once. › definition prime_product where "prime_product (n :: nat) = (∀ p. prime p ⟶ multiplicity p n ≤ 1)" text ‹ The main property is that whenever $b_1$ and $b_2$ are different prime products, then $p_1 + q_1 \sqrt{b_1} = p_2 + q_2 \sqrt{b_2}$ implies $(p_1,q_1,b_1) = (p_2,q_2,b_2)$ for all rational numbers $p_1,q_1,p_2,q_2$. This is the key property to uniquely represent numbers in $\ratsb$ by triples. In the following we develop an algorithm to decompose any natural number $n$ into $n = s^2 \cdot p$ for some $s$ and prime product $p$. › (* factor a number into square * a prime product *) function prime_product_factor_main :: "nat ⇒ nat ⇒ nat ⇒ nat ⇒ nat ⇒ nat × nat" where "prime_product_factor_main factor_sq factor_pr limit n i = (if i ≤ limit ∧ i ≥ 2 then (if i dvd n then (let n' = n div i in (if i dvd n' then let n'' = n' div i in prime_product_factor_main (factor_sq * i) factor_pr (nat (root_nat_floor 3 n'')) n'' i else (case sqrt_nat n' of Cons sn _ ⇒ (factor_sq * sn, factor_pr * i) | [] ⇒ prime_product_factor_main factor_sq (factor_pr * i) (nat (root_nat_floor 3 n')) n' (Suc i) ) ) ) else prime_product_factor_main factor_sq factor_pr limit n (Suc i)) else (factor_sq, factor_pr * n))" by pat_completeness auto termination proof - let ?m1 = "λ (factor_sq :: nat,factor_pr :: nat,limit :: nat,n :: nat,i :: nat). n" let ?m2 = "λ (factor_sq,factor_pr,limit,n,i). (Suc limit - i)" { fix i have "2 ≤ i ⟹ Suc 0 < i * i" using one_less_mult[of i i] by auto } note * = this show ?thesis using wf_measures [of "[?m1, ?m2]"] by rule (auto simp add: * elim!: dvdE split: if_splits) qed lemma prime_product_factor_main: assumes "¬ (∃ s. s * s = n)" and "limit = nat (root_nat_floor 3 n)" and "m = factor_sq * factor_sq * factor_pr * n" and "prime_product_factor_main factor_sq factor_pr limit n i = (sq, p)" and "i ≥ 2" and "⋀ j. j ≥ 2 ⟹ j < i ⟹ ¬ j dvd n" and "⋀ j. prime j ⟹ j < i ⟹ multiplicity j factor_pr ≤ 1" and "⋀ j. prime j ⟹ j ≥ i ⟹ multiplicity j factor_pr = 0" and "factor_pr > 0" shows "m = sq * sq * p ∧ prime_product p" using assms proof (induct factor_sq factor_pr limit n i rule: prime_product_factor_main.induct) case (1 factor_sq factor_pr limit n i) note IH = 1(1-3) note prems = 1(4-) note simp = prems(4)[unfolded prime_product_factor_main.simps[of factor_sq factor_pr limit n i]] show ?case proof (cases "i ≤ limit") case True note i = this with prems(5) have cond: "i ≤ limit ∧ i ≥ 2" and *: "(i ≤ limit ∧ i ≥ 2) = True" by blast+ note IH = IH[OF cond] note simp = simp[unfolded * if_True] show ?thesis proof (cases "i dvd n") case False hence *: "(i dvd n) = False" by simp note simp = simp[unfolded * if_False] note IH = IH(3)[OF False prems(1-3) simp] show ?thesis proof (rule IH) fix j assume 2: "2 ≤ j" and j: "j < Suc i" from prems(6)[OF 2] j False show "¬ j dvd n" by (cases "j = i", auto) next fix j :: nat assume j: "j < Suc i" "prime j" with prems(7-8)[of j] show "multiplicity j factor_pr ≤ 1" by (cases "j = i", auto) qed (insert prems(8-9) cond, auto) next case True note mod = this hence "(i dvd n) = True" by simp note simp = simp[unfolded this if_True Let_def] note IH = IH(1,2)[OF True refl] show ?thesis proof (cases "i dvd n div i") case True hence *: "(i dvd n div i) = True" by auto define n' where "n' = n div i div i" from mod True have n: "n = n' * i * i" by (auto simp: n'_def dvd_eq_mod_eq_0) note simp = simp[unfolded * if_True split] note IH = IH(1)[OF True refl _ refl _ simp prems(5) _ prems(7-9)] show ?thesis proof (rule IH) show "m = factor_sq * i * (factor_sq * i) * factor_pr * (n div i div i)" unfolding prems(3) n'_def[symmetric] unfolding n by (auto simp: field_simps) next fix j assume "2 ≤ j" "j < i" from prems(6)[OF this] have "¬ j dvd n" by auto thus "¬ j dvd n div i div i" by (metis dvd_mult n n'_def mult.commute) next show "¬ (∃ s. s * s = n div i div i)" proof assume "∃ s. s * s = n div i div i" then obtain s where "s * s = n div i div i" by auto hence "(s * i) * (s * i) = n" unfolding n by auto with prems(1) show False by blast qed qed next case False define n' where "n' = n div i" from mod True have n: "n = n' * i" by (auto simp: n'_def dvd_eq_mod_eq_0) have prime: "prime i" unfolding prime_nat_iff proof (intro conjI allI impI) fix m assume m: "m dvd i" hence "m dvd n" unfolding n by auto with prems(6)[of m] have choice: "m ≤ 1 ∨ m ≥ i" by arith from m prems(5) have "m > 0" by (metis dvd_0_left_iff le0 le_antisym neq0_conv zero_neq_numeral) with choice have choice: "m = 1 ∨ m ≥ i" by arith from m prems(5) have "m ≤ i" by (metis False div_by_0 dvd_refl dvd_imp_le gr0I) with choice show "m = 1 ∨ m = i" by auto qed (insert prems(5), auto) from False have "(i dvd n div i) = False" by auto note simp = simp[unfolded this if_False] note IH = IH(2)[OF False _ _ refl] from prime have "i > 0" by (simp add: prime_gt_0_nat) show ?thesis proof (cases "sqrt_nat (n div i)") case (Cons s) note simp = simp[unfolded Cons list.simps] hence sq: "sq = factor_sq * s" and p: "p = factor_pr * i" by auto from arg_cong[OF Cons, of set] have s: "s * s = n div i" by auto have pp: "prime_product (factor_pr * i)" unfolding prime_product_def proof safe fix m :: nat assume m: "prime m" consider "i < m" | "i > m" | "i = m" by force thus "multiplicity m (factor_pr * i) ≤ 1" by cases (insert prems(7)[of m] prems(8)[of m] prems(9) ‹i > 0› prime m, simp_all add: multiplicity_prime prime_elem_multiplicity_mult_distrib) qed show ?thesis unfolding sq p prems(3) n unfolding n'_def s[symmetric] using pp by auto next case Nil note simp = simp[unfolded Nil list.simps] from arg_cong[OF Nil, of set] have "¬ (∃ x. x * x = n div i)" by simp note IH = IH[OF Nil this _ simp] show ?thesis proof (rule IH) show "m = factor_sq * factor_sq * (factor_pr * i) * (n div i)" unfolding prems(3) n by auto next fix j assume *: "2 ≤ j" "j < Suc i" show "¬ j dvd n div i" proof assume j: "j dvd n div i" with False have "j ≠ i" by auto with * have "2 ≤ j" "j < i" by auto from prems(6)[OF this] j show False unfolding n by (metis dvd_mult n n'_def mult.commute) qed next fix j :: nat assume "Suc i ≤ j" and j_prime: "prime j" hence ij: "i ≤ j" and j: "j ≠ i" by auto have 0: "multiplicity j i = 0" using prime j by (rule multiplicity_prime) show "multiplicity j (factor_pr * i) = 0" unfolding prems(8)[OF j_prime ij] 0 using prime j_prime j ‹ 0 < factor_pr› ‹multiplicity j factor_pr = 0› by (subst prime_elem_multiplicity_mult_distrib) (auto simp: multiplicity_prime) next fix j assume "j < Suc i" and j_prime: "prime j" hence "j < i ∨ j = i" by auto thus "multiplicity j (factor_pr * i) ≤ 1" proof assume "j = i" with prems(8)[of i] prime j_prime ‹0 < factor_pr› show ?thesis by (subst prime_elem_multiplicity_mult_distrib) auto next assume ji: "j < i" hence "j ≠ i" by auto from prems(7)[OF j_prime ji] multiplicity_prime[OF prime this] prime j_prime ‹0 < factor_pr› show ?thesis by (subst prime_elem_multiplicity_mult_distrib) auto qed qed (insert prems(5,9), auto) qed qed qed next case False hence "(i ≤ limit ∧ i ≥ 2) = False" by auto note simp = simp[unfolded this if_False] hence sq: "sq = factor_sq" and p: "p = factor_pr * n" by auto show ?thesis proof show "m = sq * sq * p" unfolding sq p prems(3) by simp show "prime_product p" unfolding prime_product_def proof safe fix m :: nat assume m: "prime m" from prems(1) have n1: "n > 1" by (cases n, auto, case_tac nat, auto) hence n0: "n > 0" by auto have "i > limit" using False by auto from this[unfolded prems(2)] have less: "int i ≥ root_nat_floor 3 n + 1" by auto have "int n < (root_nat_floor 3 n + 1) ^ 3" by (rule root_nat_floor_upper, auto) also have "… ≤ int i ^ 3" by (rule power_mono[OF less, of 3], auto) finally have n_i3: "n < i ^ 3" by (metis of_nat_less_iff of_nat_power [symmetric]) { fix m assume m: "prime m" "multiplicity m n > 0" hence mp: "m ∈ prime_factors n" by (auto simp: prime_factors_multiplicity) hence md: "m dvd n" by auto then obtain k where n: "n = m * k" .. from mp have pm: "prime m" by auto hence m2: "m ≥ 2" and m0: "m > 0" by (auto simp: prime_nat_iff) from prems(6)[OF m2] md have mi: "m ≥ i" by force { assume "multiplicity m n ≠ 1" with m have "∃ k. multiplicity m n = 2 + k" by presburger then obtain j where mult: "multiplicity m n = 2 + j" .. from n0 n have k: "k > 0" by auto from mult m0 k n m have "multiplicity m k > 0" by (auto simp: prime_elem_multiplicity_mult_distrib) with m have mp: "m ∈ prime_factors k" by (auto simp: prime_factors_multiplicity) hence md: "m dvd k" by (auto simp: k) then obtain l where kml: "k = m * l" .. note n = n[unfolded kml] from n have "l dvd n" by auto with prems(6)[of l] have "l ≤ 1 ∨ l ≥ i" by arith with n n0 have l: "l = 1 ∨ l ≥ i" by auto from n prems(1) have "l ≠ 1" by auto with l have l: "l ≥ i" by auto from mult_le_mono[OF mult_le_mono[OF mi mi] l] have "n ≥ i^3" unfolding n by (auto simp: power3_eq_cube) with n_i3 have False by auto } with mi m have "multiplicity m n = 1 ∧ m ≥ i" by auto } note n = this have "multiplicity m p = multiplicity m factor_pr + multiplicity m n" unfolding p using prems(1,9) m ‹n > 0› by (auto simp: prime_elem_multiplicity_mult_distrib) also have "… ≤ 1" proof (cases "m < i") case True from prems(7)[of m] n[of m] True m show ?thesis by force next case False hence "m ≥ i" by auto from prems(8)[OF m(1) this] n[of m] m show ?thesis by force qed finally show "multiplicity m p ≤ 1" . qed qed qed qed definition prime_product_factor :: "nat ⇒ nat × nat" where "prime_product_factor n = (case sqrt_nat n of (Cons s _) ⇒ (s,1) | [] ⇒ prime_product_factor_main 1 1 (nat (root_nat_floor 3 n)) n 2)" lemma prime_product_one[simp, intro]: "prime_product 1" unfolding prime_product_def multiplicity_one_nat by auto lemma prime_product_factor: assumes pf: "prime_product_factor n = (sq,p)" shows "n = sq * sq * p ∧ prime_product p" proof (cases "sqrt_nat n") case (Cons s) from pf[unfolded prime_product_factor_def Cons] arg_cong[OF Cons, of set] prime_product_one show ?thesis by auto next case Nil from arg_cong[OF Nil, of set] have nsq: "¬ (∃s. s * s = n)" by auto show ?thesis by (rule prime_product_factor_main[OF nsq refl, of _ 1 1 2], unfold multiplicity_one, insert pf[unfolded prime_product_factor_def Nil], auto) qed end

# Theory Real_Impl

(* Title: Implementing field extensions of the form Q[sqrt(b)] Author: René Thiemann <rene.thiemann@uibk.ac.at> Maintainer: René Thiemann License: LGPL *) (* Copyright 2009-2014 René Thiemann This file is part of IsaFoR/CeTA. IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>. *) section ‹A representation of real numbers via triples› theory Real_Impl imports Sqrt_Babylonian.Sqrt_Babylonian begin text ‹We represent real numbers of the form $p + q \cdot \sqrt{b}$ for $p,q \in \rats$, $n \in \nats$ by triples $(p,q,b)$. However, we require the invariant that $\sqrt{b}$ is irrational. Most binary operations are implemented via partial functions where the common the restriction is that the numbers $b$ in both triples have to be identical. So, we support addition of $\sqrt{2} + \sqrt{2}$, but not $\sqrt{2} + \sqrt{3}$.› text ‹The set of natural numbers whose sqrt is irrational› definition "sqrt_irrat = { q :: nat. ¬ (∃ p. p * p = rat_of_nat q)}" lemma sqrt_irrat: assumes choice: "q = 0 ∨ b ∈ sqrt_irrat" and eq: "real_of_rat p + real_of_rat q * sqrt (of_nat b) = 0" shows "q = 0" using choice proof (cases "q = 0") case False with choice have sqrt_irrat: "b ∈ sqrt_irrat" by blast from eq have "real_of_rat q * sqrt (of_nat b) = real_of_rat (- p)" by (auto simp: of_rat_minus) then obtain p where "real_of_rat q * sqrt (of_nat b) = real_of_rat p" by blast from arg_cong[OF this, of "λ x. x * x"] have "real_of_rat (q * q) * (sqrt (of_nat b) * sqrt (of_nat b)) = real_of_rat (p * p)" by (auto simp: field_simps of_rat_mult) also have "sqrt (of_nat b) * sqrt (of_nat b) = of_nat b" by simp finally have "real_of_rat (q * q * rat_of_nat b) = real_of_rat (p * p)" by (auto simp: of_rat_mult) hence "q * q * (rat_of_nat b) = p * p" by auto from arg_cong[OF this, of "λ x. x / (q * q)"] have "(p / q) * (p / q) = rat_of_nat b" using False by (auto simp: field_simps) with sqrt_irrat show ?thesis unfolding sqrt_irrat_def by blast qed text ‹To represent numbers of the form $p + q \cdot \sqrt{b}$, use mini algebraic numbers, i.e., triples $(p,q,b)$ with irrational $\sqrt{b}$.› typedef mini_alg = "{(p,q,b) | (p :: rat) (q :: rat) (b :: nat). q = 0 ∨ b ∈ sqrt_irrat}" by auto setup_lifting type_definition_mini_alg lift_definition real_of :: "mini_alg ⇒ real" is "λ (p,q,b). of_rat p + of_rat q * sqrt (of_nat b)" . lift_definition ma_of_rat :: "rat ⇒ mini_alg" is "λ x. (x,0,0)" by auto lift_definition ma_rat :: "mini_alg ⇒ rat" is fst . lift_definition ma_base :: "mini_alg ⇒ nat" is "snd o snd" . lift_definition ma_coeff :: "mini_alg ⇒ rat" is "fst o snd" . lift_definition ma_uminus :: "mini_alg ⇒ mini_alg" is "λ (p1,q1,b1). (- p1, - q1, b1)" by auto lift_definition ma_compatible :: "mini_alg ⇒ mini_alg ⇒ bool" is "λ (p1,q1,b1) (p2,q2,b2). q1 = 0 ∨ q2 = 0 ∨ b1 = b2" . definition ma_normalize :: "rat × rat × nat ⇒ rat × rat × nat" where "ma_normalize x ≡ case x of (a,b,c) ⇒ if b = 0 then (a,0,0) else (a,b,c)" lemma ma_normalize_case[simp]: "(case ma_normalize r of (a,b,c) ⇒ real_of_rat a + real_of_rat b * sqrt (of_nat c)) = (case r of (a,b,c) ⇒ real_of_rat a + real_of_rat b * sqrt (of_nat c))" by (cases r, auto simp: ma_normalize_def) lift_definition ma_plus :: "mini_alg ⇒ mini_alg ⇒ mini_alg" is "λ (p1,q1,b1) (p2,q2,b2). if q1 = 0 then (p1 + p2, q2, b2) else ma_normalize (p1 + p2, q1 + q2, b1)" by (auto simp: ma_normalize_def) lift_definition ma_times :: "mini_alg ⇒ mini_alg ⇒ mini_alg" is "λ (p1,q1,b1) (p2,q2,b2). if q1 = 0 then ma_normalize (p1*p2, p1*q2, b2) else ma_normalize (p1*p2 + of_nat b2*q1*q2, p1*q2 + q1*p2, b1)" by (auto simp: ma_normalize_def) lift_definition ma_inverse :: "mini_alg ⇒ mini_alg" is "λ (p,q,b). let d = inverse (p * p - of_nat b * q * q) in ma_normalize (p * d, - q * d, b)" by (auto simp: Let_def ma_normalize_def) lift_definition ma_floor :: "mini_alg ⇒ int" is "λ (p,q,b). case (quotient_of p,quotient_of q) of ((z1,n1),(z2,n2)) ⇒ let z2n1 = z2 * n1; z1n2 = z1 * n2; n12 = n1 * n2; prod = z2n1 * z2n1 * int b in (z1n2 + (if z2n1 ≥ 0 then sqrt_int_floor_pos prod else - sqrt_int_ceiling_pos prod)) div n12" . lift_definition ma_sqrt :: "mini_alg ⇒ mini_alg" is "λ (p,q,b). let (a,b) = quotient_of p; aa = abs (a * b) in case sqrt_int aa of [] ⇒ (0,inverse (of_int b),nat aa) | (Cons s _) ⇒ (of_int s / of_int b,0,0)" proof (unfold Let_def) fix prod :: "rat × rat × nat" obtain p q b where prod: "prod = (p,q,b)" by (cases prod, auto) obtain a b where p: "quotient_of p = (a,b)" by force show "∃p q b. (case prod of (p, q, b) ⇒ case quotient_of p of (a, b) ⇒ (case sqrt_int ¦a * b¦ of [] ⇒ (0, inverse (of_int b), nat ¦a * b¦) | s # x ⇒ (of_int s / of_int b, 0, 0))) = (p, q, b) ∧ (q = 0 ∨ b ∈ sqrt_irrat)" proof (cases "sqrt_int (abs (a * b))") case Nil from sqrt_int[of "abs (a * b)"] Nil have irrat: "¬ (∃ y. y * y = ¦a * b¦)" by auto have "nat ¦a * b¦ ∈ sqrt_irrat" proof (rule ccontr) assume "nat ¦a * b¦ ∉ sqrt_irrat" then obtain x :: rat where "x * x = rat_of_nat (nat ¦a * b¦)" unfolding sqrt_irrat_def by auto hence "x * x = rat_of_int ¦a * b¦" by auto from sqr_rat_of_int[OF this] irrat show False by blast qed thus ?thesis using Nil by (auto simp: prod p) qed (auto simp: prod p Cons) qed lift_definition ma_equal :: "mini_alg ⇒ mini_alg ⇒ bool" is "λ (p1,q1,b1) (p2,q2,b2). p1 = p2 ∧ q1 = q2 ∧ (q1 = 0 ∨ b1 = b2)" . lift_definition ma_ge_0 :: "mini_alg ⇒ bool" is "λ (p,q,b). let bqq = of_nat b * q * q; pp = p * p in 0 ≤ p ∧ bqq ≤ pp ∨ 0 ≤ q ∧ pp ≤ bqq" . lift_definition ma_is_rat :: "mini_alg ⇒ bool" is "λ (p,q,b). q = 0" . definition ge_0 :: "real ⇒ bool" where [code del]: "ge_0 x = (x ≥ 0)" lemma ma_ge_0: "ge_0 (real_of x) = ma_ge_0 x" proof (transfer, unfold Let_def, clarsimp) fix p' q' :: rat and b' :: nat assume b': "q' = 0 ∨ b' ∈ sqrt_irrat" define b where "b = real_of_nat b'" define p where "p = real_of_rat p'" define q where "q = real_of_rat q'" from b' have b: "0 ≤ b" "q = 0 ∨ b' ∈ sqrt_irrat" unfolding b_def q_def by auto define qb where "qb = q * sqrt b" from b have sqrt: "sqrt b ≥ 0" by auto from b(2) have disj: "q = 0 ∨ b ≠ 0" unfolding sqrt_irrat_def b_def by auto have bdef: "b = real_of_rat (of_nat b')" unfolding b_def by auto have "(0 ≤ p + q * sqrt b) = (0 ≤ p + qb)" unfolding qb_def by simp also have "… ⟷ (0 ≤ p ∧ abs qb ≤ abs p ∨ 0 ≤ qb ∧ abs p ≤ abs qb)" by arith also have "… ⟷ (0 ≤ p ∧ qb * qb ≤ p * p ∨ 0 ≤ qb ∧ p * p ≤ qb * qb)" unfolding abs_lesseq_square .. also have "qb * qb = b * q * q" unfolding qb_def using b by auto also have "0 ≤ qb ⟷ 0 ≤ q" unfolding qb_def using sqrt disj by (metis le_cases mult_eq_0_iff mult_nonneg_nonneg mult_nonpos_nonneg order_class.order.antisym qb_def real_sqrt_eq_zero_cancel_iff) also have "(0 ≤ p ∧ b * q * q ≤ p * p ∨ 0 ≤ q ∧ p * p ≤ b * q * q) ⟷ (0 ≤ p' ∧ of_nat b' * q' * q' ≤ p' * p' ∨ 0 ≤ q' ∧ p' * p' ≤ of_nat b' * q' * q')" unfolding qb_def by (unfold bdef p_def q_def of_rat_mult[symmetric] of_rat_less_eq, simp) finally show "ge_0 (real_of_rat p' + real_of_rat q' * sqrt (of_nat b')) = (0 ≤ p' ∧ of_nat b' * q' * q' ≤ p' * p' ∨ 0 ≤ q' ∧ p' * p' ≤ of_nat b' * q' * q')" unfolding ge_0_def p_def b_def q_def by (auto simp: of_rat_add of_rat_mult) qed lemma ma_0: "0 = real_of (ma_of_rat 0)" by (transfer, auto) lemma ma_1: "1 = real_of (ma_of_rat 1)" by (transfer, auto) lemma ma_uminus: "- (real_of x) = real_of (ma_uminus x)" by (transfer, auto simp: of_rat_minus) lemma ma_inverse: "inverse (real_of r) = real_of (ma_inverse r)" proof (transfer, unfold Let_def, clarsimp) fix p' q' :: rat and b' :: nat assume b': "q' = 0 ∨ b' ∈ sqrt_irrat" define b where "b = real_of_nat b'" define p where "p = real_of_rat p'" define q where "q = real_of_rat q'" from b' have b: "b ≥ 0" "q = 0 ∨ b' ∈ sqrt_irrat" unfolding b_def q_def by auto have "inverse (p + q * sqrt b) = (p - q * sqrt b) * inverse (p * p - b * (q * q))" proof (cases "q = 0") case True thus ?thesis by (cases "p = 0", auto simp: field_simps) next case False from sqrt_irrat[OF b', of p'] b_def p_def q_def False have nnull: "p + q * sqrt b ≠ 0" by auto have "?thesis ⟷ (p + q * sqrt b) * inverse (p + q * sqrt b) = (p + q * sqrt b) * ((p - q * sqrt b) * inverse (p * p - b * (q * q)))" unfolding mult_left_cancel[OF nnull] by auto also have "(p + q * sqrt b) * inverse (p + q * sqrt b) = 1" using nnull by auto also have "(p + q * sqrt b) * ((p - q * sqrt b) * inverse (p * p - b * (q * q))) = (p * p - b * (q * q)) * inverse (p * p - b * (q * q))" using b by (auto simp: field_simps) also have "... = 1" proof (rule right_inverse, rule) assume eq: "p * p - b * (q * q) = 0" have "real_of_rat (p' * p' / (q' * q')) = p * p / (q * q)" unfolding p_def b_def q_def by (auto simp: of_rat_mult of_rat_divide) also have "… = b" using eq False by (auto simp: field_simps) also have "… = real_of_rat (of_nat b')" unfolding b_def by auto finally have "(p' / q') * (p' / q') = of_nat b'" unfolding of_rat_eq_iff by simp with b False show False unfolding sqrt_irrat_def by blast qed finally show ?thesis by simp qed thus "inverse (real_of_rat p' + real_of_rat q' * sqrt (of_nat b')) = real_of_rat (p' * inverse (p' * p' - of_nat b' * q' * q')) + real_of_rat (- (q' * inverse (p' * p' - of_nat b' * q' * q'))) * sqrt (of_nat b')" by (simp add: divide_simps of_rat_mult of_rat_divide of_rat_diff of_rat_minus b_def p_def q_def split: if_splits) qed lemma ma_sqrt_main: "ma_rat r ≥ 0 ⟹ ma_coeff r = 0 ⟹ sqrt (real_of r) = real_of (ma_sqrt r)" proof (transfer, unfold Let_def, clarsimp) fix p :: rat assume p: "0 ≤ p" hence abs: "abs p = p" by auto obtain a b where ab: "quotient_of p = (a,b)" by force hence pab: "p = of_int a / of_int b" by (rule quotient_of_div) from ab have b: "b > 0" by (rule quotient_of_denom_pos) with p pab have abpos: "a * b ≥ 0" by (metis of_int_0_le_iff of_int_le_0_iff zero_le_divide_iff zero_le_mult_iff) have rab: "of_nat (nat (a * b)) = real_of_int a * real_of_int b" using abpos by simp let ?lhs = "sqrt (of_int a / of_int b)" let ?rhs = "(case case quotient_of p of (a, b) ⇒ (case sqrt_int ¦a * b¦ of [] ⇒ (0, inverse (of_int b), nat ¦a * b¦) | s # x ⇒ (of_int s / of_int b, 0, 0)) of (p, q, b) ⇒ of_rat p + of_rat q * sqrt (of_nat b))" have "sqrt (real_of_rat p) = ?lhs" unfolding pab by (metis of_rat_divide of_rat_of_int_eq) also have "… = ?rhs" proof (cases "sqrt_int ¦a * b¦") case Nil define sb where "sb = sqrt (of_int b)" define sa where "sa = sqrt (of_int a)" from b sb_def have sb: "sb > 0" "real_of_int b > 0" by auto have sbb: "sb * sb = real_of_int b" unfolding sb_def by (rule sqrt_sqrt, insert b, auto) from Nil have "?thesis = (sa / sb = inverse (of_int b) * (sa * sb))" unfolding ab sa_def sb_def using abpos by (simp add: rab of_rat_divide real_sqrt_mult real_sqrt_divide of_rat_inverse) also have "… = (sa = inverse (of_int b) * sa * (sb * sb))" using sb by (metis b divide_real_def eq_divide_imp inverse_divide inverse_inverse_eq inverse_mult_distrib less_int_code(1) of_int_eq_0_iff real_sqrt_eq_zero_cancel_iff sb_def sbb times_divide_eq_right) also have "… = True" using sb(2) unfolding sbb by auto finally show "?thesis" by simp next case (Cons s x) from b have b: "real_of_int b > 0" by auto from Cons sqrt_int[of "abs (a * b)"] have "s * s = abs (a * b)" by auto with sqrt_int_pos[OF Cons] have "sqrt (real_of_int (abs (a * b))) = of_int s" by (metis abs_of_nonneg of_int_mult of_int_abs real_sqrt_abs2) with abpos have "of_int s = sqrt (real_of_int (a * b))" by auto thus ?thesis unfolding ab split using Cons b by (auto simp: of_rat_divide field_simps real_sqrt_divide real_sqrt_mult) qed finally show "sqrt (real_of_rat p) = ?rhs" . qed lemma ma_sqrt: "sqrt (real_of r) = (if ma_coeff r = 0 then (if ma_rat r ≥ 0 then real_of (ma_sqrt r) else - real_of (ma_sqrt (ma_uminus r))) else Code.abort (STR ''cannot represent sqrt of irrational number'') (λ _. sqrt (real_of r)))" proof (cases "ma_coeff r = 0") case True note 0 = this hence 00: "ma_coeff (ma_uminus r) = 0" by (transfer, auto) show ?thesis proof (cases "ma_rat r ≥ 0") case True from ma_sqrt_main[OF this 0] 0 True show ?thesis by auto next case False hence "ma_rat (ma_uminus r) ≥ 0" by (transfer, auto) from ma_sqrt_main[OF this 00, folded ma_uminus, symmetric] False 0 show ?thesis by (auto simp: real_sqrt_minus) qed qed auto lemma ma_plus: "(real_of r1 + real_of r2) = (if ma_compatible r1 r2 then real_of (ma_plus r1 r2) else Code.abort (STR ''different base'') (λ _. real_of r1 + real_of r2))" by transfer (auto split: prod.split simp: field_simps of_rat_add) lemma ma_times: "(real_of r1 * real_of r2) = (if ma_compatible r1 r2 then real_of (ma_times r1 r2) else Code.abort (STR ''different base'') (λ _. real_of r1 * real_of r2))" by transfer (auto split: prod.split simp: field_simps of_rat_mult of_rat_add) lemma ma_equal: "HOL.equal (real_of r1) (real_of r2) = (if ma_compatible r1 r2 then ma_equal r1 r2 else Code.abort (STR ''different base'') (λ _. HOL.equal (real_of r1) (real_of r2)))" proof (transfer, unfold equal_real_def, clarsimp) fix p1 q1 p2 q2 :: rat and b1 b2 :: nat assume b1: "q1 = 0 ∨ b1 ∈ sqrt_irrat" assume b2: "q2 = 0 ∨ b2 ∈ sqrt_irrat" assume base: "q1 = 0 ∨ q2 = 0 ∨ b1 = b2" let ?l = "real_of_rat p1 + real_of_rat q1 * sqrt (of_nat b1) = real_of_rat p2 + real_of_rat q2 * sqrt (of_nat b2)" let ?m = "real_of_rat q1 * sqrt (of_nat b1) = real_of_rat (p2 - p1) + real_of_rat q2 * sqrt (of_nat b2)" let ?r = "p1 = p2 ∧ q1 = q2 ∧ (q1 = 0 ∨ b1 = b2)" have "?l ⟷ real_of_rat q1 * sqrt (of_nat b1) = real_of_rat (p2 - p1) + real_of_rat q2 * sqrt (of_nat b2)" by (auto simp: of_rat_add of_rat_diff of_rat_minus) also have "… ⟷ p1 = p2 ∧ q1 = q2 ∧ (q1 = 0 ∨ b1 = b2)" proof assume ?m from base have "q1 = 0 ∨ q1 ≠ 0 ∧ q2 = 0 ∨ q1 ≠ 0 ∧ q2 ≠ 0 ∧ b1 = b2" by auto thus "p1 = p2 ∧ q1 = q2 ∧ (q1 = 0 ∨ b1 = b2)" proof assume q1: "q1 = 0" with ‹?m› have "real_of_rat (p2 - p1) + real_of_rat q2 * sqrt (of_nat b2) = 0" by auto with sqrt_irrat b2 have q2: "q2 = 0" by auto with q1 ‹?m› show ?thesis by auto next assume "q1 ≠ 0 ∧ q2 = 0 ∨ q1 ≠ 0 ∧ q2 ≠ 0 ∧ b1 = b2" thus ?thesis proof assume ass: "q1 ≠ 0 ∧ q2 = 0" with ‹?m› have "real_of_rat (p1 - p2) + real_of_rat q1 * sqrt (of_nat b1) = 0" by (auto simp: of_rat_diff) with b1 have "q1 = 0" using sqrt_irrat by auto with ass show ?thesis by auto next assume ass: "q1 ≠ 0 ∧ q2 ≠ 0 ∧ b1 = b2" with ‹?m› have *: "real_of_rat (p2 - p1) + real_of_rat (q2 - q1) * sqrt (of_nat b2) = 0" by (auto simp: field_simps of_rat_diff) have "q2 - q1 = 0" by (rule sqrt_irrat[OF _ *], insert ass b2, auto) with * ass show ?thesis by auto qed qed qed auto finally show "?l = ?r" by simp qed lemma ma_floor: "floor (real_of r) = ma_floor r" proof (transfer, unfold Let_def, clarsimp) fix p q :: rat and b :: nat obtain z1 n1 where qp: "quotient_of p = (z1,n1)" by force obtain z2 n2 where qq: "quotient_of q = (z2,n2)" by force from quotient_of_denom_pos[OF qp] have n1: "0 < n1" . from quotient_of_denom_pos[OF qq] have n2: "0 < n2" . from quotient_of_div[OF qp] have p: "p = of_int z1 / of_int n1" . from quotient_of_div[OF qq] have q: "q = of_int z2 / of_int n2" . have p: "p = of_int (z1 * n2) / of_int (n1 * n2)" unfolding p using n2 by auto have q: "q = of_int (z2 * n1) / of_int (n1 * n2)" unfolding q using n1 by auto define z1n2 where "z1n2 = z1 * n2" define z2n1 where "z2n1 = z2 * n1" define n12 where "n12 = n1 * n2" define r_add where "r_add = of_int (z2n1) * sqrt (real_of_int (int b))" from n1 n2 have n120: "n12 > 0" unfolding n12_def by simp have "floor (of_rat p + of_rat q * sqrt (real_of_nat b)) = floor ((of_int z1n2 + r_add) / of_int n12)" unfolding r_add_def n12_def z1n2_def z2n1_def unfolding p q add_divide_distrib of_rat_divide of_rat_of_int_eq of_int_of_nat_eq by simp also have "… = floor (of_int z1n2 + r_add) div n12" by (rule floor_div_pos_int[OF n120]) also have "of_int z1n2 + r_add = r_add + of_int z1n2" by simp also have "floor (…) = floor r_add + z1n2" by simp also have "… = z1n2 + floor r_add" by simp finally have id: "⌊of_rat p + of_rat q * sqrt (of_nat b)⌋ = (z1n2 + ⌊r_add⌋) div n12" . show "⌊of_rat p + of_rat q * sqrt (of_nat b)⌋ = (case quotient_of p of (z1, n1) ⇒ case quotient_of q of (z2, n2) ⇒ (z1 * n2 + (if 0 ≤ z2 * n1 then sqrt_int_floor_pos (z2 * n1 * (z2 * n1) * int b) else - sqrt_int_ceiling_pos (z2 * n1 * (z2 * n1) * int b))) div (n1 * n2))" unfolding qp qq split id n12_def z1n2_def proof (rule arg_cong[of _ _ "λ x. ((z1 * n2) + x) div (n1 * n2)"]) have ge_int: "z2 * n1 * (z2 * n1) * int b ≥ 0" by (metis mult_nonneg_nonneg zero_le_square of_nat_0_le_iff) show "⌊r_add⌋ = (if 0 ≤ z2 * n1 then sqrt_int_floor_pos (z2 * n1 * (z2 * n1) * int b) else - sqrt_int_ceiling_pos (z2 * n1 * (z2 * n1) * int b))" proof (cases "z2 * n1 ≥ 0") case True hence ge: "real_of_int (z2 * n1) ≥ 0" by (metis of_int_0_le_iff) have radd: "r_add = sqrt (of_int (z2 * n1 * (z2 * n1) * int b))" unfolding r_add_def z2n1_def using sqrt_sqrt[OF ge] by (simp add: ac_simps real_sqrt_mult) show ?thesis unfolding radd sqrt_int_floor_pos[OF ge_int] using True by simp next case False hence ge: "real_of_int (- (z2 * n1)) ≥ 0" by (metis mult_zero_left neg_0_le_iff_le of_int_0_le_iff order_refl zero_le_mult_iff) have "r_add = - sqrt (of_int (z2 * n1 * (z2 * n1) * int b))" unfolding r_add_def z2n1_def using sqrt_sqrt[OF ge] by (metis minus_minus minus_mult_commute minus_mult_right of_int_minus of_int_mult real_sqrt_minus real_sqrt_mult z2n1_def) hence radd: "floor r_add = - ceiling (sqrt (of_int (z2 * n1 * (z2 * n1) * int b)))" by (metis ceiling_def minus_minus) show ?thesis unfolding radd sqrt_int_ceiling_pos[OF ge_int] using False by simp qed qed qed lemma comparison_impl: "(x :: real) ≤ (y :: real) = ge_0 (y - x)" "(x :: real) < (y :: real) = (x ≠ y ∧ ge_0 (y - x))" by (simp_all add: ge_0_def, linarith+) lemma ma_of_rat: "real_of_rat r = real_of (ma_of_rat r)" by (transfer, auto) definition is_rat :: "real ⇒ bool" where [code_abbrev]: "is_rat x ⟷ x ∈ ℚ" lemma ma_is_rat: "is_rat (real_of x) = ma_is_rat x" proof (transfer, unfold is_rat_def, clarsimp) fix p q :: rat and b :: nat let ?p = "real_of_rat p" let ?q = "real_of_rat q" let ?b = "real_of_nat b" let ?b' = "real_of_rat (of_nat b)" assume b: "q = 0 ∨ b ∈ sqrt_irrat" show "(?p + ?q * sqrt ?b ∈ ℚ) = (q = 0)" proof (cases "q = 0") case False from False b have b: "b ∈ sqrt_irrat" by auto { assume "?p + ?q * sqrt ?b ∈ ℚ" from this[unfolded Rats_def] obtain r where r: "?p + ?q * sqrt ?b = real_of_rat r" by auto let ?r = "real_of_rat r" from r have "real_of_rat (p - r) + ?q * sqrt ?b = 0" by (simp add: of_rat_diff) from sqrt_irrat[OF disjI2[OF b] this] False have False by auto } thus ?thesis by auto qed auto qed (* compute all numbers y for which y * y = x, if x ∈ ℚ, otherwise return empty list. of course, one could have returned [-sqrt x, sqrt x], but this might result in runtime errors, e.g., if sqrt_real (sqrt 2) would be invoked. The current formulation is guaranteed to not crash. *) definition "sqrt_real x = (if x ∈ ℚ ∧ x ≥ 0 then (if x = 0 then [0] else (let sx = sqrt x in [sx,-sx])) else [])" lemma sqrt_real[simp]: assumes x: "x ∈ ℚ" shows "set (sqrt_real x) = {y . y * y = x}" proof (cases "x ≥ 0") case False hence "⋀ y. y * y ≠ x" by auto with False show ?thesis unfolding sqrt_real_def by auto next case True thus ?thesis using x by (cases "x = 0", auto simp: Let_def sqrt_real_def) qed code_datatype real_of declare [[code drop: "plus :: real ⇒ real ⇒ real" "uminus :: real ⇒ real" "times :: real ⇒ real ⇒ real" "inverse :: real ⇒ real" "floor :: real ⇒ int" sqrt "HOL.equal :: real ⇒ real ⇒ bool" ]] lemma [code]: "Ratreal = real_of ∘ ma_of_rat" by (simp add: fun_eq_iff) (transfer, simp) lemmas ma_code_eqns [code equation] = ma_ge_0 ma_floor ma_0 ma_1 ma_uminus ma_inverse ma_sqrt ma_plus ma_times ma_equal ma_is_rat comparison_impl lemma [code equation]: "(x :: real) / (y :: real) = x * inverse y" "(x :: real) - (y :: real) = x + (- y)" by (simp_all add: divide_inverse) text ‹Some tests with small numbers. To work on larger number, one should additionally import the theories for efficient calculation on numbers› value "⌊101.1 * (3 * sqrt 2 + 6 * sqrt 0.5)⌋" value "⌊606.2 * sqrt 2 + 0.001⌋" value "101.1 * (3 * sqrt 2 + 6 * sqrt 0.5) = 606.2 * sqrt 2 + 0.001" value "101.1 * (3 * sqrt 2 + 6 * sqrt 0.5) > 606.2 * sqrt 2 + 0.001" value "(sqrt 0.1 ∈ ℚ, sqrt (- 0.09) ∈ ℚ)" end

# Theory Real_Unique_Impl

(* Title: Implementing field extensions of the form Q[sqrt(b)] Author: René Thiemann <rene.thiemann@uibk.ac.at> Maintainer: René Thiemann License: LGPL *) (* Copyright 2009-2014 René Thiemann This file is part of IsaFoR/CeTA. IsaFoR/CeTA is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. IsaFoR/CeTA is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with IsaFoR/CeTA. If not, see <http://www.gnu.org/licenses/>. *) section ‹A unique representation of real numbers via triples› theory Real_Unique_Impl imports Prime_Product Real_Impl Show.Show_Instances Show.Show_Real begin text ‹We implement the real numbers again using triples, but now we require an additional invariant on the triples, namely that the base has to be a prime product. This has the consequence that the mapping of triples into $\reals$ is injective. Hence, equality on reals is now equality on triples, which can even be executed in case of different bases. Similarly, we now also allow different basis in comparisons. Ultimately, injectivity allows us to define a show-function for real numbers, which pretty prints real numbers into strings.› typedef mini_alg_unique = "{ r :: mini_alg . ma_coeff r = 0 ∧ ma_base r = 0 ∨ ma_coeff r ≠ 0 ∧ prime_product (ma_base r)}" by (transfer, auto) setup_lifting type_definition_mini_alg_unique lift_definition real_of_u :: "mini_alg_unique ⇒ real" is real_of . lift_definition mau_floor :: "mini_alg_unique ⇒ int" is ma_floor . lift_definition mau_of_rat :: "rat ⇒ mini_alg_unique" is ma_of_rat by (transfer, auto) lift_definition mau_rat :: "mini_alg_unique ⇒ rat" is ma_rat . lift_definition mau_base :: "mini_alg_unique ⇒ nat" is ma_base . lift_definition mau_coeff :: "mini_alg_unique ⇒ rat" is ma_coeff . lift_definition mau_uminus :: "mini_alg_unique ⇒ mini_alg_unique" is ma_uminus by (transfer, auto) lift_definition mau_compatible :: "mini_alg_unique ⇒ mini_alg_unique ⇒ bool" is ma_compatible . lift_definition mau_ge_0 :: "mini_alg_unique ⇒ bool" is ma_ge_0 . lift_definition mau_inverse :: "mini_alg_unique ⇒ mini_alg_unique" is ma_inverse by (transfer, auto simp: ma_normalize_def Let_def split: if_splits) lift_definition mau_plus :: "mini_alg_unique ⇒ mini_alg_unique ⇒ mini_alg_unique" is ma_plus by (transfer, auto simp: ma_normalize_def split: if_splits) lift_definition mau_times :: "mini_alg_unique ⇒ mini_alg_unique ⇒ mini_alg_unique" is ma_times by (transfer, auto simp: ma_normalize_def split: if_splits) lift_definition ma_identity :: "mini_alg ⇒ mini_alg ⇒ bool" is "(=)" . lift_definition mau_equal :: "mini_alg_unique ⇒ mini_alg_unique ⇒ bool" is ma_identity . lift_definition mau_is_rat :: "mini_alg_unique ⇒ bool" is ma_is_rat . lemma Ratreal_code[code]: "Ratreal = real_of_u ∘ mau_of_rat" by (simp add: fun_eq_iff) (transfer, transfer, simp) lemma mau_floor: "floor (real_of_u r) = mau_floor r" using ma_floor by (transfer, auto) lemma mau_inverse: "inverse (real_of_u r) = real_of_u (mau_inverse r)" using ma_inverse by (transfer, auto) lemma mau_uminus: "- (real_of_u r) = real_of_u (mau_uminus r)" using ma_uminus by (transfer, auto) lemma mau_times: "(real_of_u r1 * real_of_u r2) = (if mau_compatible r1 r2 then real_of_u (mau_times r1 r2) else Code.abort (STR ''different base'') (λ _. real_of_u r1 * real_of_u r2))" using ma_times by (transfer, auto) lemma mau_plus: "(real_of_u r1 + real_of_u r2) = (if mau_compatible r1 r2 then real_of_u (mau_plus r1 r2) else Code.abort (STR ''different base'') (λ _. real_of_u r1 + real_of_u r2))" using ma_plus by (transfer, auto) lemma real_of_u_inj[simp]: "real_of_u x = real_of_u y ⟷ x = y" proof note field_simps[simp] of_rat_diff[simp] assume "real_of_u x = real_of_u y" thus "x = y" proof (transfer) fix x y assume "ma_coeff x = 0 ∧ ma_base x = 0 ∨ ma_coeff x ≠ 0 ∧ prime_product (ma_base x)" and "ma_coeff y = 0 ∧ ma_base y = 0 ∨ ma_coeff y ≠ 0 ∧ prime_product (ma_base y)" and "real_of x = real_of y" thus "x = y" proof (transfer, clarsimp) fix p1 q1 p2 q2 :: rat and b1 b2 let ?p1 = "real_of_rat p1" let ?p2 = "real_of_rat p2" let ?q1 = "real_of_rat q1" let ?q2 = "real_of_rat q2" let ?b1 = "real_of_nat b1" let ?b2 = "real_of_nat b2" assume q1: "q1 = 0 ∧ b1 = 0 ∨ q1 ≠ 0 ∧ prime_product b1" and q2: "q2 = 0 ∧ b2 = 0 ∨ q2 ≠ 0 ∧ prime_product b2" and i1: "q1 = 0 ∨ b1 ∈ sqrt_irrat" and i2: "q2 = 0 ∨ b2 ∈ sqrt_irrat" and eq: "?p1 + ?q1 * sqrt ?b1 = ?p2 + ?q2 * sqrt ?b2" show "p1 = p2 ∧ q1 = q2 ∧ b1 = b2" proof (cases "q1 = 0") case True have "q2 = 0" by (rule sqrt_irrat[OF i2, of "p2 - p1"], insert eq True q1, auto) with True q1 q2 eq show ?thesis by auto next case False hence 1: "q1 ≠ 0" "prime_product b1" using q1 by auto { assume *: "q2 = 0" have "q1 = 0" by (rule sqrt_irrat[OF i1, of "p1 - p2"], insert eq * q2, auto) with False have False by auto } hence 2: "q2 ≠ 0" "prime_product b2" using q2 by auto from 1 i1 have b1: "b1 ≠ 0" unfolding sqrt_irrat_def by (cases b1, auto) from 2 i2 have b2: "b2 ≠ 0" unfolding sqrt_irrat_def by (cases b2, auto) let ?sq = "λ x. x * x" define q3 where "q3 = p2 - p1" let ?q3 = "real_of_rat q3" let ?e = "of_rat (q2 * q2 * of_nat b2 + ?sq q3 - ?sq q1 * of_nat b1) + of_rat (2 * q2 * q3) * sqrt ?b2" from eq have *: "?q1 * sqrt ?b1 = ?q2 * sqrt ?b2 + ?q3" by (simp add: q3_def) from arg_cong[OF this, of ?sq] have "0 = (real_of_rat 2 * ?q2 * ?q3) * sqrt ?b2 + (?sq ?q2 * ?b2 + ?sq ?q3 - ?sq ?q1 * ?b1)" by auto also have "… = ?e" by (simp add: of_rat_mult of_rat_add of_rat_minus) finally have eq: "?e = 0" by simp from sqrt_irrat[OF _ this] 2 i2 have q3: "q3 = 0" by auto hence p: "p1 = p2" unfolding q3_def by simp let ?b1 = "rat_of_nat b1" let ?b2 = "rat_of_nat b2" from eq[unfolded q3] have eq: "?sq q2 * ?b2 = ?sq q1 * ?b1" by auto obtain z1 n1 where d1: "quotient_of q1 = (z1,n1)" by force obtain z2 n2 where d2: "quotient_of q2 = (z2,n2)" by force note id = quotient_of_div[OF d1] quotient_of_div[OF d2] note pos = quotient_of_denom_pos[OF d1] quotient_of_denom_pos[OF d2] from id(1) 1(1) pos(1) have z1: "z1 ≠ 0" by auto from id(2) 2(1) pos(2) have z2: "z2 ≠ 0" by auto let ?n1 = "rat_of_int n1" let ?n2 = "rat_of_int n2" let ?z1 = "rat_of_int z1" let ?z2 = "rat_of_int z2" from arg_cong[OF eq[simplified id], of "λ x. x * ?sq ?n1 * ?sq ?n2", simplified field_simps] have "?sq (?n1 * ?z2) * ?b2 = ?sq (?n2 * ?z1) * ?b1" using pos by auto moreover have "?n1 * ?z2 ≠ 0" "?n2 * ?z1 ≠ 0" using z1 z2 pos by auto ultimately obtain i1 i2 where 0: "rat_of_int i1 ≠ 0" "rat_of_int i2 ≠ 0" and eq: "?sq (rat_of_int i2) * ?b2 = ?sq (rat_of_int i1) * ?b1" unfolding of_int_mult[symmetric] by blast+ let ?b1 = "int b1" let ?b2 = "int b2" from eq have eq: "?sq i1 * ?b1 = ?sq i2 * ?b2" by (metis (hide_lams, no_types) of_int_eq_iff of_int_mult of_int_of_nat_eq) from 0 have 0: "i1 ≠ 0" "i2 ≠ 0" by auto from arg_cong[OF eq, of nat] have "?sq (nat (abs i1)) * b1 = ?sq (nat (abs i2)) * b2" by (metis abs_of_nat eq nat_abs_mult_distrib nat_int) moreover have "nat (abs i1) > 0" "nat (abs i2) > 0" using 0 by auto ultimately obtain n1 n2 where 0: "n1 > 0" "n2 > 0" and eq: "?sq n1 * b1 = ?sq n2 * b2" by blast from b1 0 have b1: "b1 > 0" "n1 > 0" "n1 * n1 > 0" by auto from b2 0 have b2: "b2 > 0" "n2 > 0" "n2 * n2 > 0" by auto { fix p :: nat assume p: "prime p" have "multiplicity p (?sq n1 * b1) = multiplicity p b1 + 2 * multiplicity p n1" using b1 p by (auto simp: prime_elem_multiplicity_mult_distrib) also have "… mod 2 = multiplicity p b1 mod 2" by presburger finally have id1: "multiplicity p (?sq n1 * b1) mod 2 = multiplicity p b1 mod 2" . have "multiplicity p (?sq n2 * b2) = multiplicity p b2 + 2 * multiplicity p n2" using b2 p by (auto simp: prime_elem_multiplicity_mult_distrib) also have "… mod 2 = multiplicity p b2 mod 2" by presburger finally have id2: "multiplicity p (?sq n2 * b2) mod 2 = multiplicity p b2 mod 2" . from id1 id2 eq have eq: "multiplicity p b1 mod 2 = multiplicity p b2 mod 2" by simp from 1(2) 2(2) p have "multiplicity p b1 ≤ 1" "multiplicity p b2 ≤ 1" unfolding prime_product_def by auto with eq have "multiplicity p b1 = multiplicity p b2" by simp } with b1(1) b2(1) have b: "b1 = b2" by (rule multiplicity_eq_nat) from *[unfolded b q3] b1(1) b2(1) have q: "q1 = q2" by simp from p q b show ?thesis by blast qed qed qed qed simp lift_definition mau_sqrt :: "mini_alg_unique ⇒ mini_alg_unique" is "λ ma. let (a,b) = quotient_of (ma_rat ma); (sq,fact) = prime_product_factor (nat (abs a * b)); ma' = ma_of_rat (of_int (sgn(a)) * of_nat sq / of_int b) in ma_times ma' (ma_sqrt (ma_of_rat (of_nat fact)))" proof - fix ma :: mini_alg let ?num = "let (a, b) = quotient_of (ma_rat ma); (sq, fact) = prime_product_factor (nat (¦a¦ * b)); ma' = ma_of_rat (rat_of_int (sgn a) * rat_of_nat sq / of_int b) in ma_times ma' (ma_sqrt (ma_of_rat (rat_of_nat fact)))" obtain a b where q: "quotient_of (ma_rat ma) = (a,b)" by force obtain sq fact where ppf: "prime_product_factor (nat (abs a * b)) = (sq,fact)" by force define asq where "asq = rat_of_int (sgn a) * of_nat sq / of_int b" define ma' where "ma' = ma_of_rat asq" define sqrt where "sqrt = ma_sqrt (ma_of_rat (rat_of_nat fact))" have num: "?num = ma_times ma' sqrt" unfolding q ppf asq_def Let_def split ma'_def sqrt_def .. let ?inv = "λ ma. ma_coeff ma = 0 ∧ ma_base ma = 0 ∨ ma_coeff ma ≠ 0 ∧ prime_product (ma_base ma)" have ma': "?inv ma'" unfolding ma'_def by (transfer, auto) have id: "⋀ i. int i * 1 = i" "⋀ i :: rat. i / 1 = i" "rat_of_int 1 = 1" "inverse (1 :: rat) = 1" "⋀ n. nat ¦int n¦ = n" by auto from prime_product_factor[OF ppf] have "prime_product fact" by auto hence sqrt: "?inv sqrt" unfolding sqrt_def by (transfer, unfold split quotient_of_nat Let_def id, case_tac "sqrt_int ¦int facta¦", auto) show "?inv ?num" unfolding num using ma' sqrt by (transfer, auto simp: ma_normalize_def split: if_splits) (* slow *) qed lemma sqrt_sgn[simp]: "sqrt (of_int (sgn a)) = of_int (sgn a)" by (cases "a ≥ 0", cases "a = 0", auto simp: real_sqrt_minus) lemma mau_sqrt_main: "mau_coeff r = 0 ⟹ sqrt (real_of_u r) = real_of_u (mau_sqrt r)" proof (transfer) fix r assume "ma_coeff r = 0" hence rr: "real_of r = of_rat (ma_rat r)" by (transfer, auto) obtain a b where q: "quotient_of (ma_rat r) = (a,b)" by force from quotient_of_div[OF q] have r: "ma_rat r = of_int a / of_int b" by auto from quotient_of_denom_pos[OF q] have b: "b > 0" by auto obtain sq fact where ppf: "prime_product_factor (nat (¦a¦ * b)) = (sq, fact)" by force from prime_product_factor[OF ppf] have ab: "nat (¦a¦ * b) = sq * sq * fact" .. have "sqrt (real_of r) = sqrt(of_int a / of_int b)" unfolding rr r by (metis of_rat_divide of_rat_of_int_eq) also have "real_of_int a / of_int b = of_int a * of_int b / (of_int b * of_int b)" using b by auto also have "sqrt (…) = sqrt (of_int a * of_int b) / of_int b" using sqrt_sqrt[of "real_of_int b"] b by (metis less_eq_real_def of_int_0_less_iff real_sqrt_divide real_sqrt_mult) also have "real_of_int a * of_int b = real_of_int (a * b)" by auto also have "a * b = sgn a * (abs a * b)" by (simp, metis mult_sgn_abs) also have "real_of_int (…) = of_int (sgn a) * real_of_int (¦a¦ * b)" unfolding of_int_mult[of "sgn a"] .. also have "real_of_int (¦a¦ * b) = of_nat (nat (abs a * b))" using b by (metis abs_sgn mult_pos_pos mult_zero_left nat_int of_int_of_nat_eq of_nat_0 zero_less_abs_iff zero_less_imp_eq_int) also have "… = of_nat fact * (of_nat sq * of_nat sq)" unfolding ab of_nat_mult by simp also have "sqrt (of_int (sgn a) * (of_nat fact * (of_nat sq * of_nat sq))) = of_int (sgn a) * sqrt (of_nat fact) * of_nat sq" unfolding real_sqrt_mult by simp finally have r: "sqrt (real_of r) = real_of_int (sgn a) * real_of_nat sq / real_of_int b * sqrt (real_of_nat fact)" by simp let ?asqb = "ma_of_rat (rat_of_int (sgn a) * rat_of_nat sq / rat_of_int b)" let ?f = "ma_of_rat (rat_of_nat fact)" let ?sq = "ma_sqrt ?f" have sq: "0 ≤ ma_rat ?f" "ma_coeff ?f = 0" by ((transfer, simp)+) have compat: "⋀ m. (ma_compatible ?asqb m) = True" by (transfer, auto) show "sqrt (real_of r) = real_of (let (a, b) = quotient_of (ma_rat r); (sq, fact) = prime_product_factor (nat (¦a¦ * b)); ma' = ma_of_rat (rat_of_int (sgn a) * rat_of_nat sq / rat_of_int b) in ma_times ma' (ma_sqrt (ma_of_rat (rat_of_nat fact))))" unfolding q ppf Let_def split unfolding r unfolding ma_times[symmetric, of ?asqb, unfolded compat if_True] unfolding ma_sqrt_main[OF sq, symmetric] unfolding ma_of_rat[symmetric] by (simp add: of_rat_divide of_rat_mult) qed lemma mau_sqrt: "sqrt (real_of_u r) = (if mau_coeff r = 0 then real_of_u (mau_sqrt r) else Code.abort (STR ''cannot represent sqrt of irrational number'') (λ _. sqrt (real_of_u r)))" using mau_sqrt_main[of r] by (cases "mau_coeff r = 0", auto) lemma mau_0: "0 = real_of_u (mau_of_rat 0)" using ma_0 by (transfer, auto) lemma mau_1: "1 = real_of_u (mau_of_rat 1)" using ma_1 by (transfer, auto) lemma mau_equal: "HOL.equal (real_of_u r1) (real_of_u r2) = mau_equal r1 r2" unfolding equal_real_def using real_of_u_inj[of r1 r2] by (transfer, transfer, auto) lemma mau_ge_0: "ge_0 (real_of_u x) = mau_ge_0 x" using ma_ge_0 by (transfer, auto) definition real_lt :: "real ⇒ real ⇒ bool" where "real_lt = (<)" text‹The following code equation terminates if it is started on two different inputs.› lemma real_lt [code equation]: "real_lt x y = (let fx = floor x; fy = floor y in (if fx < fy then True else if fx > fy then False else real_lt (x * 1024) (y * 1024)))" proof (cases "floor x < floor y") case True thus ?thesis by (auto simp: real_lt_def floor_less_cancel) next case False note nless = this show ?thesis proof (cases "floor x > floor y") case True from floor_less_cancel[OF this] True nless show ?thesis by (simp add: real_lt_def) next case False with nless show ?thesis unfolding real_lt_def by auto qed qed text ‹For comparisons we first check for equality. Then, if the bases are compatible we can just compare the differences with 0. Otherwise, we start the recursive algorithm @{const real_lt} which works on arbitrary bases. In this way, we have an implementation of comparisons which can compare all representable numbers. Note that in @{theory Real_Impl.Real_Impl} we did not use @{const real_lt} as there the code-equations for equality already require identical bases. › lemma comparison_impl: "real_of_u x ≤ real_of_u y ⟷ real_of_u x = real_of_u y ∨ (if mau_compatible x y then ge_0 (real_of_u y - real_of_u x) else real_lt (real_of_u x) (real_of_u y))" "real_of_u x < real_of_u y ⟷ real_of_u x ≠ real_of_u y ∧ (if mau_compatible x y then ge_0 (real_of_u y - real_of_u x) else real_lt (real_of_u x) (real_of_u y))" unfolding ge_0_def real_lt_def by (auto simp del: real_of_u_inj) lemma mau_is_rat: "is_rat (real_of_u x) = mau_is_rat x" using ma_is_rat by (transfer, auto) lift_definition ma_show_real :: "mini_alg ⇒ string" is "λ (p,q,b). let sb = shows ''sqrt('' ∘ shows b ∘ shows '')''; qb = (if q = 1 then sb else if q = -1 then shows ''-'' ∘ sb else shows q ∘ shows ''*'' ∘ sb) in if q = 0 then shows p [] else if p = 0 then qb [] else if q < 0 then ((shows p ∘ qb) []) else ((shows p ∘ shows ''+'' ∘ qb) [])" . lift_definition mau_show_real :: "mini_alg_unique ⇒ string" is ma_show_real . overloading show_real ≡ show_real begin definition show_real where "show_real x ≡ (if (∃ y. x = real_of_u y) then mau_show_real (THE y. x = real_of_u y) else [])" end lemma mau_show_real: "show_real (real_of_u x) = mau_show_real x" unfolding show_real_def by simp code_datatype real_of_u declare [[code drop: "plus :: real ⇒ real ⇒ real" "uminus :: real ⇒ real" "times :: real ⇒ real ⇒ real" "inverse :: real ⇒ real" "floor :: real ⇒ int" sqrt "HOL.equal :: real ⇒ real ⇒ bool" ge_0 is_rat "less :: real ⇒ real ⇒ bool" "less_eq :: real ⇒ real ⇒ bool" ]] lemmas mau_code_eqns [code] = mau_floor mau_0 mau_1 mau_uminus mau_inverse mau_sqrt mau_plus mau_times mau_equal mau_ge_0 mau_is_rat mau_show_real comparison_impl text ‹Some tests with small numbers. To work on larger number, one should additionally import the theories for efficient calculation on numbers› value "⌊101.1 * (sqrt 18 + 6 * sqrt 0.5)⌋" value "⌊324 * sqrt 7 + 0.001⌋" value "101.1 * (sqrt 18 + 6 * sqrt 0.5) = 324 * sqrt 7 + 0.001" value "101.1 * (sqrt 18 + 6 * sqrt 0.5) > 324 * sqrt 7 + 0.001" value "show (101.1 * (sqrt 18 + 6 * sqrt 0.5))" value "(sqrt 0.1 ∈ ℚ, sqrt (- 0.09) ∈ ℚ)" end