Session Real_Impl

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