(* Author:     Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 *)

 header {* Prove Real Valued Inequalities by Computation *}

 theory Approximation

 imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat

 begin

 section "Horner Scheme"

 subsection {* Define auxiliary helper @{text horner} function *}

 primrec horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where

 "horner F G 0 i k x       = 0" |

 "horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"

 lemma horner_schema': fixes x :: real  and a :: "nat \<Rightarrow> real"

   shows "a 0 - x * (\<Sum> i=0..<n. (-1)^i * a (Suc i) * x^i) = (\<Sum> i=0..<Suc n. (-1)^i * a i * x^i)"

 proof -

   have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto

   show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]

     setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n  *a n * x^n"] by auto

 qed

 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"

   assumes f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

   shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / real (f (j' + j))) * x ^ j)"

 proof (induct n arbitrary: i k j')

   case (Suc n)

   show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]

     using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto

 qed auto

 lemma horner_bounds':

   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

   and lb_0: "\<And> i k x. lb 0 i k x = 0"

   and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"

   and ub_0: "\<And> i k x. ub 0 i k x = 0"

   and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"

   shows "real (lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') (real x) \<and> 

          horner F G n ((F ^^ j') s) (f j') (real x) \<le> real (ub n ((F ^^ j') s) (f j') x)"

   (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")

 proof (induct n arbitrary: j')

   case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto

 next

   case (Suc n)

   have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps real_of_float_sub diff_def

   proof (rule add_mono)

     show "real (lapprox_rat prec 1 (int (f j'))) \<le> 1 / real (f j')" using lapprox_rat[of prec 1  "int (f j')"] by auto

     from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> real x`

     show "- real (x * ub n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x) \<le> - (real x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) (real x))"

       unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)

   qed

   moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps real_of_float_sub diff_def

   proof (rule add_mono)

     show "1 / real (f j') \<le> real (rapprox_rat prec 1 (int (f j')))" using rapprox_rat[of 1 "int (f j')" prec] by auto

     from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> real x`

     show "- (real x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) (real x)) \<le> 

           - real (x * lb n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x)"

       unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)

   qed

   ultimately show ?case by blast

 qed

 subsection "Theorems for floating point functions implementing the horner scheme"

 text {*

 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are

 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.

 *}

 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"

   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

   and lb_0: "\<And> i k x. lb 0 i k x = 0"

   and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"

   and ub_0: "\<And> i k x. ub 0 i k x = 0"

   and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"

   shows "real (lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * real x ^ j)" (is "?lb") and 

     "(\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * (real x ^ j)) \<le> real (ub n ((F ^^ j') s) (f j') x)" (is "?ub")

 proof -

   have "?lb  \<and> ?ub" 

     using horner_bounds'[where lb=lb, OF `0 \<le> real x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]

     unfolding horner_schema[where f=f, OF f_Suc] .

   thus "?lb" and "?ub" by auto

 qed

 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"

   assumes "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

   and lb_0: "\<And> i k x. lb 0 i k x = 0"

   and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)"

   and ub_0: "\<And> i k x. ub 0 i k x = 0"

   and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)"

   shows "real (lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / real (f (j' + j))) * real x ^ j)" (is "?lb") and 

     "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (real x ^ j)) \<le> real (ub n ((F ^^ j') s) (f j') x)" (is "?ub")

 proof -

   { fix x y z :: float have "x - y * z = x + - y * z"

       by (cases x, cases y, cases z, simp add: plus_float.simps minus_float_def uminus_float.simps times_float.simps algebra_simps)

   } note diff_mult_minus = this

   { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this

   have move_minus: "real (-x) = -1 * real x" by auto

   have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * real x ^ j) = 

     (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j)"

   proof (rule setsum_cong, simp)

     fix j assume "j \<in> {0 ..< n}"

     show "1 / real (f (j' + j)) * real x ^ j = -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j"

       unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]

       unfolding real_mult_commute unfolding real_mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]

       by auto

   qed

   have "0 \<le> real (-x)" using assms by auto

   from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec

     and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", unfolded lb_Suc ub_Suc diff_mult_minus,

     OF this f_Suc lb_0 refl ub_0 refl]

   show "?lb" and "?ub" unfolding minus_minus sum_eq

     by auto

 qed

 subsection {* Selectors for next even or odd number *}

 text {*

 The horner scheme computes alternating series. To get the upper and lower bounds we need to

 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.

 *}

 definition get_odd :: "nat \<Rightarrow> nat" where

   "get_odd n = (if odd n then n else (Suc n))"

 definition get_even :: "nat \<Rightarrow> nat" where

   "get_even n = (if even n then n else (Suc n))"

 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)

 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)

 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"

 proof (cases "odd n")

   case True hence "0 < n" by (rule odd_pos)

   from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto

   thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast

 next

   case False hence "odd (Suc n)" by auto

   thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast

 qed

 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .

 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto

 section "Power function"

 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where

 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)

                       else if u < 0         then (u ^ n, l ^ n)

                                             else (0, (max (-l) u) ^ n))"

 lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {real l .. real u}"

   shows "x ^ n \<in> {real l1..real u1}"

 proof (cases "even n")

   case True

   show ?thesis

   proof (cases "0 < l")

     case True hence "odd n \<or> 0 < l" and "0 \<le> real l" unfolding less_float_def by auto

     have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto

     have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using `0 \<le> real l` and assms unfolding atLeastAtMost_iff using power_mono[of "real l" x] power_mono[of x "real u"] by auto

     thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto

   next

     case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto

     show ?thesis

     proof (cases "u < 0")

       case True hence "0 \<le> - real u" and "- real u \<le> - x" and "0 \<le> - x" and "-x \<le> - real l" using assms unfolding less_float_def by auto

       hence "real u ^ n \<le> x ^ n" and "x ^ n \<le> real l ^ n" using power_mono[of  "-x" "-real l" n] power_mono[of "-real u" "-x" n] 

 	unfolding power_minus_even[OF `even n`] by auto

       moreover have u1: "u1 = l ^ n" and l1: "l1 = u ^ n" using assms unfolding float_power_bnds_def if_not_P[OF P] if_P[OF True] by auto

       ultimately show ?thesis using float_power by auto

     next

       case False

       have "\<bar>x\<bar> \<le> real (max (-l) u)"

       proof (cases "-l \<le> u")

 	case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto

       next

 	case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto

       qed

       hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto

       have u1: "u1 = (max (-l) u) ^ n" and l1: "l1 = 0" using assms unfolding float_power_bnds_def if_not_P[OF P] if_not_P[OF False] by auto

       show ?thesis unfolding atLeastAtMost_iff l1 u1 float_power using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto

     qed

   qed

 next

   case False hence "odd n \<or> 0 < l" by auto

   have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto

   have "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto

   thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto

 qed

 lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {real l .. real u} \<longrightarrow> real l1 \<le> x ^ n \<and> x ^ n \<le> real u1"

   using float_power_bnds by auto

 section "Square root"

 text {*

 The square root computation is implemented as newton iteration. As first first step we use the

 nearest power of two greater than the square root.

 *}

 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |

 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x

                                   in Float 1 -1 * (y + float_divr prec x y))"

 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where

 "ub_sqrt prec x = (if 0 < x then (sqrt_iteration prec prec x)

               else if x < 0 then - lb_sqrt prec (- x)

                             else 0)" |

 "lb_sqrt prec x = (if 0 < x then (float_divl prec x (sqrt_iteration prec prec x))

               else if x < 0 then - ub_sqrt prec (- x)

                             else 0)"

 by pat_completeness auto

 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def)

 declare lb_sqrt.simps[simp del]

 declare ub_sqrt.simps[simp del]

 lemma sqrt_ub_pos_pos_1:

   assumes "sqrt x < b" and "0 < b" and "0 < x"

   shows "sqrt x < (b + x / b)/2"

 proof -

   from assms have "0 < (b - sqrt x) ^ 2 " by simp

   also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra

   also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)

   finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption

   hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms

     by (simp add: field_simps power2_eq_square)

   thus ?thesis by (simp add: field_simps)

 qed

 lemma sqrt_iteration_bound: assumes "0 < real x"

   shows "sqrt (real x) < real (sqrt_iteration prec n x)"

 proof (induct n)

   case 0

   show ?case

   proof (cases x)

     case (Float m e)

     hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto

     hence "0 < sqrt (real m)" by auto

     have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto

     have "real x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"

       unfolding pow2_add pow2_int Float real_of_float_simp by auto

     also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"

     proof (rule mult_strict_right_mono, auto)

       show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]

 	unfolding real_of_int_less_iff[of m, symmetric] by auto

     qed

     finally have "sqrt (real x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto

     also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"

     proof -

       let ?E = "e + bitlen m"

       have E_mod_pow: "pow2 (?E mod 2) < 4"

       proof (cases "?E mod 2 = 1")

 	case True thus ?thesis by auto

       next

 	case False

 	have "0 \<le> ?E mod 2" by auto

 	have "?E mod 2 < 2" by auto

 	from this[THEN zless_imp_add1_zle]

 	have "?E mod 2 \<le> 0" using False by auto

 	from xt1(5)[OF `0 \<le> ?E mod 2` this]

 	show ?thesis by auto

       qed

       hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto

       hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto

       have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto

       have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"

 	unfolding E_eq unfolding pow2_add ..

       also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"

 	unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto

       also have "\<dots> < pow2 (?E div 2) * 2"

 	by (rule mult_strict_left_mono, auto intro: E_mod_pow)

       also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto

       finally show ?thesis by auto

     qed

     finally show ?thesis

       unfolding Float sqrt_iteration.simps real_of_float_simp by auto

   qed

 next

   case (Suc n)

   let ?b = "sqrt_iteration prec n x"

   have "0 < sqrt (real x)" using `0 < real x` by auto

   also have "\<dots> < real ?b" using Suc .

   finally have "sqrt (real x) < (real ?b + real x / real ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < real x`] by auto

   also have "\<dots> \<le> (real ?b + real (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)

   also have "\<dots> = real (Float 1 -1) * (real ?b + real (float_divr prec x ?b))" by auto

   finally show ?case unfolding sqrt_iteration.simps Let_def real_of_float_mult real_of_float_add right_distrib .

 qed

 lemma sqrt_iteration_lower_bound: assumes "0 < real x"

   shows "0 < real (sqrt_iteration prec n x)" (is "0 < ?sqrt")

 proof -

   have "0 < sqrt (real x)" using assms by auto

   also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .

   finally show ?thesis .

 qed

 lemma lb_sqrt_lower_bound: assumes "0 \<le> real x"

   shows "0 \<le> real (lb_sqrt prec x)"

 proof (cases "0 < x")

   case True hence "0 < real x" and "0 \<le> x" using `0 \<le> real x` unfolding less_float_def le_float_def by auto

   hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto 

   hence "0 \<le> real (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF `0 \<le> x`] unfolding le_float_def by auto

   thus ?thesis unfolding lb_sqrt.simps using True by auto

 next

   case False with `0 \<le> real x` have "real x = 0" unfolding less_float_def by auto

   thus ?thesis unfolding lb_sqrt.simps less_float_def by auto

 qed

 lemma bnds_sqrt':

   shows "sqrt (real x) \<in> { real (lb_sqrt prec x) .. real (ub_sqrt prec x) }"

 proof -

   { fix x :: float assume "0 < x"

     hence "0 < real x" and "0 \<le> real x" unfolding less_float_def by auto

     hence sqrt_gt0: "0 < sqrt (real x)" by auto

     hence sqrt_ub: "sqrt (real x) < real (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto

     have "real (float_divl prec x (sqrt_iteration prec prec x)) \<le>

           real x / real (sqrt_iteration prec prec x)" by (rule float_divl)

     also have "\<dots> < real x / sqrt (real x)"

       by (rule divide_strict_left_mono[OF sqrt_ub `0 < real x`

                mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])

     also have "\<dots> = sqrt (real x)" 

       unfolding inverse_eq_iff_eq[of _ "sqrt (real x)", symmetric]

 	        sqrt_divide_self_eq[OF `0 \<le> real x`, symmetric] by auto

     finally have "real (lb_sqrt prec x) \<le> sqrt (real x)"

       unfolding lb_sqrt.simps if_P[OF `0 < x`] by auto }

   note lb = this

   { fix x :: float assume "0 < x"

     hence "0 < real x" unfolding less_float_def by auto

     hence "0 < sqrt (real x)" by auto

     hence "sqrt (real x) < real (sqrt_iteration prec prec x)"

       using sqrt_iteration_bound by auto

     hence "sqrt (real x) \<le> real (ub_sqrt prec x)"

       unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }

   note ub = this

   show ?thesis

   proof (cases "0 < x")

     case True with lb ub show ?thesis by auto

   next case False show ?thesis

   proof (cases "real x = 0")

     case True thus ?thesis 

       by (auto simp add: less_float_def lb_sqrt.simps ub_sqrt.simps)

   next

     case False with `\<not> 0 < x` have "x < 0" and "0 < -x"

       by (auto simp add: less_float_def)

     with `\<not> 0 < x`

     show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]

       by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)

   qed qed

 qed

 lemma bnds_sqrt: "\<forall> x lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> sqrt x \<and> sqrt x \<le> real u"

 proof ((rule allI) +, rule impI, erule conjE, rule conjI)

   fix x lx ux

   assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)"

     and x: "x \<in> {real lx .. real ux}"

   hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto

   have "sqrt (real lx) \<le> sqrt x" using x by auto

   from order_trans[OF _ this]

   show "real l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto

   have "sqrt x \<le> sqrt (real ux)" using x by auto

   from order_trans[OF this]

   show "sqrt x \<le> real u" unfolding u using bnds_sqrt'[of ux prec] by auto

 qed

 section "Arcus tangens and \<pi>"

 subsection "Compute arcus tangens series"

 text {*

 As first step we implement the computation of the arcus tangens series. This is only valid in the range

 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.

 *}

 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"

 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

   "ub_arctan_horner prec 0 k x = 0"

 | "ub_arctan_horner prec (Suc n) k x = 

     (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"

 | "lb_arctan_horner prec 0 k x = 0"

 | "lb_arctan_horner prec (Suc n) k x = 

     (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"

 lemma arctan_0_1_bounds': assumes "0 \<le> real x" "real x \<le> 1" and "even n"

   shows "arctan (real x) \<in> {real (x * lb_arctan_horner prec n 1 (x * x)) .. real (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"

 proof -

   let "?c i" = "-1^i * (1 / real (i * 2 + 1) * real x ^ (i * 2 + 1))"

   let "?S n" = "\<Sum> i=0..<n. ?c i"

   have "0 \<le> real (x * x)" by auto

   from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto

   have "arctan (real x) \<in> { ?S n .. ?S (Suc n) }"

   proof (cases "real x = 0")

     case False

     hence "0 < real x" using `0 \<le> real x` by auto

     hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * real x ^ (0 * 2 + 1)" by auto 

     have "\<bar> real x \<bar> \<le> 1"  using `0 \<le> real x` `real x \<le> 1` by auto

     from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`]

     show ?thesis unfolding arctan_series[OF `\<bar> real x \<bar> \<le> 1`] Suc_plus1  .

   qed auto

   note arctan_bounds = this[unfolded atLeastAtMost_iff]

   have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto

   note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 

     and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"

     and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", 

     OF `0 \<le> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]

   { have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"

       using bounds(1) `0 \<le> real x`

       unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

       unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]

       by (auto intro!: mult_left_mono)

     also have "\<dots> \<le> arctan (real x)" using arctan_bounds ..

     finally have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (real x)" . }

   moreover

   { have "arctan (real x) \<le> ?S (Suc n)" using arctan_bounds ..

     also have "\<dots> \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))"

       using bounds(2)[of "Suc n"] `0 \<le> real x`

       unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

       unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]

       by (auto intro!: mult_left_mono)

     finally have "arctan (real x) \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }

   ultimately show ?thesis by auto

 qed

 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"

   shows "arctan (real x) \<in> {real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"

 proof (cases "even n")

   case True

   obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto

   hence "even n'" unfolding even_Suc by auto

   have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"

     unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto

   moreover

   have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"

     unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n`] by auto

   ultimately show ?thesis by auto

 next

   case False hence "0 < n" by (rule odd_pos)

   from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..

   from False[unfolded this even_Suc]

   have "even n'" and "even (Suc (Suc n'))" by auto

   have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .

   have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"

     unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto

   moreover

   have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"

     unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even (Suc (Suc n'))`] by auto

   ultimately show ?thesis by auto

 qed

 subsection "Compute \<pi>"

 definition ub_pi :: "nat \<Rightarrow> float" where

   "ub_pi prec = (let A = rapprox_rat prec 1 5 ; 

                      B = lapprox_rat prec 1 239

                  in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) - 

                                                   B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"

 definition lb_pi :: "nat \<Rightarrow> float" where

   "lb_pi prec = (let A = lapprox_rat prec 1 5 ; 

                      B = rapprox_rat prec 1 239

                  in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) - 

                                                   B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"

 lemma pi_boundaries: "pi \<in> {real (lb_pi n) .. real (ub_pi n)}"

 proof -

   have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto

   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto

     let ?k = "rapprox_rat prec 1 k"

     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto

     have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)

     have "real ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]

       by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)

     have "1 / real k \<le> real ?k" using rapprox_rat[where x=1 and y=k] by auto

     hence "arctan (1 / real k) \<le> arctan (real ?k)" by (rule arctan_monotone')

     also have "\<dots> \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"

       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto

     finally have "arctan (1 / (real k)) \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .

   } note ub_arctan = this

   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto

     let ?k = "lapprox_rat prec 1 k"

     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto

     have "1 / real k \<le> 1" using `1 < k` by auto

     have "\<And>n. 0 \<le> real ?k" using lapprox_rat_bottom[where x=1 and y=k, OF zero_le_one `0 < k`] by (auto simp add: `1 div k = 0`)

     have "\<And>n. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)

     have "real ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto

     have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (real ?k)"

       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto

     also have "\<dots> \<le> arctan (1 / real k)" using `real ?k \<le> 1 / real k` by (rule arctan_monotone')

     finally have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .

   } note lb_arctan = this

   have "pi \<le> real (ub_pi n)"

     unfolding ub_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub unfolding Float_num

     using lb_arctan[of 239] ub_arctan[of 5]

     by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)

   moreover

   have "real (lb_pi n) \<le> pi"

     unfolding lb_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub Float_num

     using lb_arctan[of 5] ub_arctan[of 239]

     by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)

   ultimately show ?thesis by auto

 qed

 subsection "Compute arcus tangens in the entire domain"

 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where

   "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;

                            lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)

     in (if x < 0          then - ub_arctan prec (-x) else

         if x \<le> Float 1 -1 then lb_horner x else

         if x \<le> Float 1 1  then Float 1 1 * lb_horner (float_divl prec x (1 + ub_sqrt prec (1 + x * x)))

                           else (let inv = float_divr prec 1 x

                                 in if inv > 1 then 0

                                               else lb_pi prec * Float 1 -1 - ub_horner inv)))"

 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;

                            ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)

     in (if x < 0          then - lb_arctan prec (-x) else

         if x \<le> Float 1 -1 then ub_horner x else

         if x \<le> Float 1 1  then let y = float_divr prec x (1 + lb_sqrt prec (1 + x * x))

                                in if y > 1 then ub_pi prec * Float 1 -1

                                            else Float 1 1 * ub_horner y

                           else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"

 by pat_completeness auto

 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def)

 declare ub_arctan_horner.simps[simp del]

 declare lb_arctan_horner.simps[simp del]

 lemma lb_arctan_bound': assumes "0 \<le> real x"

   shows "real (lb_arctan prec x) \<le> arctan (real x)"

 proof -

   have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto

   let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"

     and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"

   show ?thesis

   proof (cases "x \<le> Float 1 -1")

     case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto

     show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]

       using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto

   next

     case False hence "0 < real x" unfolding le_float_def Float_num by auto

     let ?R = "1 + sqrt (1 + real x * real x)"

     let ?fR = "1 + ub_sqrt prec (1 + x * x)"

     let ?DIV = "float_divl prec x ?fR"

     have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto

     hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)

     have "sqrt (real (1 + x * x)) \<le> real (ub_sqrt prec (1 + x * x))"

       using bnds_sqrt'[of "1 + x * x"] by auto

     hence "?R \<le> real ?fR" by auto

     hence "0 < ?fR" and "0 < real ?fR" unfolding less_float_def using `0 < ?R` by auto

     have monotone: "real (float_divl prec x ?fR) \<le> real x / ?R"

     proof -

       have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)

       also have "\<dots> \<le> real x / ?R" by (rule divide_left_mono[OF `?R \<le> real ?fR` `0 \<le> real x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> real ?fR`] divisor_gt0]])

       finally show ?thesis .

     qed

     show ?thesis

     proof (cases "x \<le> Float 1 1")

       case True

       have "real x \<le> sqrt (real (1 + x * x))" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto

       also have "\<dots> \<le> real (ub_sqrt prec (1 + x * x))"

 	using bnds_sqrt'[of "1 + x * x"] by auto

       finally have "real x \<le> real ?fR" by auto

       moreover have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)

       ultimately have "real ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < real ?fR`, symmetric] by auto

       have "0 \<le> real ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto

       have "real (Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (real (float_divl prec x ?fR))" unfolding real_of_float_mult[of "Float 1 1"] Float_num

 	using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto

       also have "\<dots> \<le> 2 * arctan (real x / ?R)"

 	using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)

       also have "2 * arctan (real x / ?R) = arctan (real x)" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 . 

       finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF True] .

     next

       case False

       hence "2 < real x" unfolding le_float_def Float_num by auto

       hence "1 \<le> real x" by auto

       let "?invx" = "float_divr prec 1 x"

       have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto

       show ?thesis

       proof (cases "1 < ?invx")

 	case True

 	show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF False] if_P[OF True] 

 	  using `0 \<le> arctan (real x)` by auto

       next

 	case False

 	hence "real ?invx \<le> 1" unfolding less_float_def by auto

 	have "0 \<le> real ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> real x`)

 	have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto

 	have "arctan (1 / real x) \<le> arctan (real ?invx)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divr)

 	also have "\<dots> \<le> real (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto

 	finally have "pi / 2 - real (?ub_horner ?invx) \<le> arctan (real x)"

 	  using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`]

 	  unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto

 	moreover

 	have "real (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding real_of_float_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto

 	ultimately

 	show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]

 	  by auto

       qed

     qed

   qed

 qed

 lemma ub_arctan_bound': assumes "0 \<le> real x"

   shows "arctan (real x) \<le> real (ub_arctan prec x)"

 proof -

   have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto

   let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"

     and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"

   show ?thesis

   proof (cases "x \<le> Float 1 -1")

     case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto

     show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]

       using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto

   next

     case False hence "0 < real x" unfolding le_float_def Float_num by auto

     let ?R = "1 + sqrt (1 + real x * real x)"

     let ?fR = "1 + lb_sqrt prec (1 + x * x)"

     let ?DIV = "float_divr prec x ?fR"

     have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto

     hence "0 \<le> real (1 + x*x)" by auto

     hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)

     have "real (lb_sqrt prec (1 + x * x)) \<le> sqrt (real (1 + x * x))"

       using bnds_sqrt'[of "1 + x * x"] by auto

     hence "real ?fR \<le> ?R" by auto

     have "0 < real ?fR" unfolding real_of_float_add real_of_float_1 by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> real (1 + x*x)`])

     have monotone: "real x / ?R \<le> real (float_divr prec x ?fR)"

     proof -

       from divide_left_mono[OF `real ?fR \<le> ?R` `0 \<le> real x` mult_pos_pos[OF divisor_gt0 `0 < real ?fR`]]

       have "real x / ?R \<le> real x / real ?fR" .

       also have "\<dots> \<le> real ?DIV" by (rule float_divr)

       finally show ?thesis .

     qed

     show ?thesis

     proof (cases "x \<le> Float 1 1")

       case True

       show ?thesis

       proof (cases "?DIV > 1")

 	case True

 	have "pi / 2 \<le> real (ub_pi prec * Float 1 -1)" unfolding real_of_float_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto

 	from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]

 	show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_P[OF True] .

       next

 	case False

 	hence "real ?DIV \<le> 1" unfolding less_float_def by auto

 	have "0 \<le> real x / ?R" using `0 \<le> real x` `0 < ?R` unfolding real_0_le_divide_iff by auto

 	hence "0 \<le> real ?DIV" using monotone by (rule order_trans)

 	have "arctan (real x) = 2 * arctan (real x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .

 	also have "\<dots> \<le> 2 * arctan (real ?DIV)"

 	  using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)

 	also have "\<dots> \<le> real (Float 1 1 * ?ub_horner ?DIV)" unfolding real_of_float_mult[of "Float 1 1"] Float_num

 	  using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto

 	finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] .

       qed

     next

       case False

       hence "2 < real x" unfolding le_float_def Float_num by auto

       hence "1 \<le> real x" by auto

       hence "0 < real x" by auto

       hence "0 < x" unfolding less_float_def by auto

       let "?invx" = "float_divl prec 1 x"

       have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto

       have "real ?invx \<le> 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: `1 \<le> real x` divide_le_eq_1_pos[OF `0 < real x`])

       have "0 \<le> real ?invx" unfolding real_of_float_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`)

       have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto

       have "real (?lb_horner ?invx) \<le> arctan (real ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto

       also have "\<dots> \<le> arctan (1 / real x)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divl)

       finally have "arctan (real x) \<le> pi / 2 - real (?lb_horner ?invx)"

 	using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`] 

 	unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto

       moreover

       have "pi / 2 \<le> real (ub_pi prec * Float 1 -1)" unfolding real_of_float_mult Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto

       ultimately

       show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]

 	by auto

     qed

   qed

 qed

 lemma arctan_boundaries:

   "arctan (real x) \<in> {real (lb_arctan prec x) .. real (ub_arctan prec x)}"

 proof (cases "0 \<le> x")

   case True hence "0 \<le> real x" unfolding le_float_def by auto

   show ?thesis using ub_arctan_bound'[OF `0 \<le> real x`] lb_arctan_bound'[OF `0 \<le> real x`] unfolding atLeastAtMost_iff by auto

 next

   let ?mx = "-x"

   case False hence "x < 0" and "0 \<le> real ?mx" unfolding le_float_def less_float_def by auto

   hence bounds: "real (lb_arctan prec ?mx) \<le> arctan (real ?mx) \<and> arctan (real ?mx) \<le> real (ub_arctan prec ?mx)"

     using ub_arctan_bound'[OF `0 \<le> real ?mx`] lb_arctan_bound'[OF `0 \<le> real ?mx`] by auto

   show ?thesis unfolding real_of_float_minus arctan_minus lb_arctan.simps[where x=x] ub_arctan.simps[where x=x] Let_def if_P[OF `x < 0`]

     unfolding atLeastAtMost_iff using bounds[unfolded real_of_float_minus arctan_minus] by auto

 qed

 lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> arctan x \<and> arctan x \<le> real u"

 proof (rule allI, rule allI, rule allI, rule impI)

   fix x lx ux

   assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {real lx .. real ux}"

   hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto

   { from arctan_boundaries[of lx prec, unfolded l]

     have "real l \<le> arctan (real lx)" by (auto simp del: lb_arctan.simps)

     also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')

     finally have "real l \<le> arctan x" .

   } moreover

   { have "arctan x \<le> arctan (real ux)" using x by (auto intro: arctan_monotone')

     also have "\<dots> \<le> real u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)

     finally have "arctan x \<le> real u" .

   } ultimately show "real l \<le> arctan x \<and> arctan x \<le> real u" ..

 qed

 section "Sinus and Cosinus"

 subsection "Compute the cosinus and sinus series"

 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"

 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

   "ub_sin_cos_aux prec 0 i k x = 0"

 | "ub_sin_cos_aux prec (Suc n) i k x = 

     (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"

 | "lb_sin_cos_aux prec 0 i k x = 0"

 | "lb_sin_cos_aux prec (Suc n) i k x = 

     (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"

 lemma cos_aux:

   shows "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (real x)^(2 * i))" (is "?lb")

   and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (real x)^(2 * i)) \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")

 proof -

   have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto

   let "?f n" = "fact (2 * n)"

   { fix n 

     have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)

     have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ n) 1 + 1)"

       unfolding F by auto } note f_eq = this

   from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 

     OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]

   show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"])

 qed

 lemma cos_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"

   shows "cos (real x) \<in> {real (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. real (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"

 proof (cases "real x = 0")

   case False hence "real x \<noteq> 0" by auto

   hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto

   have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0

     using mult_pos_pos[where a="real x" and b="real x"] by auto

   { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^ (2 * i))

     = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")

   proof -

     have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto

     also have "\<dots> = 

       (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto

     also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"

       unfolding sum_split_even_odd ..

     also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"

       by (rule setsum_cong2) auto

     finally show ?thesis by assumption

   qed } note morph_to_if_power = this

   { fix n :: nat assume "0 < n"

     hence "0 < 2 * n" by auto

     obtain t where "0 < t" and "t < real x" and

       cos_eq: "cos (real x) = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * (real x) ^ i) 

       + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (real x)^(2*n)" 

       (is "_ = ?SUM + ?rest / ?fact * ?pow")

       using Maclaurin_cos_expansion2[OF `0 < real x` `0 < 2 * n`] by auto

     have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto

     also have "\<dots> = cos (t + real n * pi)"  using cos_add by auto

     also have "\<dots> = ?rest" by auto

     finally have "cos t * -1^n = ?rest" .

     moreover

     have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto

     hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto

     ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto

     have "0 < ?fact" by auto

     have "0 < ?pow" using `0 < real x` by auto

     {

       assume "even n"

       have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"

 	unfolding morph_to_if_power[symmetric] using cos_aux by auto 

       also have "\<dots> \<le> cos (real x)"

       proof -

 	from even[OF `even n`] `0 < ?fact` `0 < ?pow`

 	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

 	thus ?thesis unfolding cos_eq by auto

       qed

       finally have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (real x)" .

     } note lb = this

     {

       assume "odd n"

       have "cos (real x) \<le> ?SUM"

       proof -

 	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]

 	have "0 \<le> (- ?rest) / ?fact * ?pow"

 	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

 	thus ?thesis unfolding cos_eq by auto

       qed

       also have "\<dots> \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))"

 	unfolding morph_to_if_power[symmetric] using cos_aux by auto

       finally have "cos (real x) \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))" .

     } note ub = this and lb

   } note ub = this(1) and lb = this(2)

   have "cos (real x) \<le> real (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .

   moreover have "real (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (real x)" 

   proof (cases "0 < get_even n")

     case True show ?thesis using lb[OF True get_even] .

   next

     case False

     hence "get_even n = 0" by auto

     have "- (pi / 2) \<le> real x" by (rule order_trans[OF _ `0 < real x`[THEN less_imp_le]], auto)

     with `real x \<le> pi / 2`

     show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps real_of_float_minus real_of_float_0 using cos_ge_zero by auto

   qed

   ultimately show ?thesis by auto

 next

   case True

   show ?thesis

   proof (cases "n = 0")

     case True 

     thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `real x = 0` lapprox_rat[where x="-1" and y=1] by auto

   next

     case False with not0_implies_Suc obtain m where "n = Suc m" by blast

     thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `real x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)

   qed

 qed

 lemma sin_aux: assumes "0 \<le> real x"

   shows "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (real x)^(2 * i + 1))" (is "?lb")

   and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (real x)^(2 * i + 1)) \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")

 proof -

   have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto

   let "?f n" = "fact (2 * n + 1)"

   { fix n 

     have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)

     have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ n) 2 + 1)"

       unfolding F by auto } note f_eq = this

   from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,

     OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]

   show "?lb" and "?ub" using `0 \<le> real x` unfolding real_of_float_mult

     unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

     unfolding real_mult_commute

     by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real x"])

 qed

 lemma sin_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"

   shows "sin (real x) \<in> {real (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. real (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"

 proof (cases "real x = 0")

   case False hence "real x \<noteq> 0" by auto

   hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto

   have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0

     using mult_pos_pos[where a="real x" and b="real x"] by auto

   { fix x n have "(\<Sum> j = 0 ..< n. -1 ^ (((2 * j + 1) - Suc 0) div 2) / (real (fact (2 * j + 1))) * x ^(2 * j + 1))

     = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")

     proof -

       have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto

       have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto

       also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i)) * x ^ i)"

 	unfolding sum_split_even_odd ..

       also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i))) * x ^ i)"

 	by (rule setsum_cong2) auto

       finally show ?thesis by assumption

     qed } note setsum_morph = this

   { fix n :: nat assume "0 < n"

     hence "0 < 2 * n + 1" by auto

     obtain t where "0 < t" and "t < real x" and

       sin_eq: "sin (real x) = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i) 

       + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (real x)^(2*n + 1)" 

       (is "_ = ?SUM + ?rest / ?fact * ?pow")

       using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < real x`] by auto

     have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto

     moreover

     have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto

     hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto

     ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto

     have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)

     have "0 < ?pow" using `0 < real x` by (rule zero_less_power)

     {

       assume "even n"

       have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> 

             (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"

 	using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto

       also have "\<dots> \<le> ?SUM" by auto

       also have "\<dots> \<le> sin (real x)"

       proof -

 	from even[OF `even n`] `0 < ?fact` `0 < ?pow`

 	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

 	thus ?thesis unfolding sin_eq by auto

       qed

       finally have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (real x)" .

     } note lb = this

     {

       assume "odd n"

       have "sin (real x) \<le> ?SUM"

       proof -

 	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]

 	have "0 \<le> (- ?rest) / ?fact * ?pow"

 	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

 	thus ?thesis unfolding sin_eq by auto

       qed

       also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"

 	 by auto

       also have "\<dots> \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" 

 	using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto

       finally have "sin (real x) \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" .

     } note ub = this and lb

   } note ub = this(1) and lb = this(2)

   have "sin (real x) \<le> real (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .

   moreover have "real (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (real x)" 

   proof (cases "0 < get_even n")

     case True show ?thesis using lb[OF True get_even] .

   next

     case False

     hence "get_even n = 0" by auto

     with `real x \<le> pi / 2` `0 \<le> real x`

     show ?thesis unfolding `get_even n = 0` ub_sin_cos_aux.simps real_of_float_minus real_of_float_0 using sin_ge_zero by auto

   qed

   ultimately show ?thesis by auto

 next

   case True

   show ?thesis

   proof (cases "n = 0")

     case True 

     thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `real x = 0` lapprox_rat[where x="-1" and y=1] by auto

   next

     case False with not0_implies_Suc obtain m where "n = Suc m" by blast

     thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `real x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)

   qed

 qed

 subsection "Compute the cosinus in the entire domain"

 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where

 "lb_cos prec x = (let

     horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;

     half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1

   in if x < Float 1 -1 then horner x

 else if x < 1          then half (horner (x * Float 1 -1))

                        else half (half (horner (x * Float 1 -2))))"

 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where

 "ub_cos prec x = (let

     horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;

     half = \<lambda> x. Float 1 1 * x * x - 1

   in if x < Float 1 -1 then horner x

 else if x < 1          then half (horner (x * Float 1 -1))

                        else half (half (horner (x * Float 1 -2))))"

 lemma lb_cos: assumes "0 \<le> real x" and "real x \<le> pi"

   shows "cos (real x) \<in> {real (lb_cos prec x) .. real (ub_cos prec x)}" (is "?cos x \<in> { real (?lb x) .. real (?ub x) }")

 proof -

   { fix x :: real

     have "cos x = cos (x / 2 + x / 2)" by auto

     also have "\<dots> = cos (x / 2) * cos (x / 2) + sin (x / 2) * sin (x / 2) - sin (x / 2) * sin (x / 2) + cos (x / 2) * cos (x / 2) - 1"

       unfolding cos_add by auto

     also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra

     finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .

   } note x_half = this[symmetric]

   have "\<not> x < 0" using `0 \<le> real x` unfolding less_float_def by auto

   let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"

   let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"

   let "?ub_half x" = "Float 1 1 * x * x - 1"

   let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"

   show ?thesis

   proof (cases "x < Float 1 -1")

     case True hence "real x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto

     show ?thesis unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_P[OF `x < Float 1 -1`] Let_def

       using cos_boundaries[OF `0 \<le> real x` `real x \<le> pi / 2`] .

   next

     case False

     { fix y x :: float let ?x2 = "real (x * Float 1 -1)"

       assume "real y \<le> cos ?x2" and "-pi \<le> real x" and "real x \<le> pi"

       hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto

       hence "0 \<le> cos ?x2" by (rule cos_ge_zero)

       have "real (?lb_half y) \<le> cos (real x)"

       proof (cases "y < 0")

 	case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto

       next

 	case False

 	hence "0 \<le> real y" unfolding less_float_def by auto

 	from mult_mono[OF `real y \<le> cos ?x2` `real y \<le> cos ?x2` `0 \<le> cos ?x2` this]

 	have "real y * real y \<le> cos ?x2 * cos ?x2" .

 	hence "2 * real y * real y \<le> 2 * cos ?x2 * cos ?x2" by auto

 	hence "2 * real y * real y - 1 \<le> 2 * cos (real x / 2) * cos (real x / 2) - 1" unfolding Float_num real_of_float_mult by auto

 	thus ?thesis unfolding if_not_P[OF False] x_half Float_num real_of_float_mult real_of_float_sub by auto

       qed

     } note lb_half = this

     { fix y x :: float let ?x2 = "real (x * Float 1 -1)"

       assume ub: "cos ?x2 \<le> real y" and "- pi \<le> real x" and "real x \<le> pi"

       hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto

       hence "0 \<le> cos ?x2" by (rule cos_ge_zero)

       have "cos (real x) \<le> real (?ub_half y)"

       proof -

 	have "0 \<le> real y" using `0 \<le> cos ?x2` ub by (rule order_trans)

 	from mult_mono[OF ub ub this `0 \<le> cos ?x2`]

 	have "cos ?x2 * cos ?x2 \<le> real y * real y" .

 	hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real y * real y" by auto

 	hence "2 * cos (real x / 2) * cos (real x / 2) - 1 \<le> 2 * real y * real y - 1" unfolding Float_num real_of_float_mult by auto

 	thus ?thesis unfolding x_half real_of_float_mult Float_num real_of_float_sub by auto

       qed

     } note ub_half = this

     let ?x2 = "x * Float 1 -1"

     let ?x4 = "x * Float 1 -1 * Float 1 -1"

     have "-pi \<le> real x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] `0 \<le> real x` by (rule order_trans)

     show ?thesis

     proof (cases "x < 1")

       case True hence "real x \<le> 1" unfolding less_float_def by auto

       have "0 \<le> real ?x2" and "real ?x2 \<le> pi / 2" using pi_ge_two `0 \<le> real x` unfolding real_of_float_mult Float_num using assms by auto

       from cos_boundaries[OF this]

       have lb: "real (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> real (?ub_horner ?x2)" by auto

       have "real (?lb x) \<le> ?cos x"

       proof -

 	from lb_half[OF lb `-pi \<le> real x` `real x \<le> pi`]

 	show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto

       qed

       moreover have "?cos x \<le> real (?ub x)"

       proof -

 	from ub_half[OF ub `-pi \<le> real x` `real x \<le> pi`]

 	show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto 

       qed

       ultimately show ?thesis by auto

     next

       case False

       have "0 \<le> real ?x4" and "real ?x4 \<le> pi / 2" using pi_ge_two `0 \<le> real x` `real x \<le> pi` unfolding real_of_float_mult Float_num by auto

       from cos_boundaries[OF this]

       have lb: "real (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> real (?ub_horner ?x4)" by auto

       have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)

       have "real (?lb x) \<le> ?cos x"

       proof -

 	have "-pi \<le> real ?x2" and "real ?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` `real x \<le> pi` by auto

 	from lb_half[OF lb_half[OF lb this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]

 	show ?thesis unfolding lb_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .

       qed

       moreover have "?cos x \<le> real (?ub x)"

       proof -

 	have "-pi \<le> real ?x2" and "real ?x2 \<le> pi" unfolding real_of_float_mult Float_num using pi_ge_two `0 \<le> real x` `real x \<le> pi` by auto

 	from ub_half[OF ub_half[OF ub this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]

 	show ?thesis unfolding ub_cos_def[where x=x] if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x < Float 1 -1`] if_not_P[OF `\<not> x < 1`] Let_def .

       qed

       ultimately show ?thesis by auto

     qed

   qed

 qed

 lemma lb_cos_minus: assumes "-pi \<le> real x" and "real x \<le> 0"

   shows "cos (real (-x)) \<in> {real (lb_cos prec (-x)) .. real (ub_cos prec (-x))}"

 proof -

   have "0 \<le> real (-x)" and "real (-x) \<le> pi" using `-pi \<le> real x` `real x \<le> 0` by auto

   from lb_cos[OF this] show ?thesis .

 qed

 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where

 "bnds_cos prec lx ux = (let

     lpi = round_down prec (lb_pi prec) ;

     upi = round_up prec (ub_pi prec) ;

     k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ;

     lx = lx - k * 2 * (if k < 0 then lpi else upi) ;

     ux = ux - k * 2 * (if k < 0 then upi else lpi)

   in   if - lpi \<le> lx \<and> ux \<le> 0    then (lb_cos prec (-lx), ub_cos prec (-ux))

   else if 0 \<le> lx \<and> ux \<le> lpi      then (lb_cos prec ux, ub_cos prec lx)

   else if - lpi \<le> lx \<and> ux \<le> lpi  then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0)

   else if 0 \<le> lx \<and> ux \<le> 2 * lpi  then (Float -1 0, max (ub_cos prec lx) (ub_cos prec (- (ux - 2 * lpi))))

                                  else (Float -1 0, Float 1 0))"

 lemma floor_int:

   obtains k :: int where "real k = real (floor_fl f)"

 proof -

   assume *: "\<And> k :: int. real k = real (floor_fl f) \<Longrightarrow> thesis"

   obtain m e where fl: "Float m e = floor_fl f" by (cases "floor_fl f", auto)

   from floor_pos_exp[OF this]

   have "real (m* 2^(nat e)) = real (floor_fl f)"

     by (auto simp add: fl[symmetric] real_of_float_def pow2_def)

   from *[OF this] show thesis by blast

 qed

 lemma float_remove_real_numeral[simp]: "real (number_of k :: float) = number_of k"

 proof -

   have "real (number_of k :: float) = real k"

     unfolding number_of_float_def real_of_float_def pow2_def by auto

   also have "\<dots> = real (number_of k :: int)"

     by (simp add: number_of_is_id)

   finally show ?thesis by auto

 qed

 lemma cos_periodic_nat[simp]: fixes n :: nat shows "cos (x + real n * 2 * pi) = cos x"

 proof (induct n arbitrary: x)

   case (Suc n)

   have split_pi_off: "x + real (Suc n) * 2 * pi = (x + real n * 2 * pi) + 2 * pi"

     unfolding Suc_plus1 real_of_nat_add real_of_one real_add_mult_distrib by auto

   show ?case unfolding split_pi_off using Suc by auto

 qed auto

 lemma cos_periodic_int[simp]: fixes i :: int shows "cos (x + real i * 2 * pi) = cos x"

 proof (cases "0 \<le> i")

   case True hence i_nat: "real i = real (nat i)" by auto

   show ?thesis unfolding i_nat by auto

 next

   case False hence i_nat: "real i = - real (nat (-i))" by auto

   have "cos x = cos (x + real i * 2 * pi - real i * 2 * pi)" by auto

   also have "\<dots> = cos (x + real i * 2 * pi)"

     unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat)

   finally show ?thesis by auto

 qed

 lemma bnds_cos: "\<forall> x lx ux. (l, u) = bnds_cos prec lx ux \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> cos x \<and> cos x \<le> real u"

 proof ((rule allI | rule impI | erule conjE) +)

   fix x lx ux

   assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {real lx .. real ux}"

   let ?lpi = "round_down prec (lb_pi prec)"

   let ?upi = "round_up prec (ub_pi prec)"

   let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))"

   let ?lx = "lx - ?k * 2 * (if ?k < 0 then ?lpi else ?upi)"

   let ?ux = "ux - ?k * 2 * (if ?k < 0 then ?upi else ?lpi)"

   obtain k :: int where k: "real k = real ?k" using floor_int .

   have upi: "pi \<le> real ?upi" and lpi: "real ?lpi \<le> pi"

     using round_up[of "ub_pi prec" prec] pi_boundaries[of prec]

           round_down[of prec "lb_pi prec"] by auto

   hence "real ?lx \<le> x - real k * 2 * pi \<and> x - real k * 2 * pi \<le> real ?ux"

     using x by (cases "k = 0") (auto intro!: add_mono

                 simp add: real_diff_def k[symmetric] less_float_def)

   note lx = this[THEN conjunct1] and ux = this[THEN conjunct2]

   hence lx_less_ux: "real ?lx \<le> real ?ux" by (rule order_trans)

   { assume "- ?lpi \<le> ?lx" and x_le_0: "x - real k * 2 * pi \<le> 0"

     with lpi[THEN le_imp_neg_le] lx

     have pi_lx: "- pi \<le> real ?lx" and lx_0: "real ?lx \<le> 0"

       by (simp_all add: le_float_def)

     have "real (lb_cos prec (- ?lx)) \<le> cos (real (- ?lx))"

       using lb_cos_minus[OF pi_lx lx_0] by simp

     also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"

       using cos_monotone_minus_pi_0'[OF pi_lx lx x_le_0]

       by (simp only: real_of_float_minus real_of_int_minus

 	cos_minus real_diff_def mult_minus_left)

     finally have "real (lb_cos prec (- ?lx)) \<le> cos x"

       unfolding cos_periodic_int . }

   note negative_lx = this

   { assume "0 \<le> ?lx" and pi_x: "x - real k * 2 * pi \<le> pi"

     with lx

     have pi_lx: "real ?lx \<le> pi" and lx_0: "0 \<le> real ?lx"

       by (auto simp add: le_float_def)

     have "cos (x + real (-k) * 2 * pi) \<le> cos (real ?lx)"

       using cos_monotone_0_pi'[OF lx_0 lx pi_x]

       by (simp only: real_of_float_minus real_of_int_minus

 	cos_minus real_diff_def mult_minus_left)

     also have "\<dots> \<le> real (ub_cos prec ?lx)"

       using lb_cos[OF lx_0 pi_lx] by simp

     finally have "cos x \<le> real (ub_cos prec ?lx)"

       unfolding cos_periodic_int . }

   note positive_lx = this

   { assume pi_x: "- pi \<le> x - real k * 2 * pi" and "?ux \<le> 0"

     with ux

     have pi_ux: "- pi \<le> real ?ux" and ux_0: "real ?ux \<le> 0"

       by (simp_all add: le_float_def)

     have "cos (x + real (-k) * 2 * pi) \<le> cos (real (- ?ux))"

       using cos_monotone_minus_pi_0'[OF pi_x ux ux_0]

       by (simp only: real_of_float_minus real_of_int_minus

 	  cos_minus real_diff_def mult_minus_left)

     also have "\<dots> \<le> real (ub_cos prec (- ?ux))"

       using lb_cos_minus[OF pi_ux ux_0, of prec] by simp

     finally have "cos x \<le> real (ub_cos prec (- ?ux))"

       unfolding cos_periodic_int . }

   note negative_ux = this

   { assume "?ux \<le> ?lpi" and x_ge_0: "0 \<le> x - real k * 2 * pi"

     with lpi ux

     have pi_ux: "real ?ux \<le> pi" and ux_0: "0 \<le> real ?ux"

       by (simp_all add: le_float_def)

     have "real (lb_cos prec ?ux) \<le> cos (real ?ux)"

       using lb_cos[OF ux_0 pi_ux] by simp

     also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"

       using cos_monotone_0_pi'[OF x_ge_0 ux pi_ux]

       by (simp only: real_of_float_minus real_of_int_minus

 	cos_minus real_diff_def mult_minus_left)

     finally have "real (lb_cos prec ?ux) \<le> cos x"

       unfolding cos_periodic_int . }

   note positive_ux = this

   show "real l \<le> cos x \<and> cos x \<le> real u"

   proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> 0")

     case True with bnds

     have l: "l = lb_cos prec (-?lx)"

       and u: "u = ub_cos prec (-?ux)"

       by (auto simp add: bnds_cos_def Let_def)

     from True lpi[THEN le_imp_neg_le] lx ux

     have "- pi \<le> x - real k * 2 * pi"

       and "x - real k * 2 * pi \<le> 0"

       by (auto simp add: le_float_def)

     with True negative_ux negative_lx

     show ?thesis unfolding l u by simp

   next case False note 1 = this show ?thesis

   proof (cases "0 \<le> ?lx \<and> ?ux \<le> ?lpi")

     case True with bnds 1

     have l: "l = lb_cos prec ?ux"

       and u: "u = ub_cos prec ?lx"

       by (auto simp add: bnds_cos_def Let_def)

     from True lpi lx ux

     have "0 \<le> x - real k * 2 * pi"

       and "x - real k * 2 * pi \<le> pi"

       by (auto simp add: le_float_def)

     with True positive_ux positive_lx

     show ?thesis unfolding l u by simp

   next case False note 2 = this show ?thesis

   proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> ?lpi")

     case True note Cond = this with bnds 1 2

     have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)"

       and u: "u = Float 1 0"

       by (auto simp add: bnds_cos_def Let_def)

     show ?thesis unfolding u l using negative_lx positive_ux Cond

       by (cases "x - real k * 2 * pi < 0", simp_all add: real_of_float_min)

   next case False note 3 = this show ?thesis

   proof (cases "0 \<le> ?lx \<and> ?ux \<le> 2 * ?lpi")

     case True note Cond = this with bnds 1 2 3

     have l: "l = Float -1 0"

       and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))"

       by (auto simp add: bnds_cos_def Let_def)

     have "cos x \<le> real u"

     proof (cases "x - real k * 2 * pi < pi")

       case True hence "x - real k * 2 * pi \<le> pi" by simp

       from positive_lx[OF Cond[THEN conjunct1] this]

       show ?thesis unfolding u by (simp add: real_of_float_max)

     next

       case False hence "pi \<le> x - real k * 2 * pi" by simp

       hence pi_x: "- pi \<le> x - real k * 2 * pi - 2 * pi" by simp

       have "real ?ux \<le> 2 * pi" using Cond lpi by (auto simp add: le_float_def)

       hence "x - real k * 2 * pi - 2 * pi \<le> 0" using ux by simp

       have ux_0: "real (?ux - 2 * ?lpi) \<le> 0"

 	using Cond by (auto simp add: le_float_def)

       from 2 and Cond have "\<not> ?ux \<le> ?lpi" by auto

       hence "- ?lpi \<le> ?ux - 2 * ?lpi" by (auto simp add: le_float_def)

       hence pi_ux: "- pi \<le> real (?ux - 2 * ?lpi)"

 	using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)

       have x_le_ux: "x - real k * 2 * pi - 2 * pi \<le> real (?ux - 2 * ?lpi)"

 	using ux lpi by auto

       have "cos x = cos (x + real (-k) * 2 * pi + real (-1 :: int) * 2 * pi)"

 	unfolding cos_periodic_int ..

       also have "\<dots> \<le> cos (real (?ux - 2 * ?lpi))"

 	using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0]

 	by (simp only: real_of_float_minus real_of_int_minus real_of_one

 	    number_of_Min real_diff_def mult_minus_left real_mult_1)

       also have "\<dots> = cos (real (- (?ux - 2 * ?lpi)))"

 	unfolding real_of_float_minus cos_minus ..

       also have "\<dots> \<le> real (ub_cos prec (- (?ux - 2 * ?lpi)))"

 	using lb_cos_minus[OF pi_ux ux_0] by simp

       finally show ?thesis unfolding u by (simp add: real_of_float_max)

     qed

     thus ?thesis unfolding l by auto

   next

     case False with bnds 1 2 3 show ?thesis by (auto simp add: bnds_cos_def Let_def)

   qed qed qed qed

 qed

 section "Exponential function"

 subsection "Compute the series of the exponential function"

 fun ub_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and lb_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

 "ub_exp_horner prec 0 i k x       = 0" |

 "ub_exp_horner prec (Suc n) i k x = rapprox_rat prec 1 (int k) + x * lb_exp_horner prec n (i + 1) (k * i) x" |

 "lb_exp_horner prec 0 i k x       = 0" |

 "lb_exp_horner prec (Suc n) i k x = lapprox_rat prec 1 (int k) + x * ub_exp_horner prec n (i + 1) (k * i) x"

 lemma bnds_exp_horner: assumes "real x \<le> 0"

   shows "exp (real x) \<in> { real (lb_exp_horner prec (get_even n) 1 1 x) .. real (ub_exp_horner prec (get_odd n) 1 1 x) }"

 proof -

   { fix n

     have F: "\<And> m. ((\<lambda>i. i + 1) ^^ n) m = n + m" by (induct n, auto)

     have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^^ n) 1" unfolding F by auto } note f_eq = this

   note bounds = horner_bounds_nonpos[where f="fact" and lb="lb_exp_horner prec" and ub="ub_exp_horner prec" and j'=0 and s=1,

     OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]

   { have "real (lb_exp_horner prec (get_even n) 1 1 x) \<le> (\<Sum>j = 0..<get_even n. 1 / real (fact j) * real x ^ j)"

       using bounds(1) by auto

     also have "\<dots> \<le> exp (real x)"

     proof -

       obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp (real x) = (\<Sum>m = 0..<get_even n. (real x) ^ m / real (fact m)) + exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"

 	using Maclaurin_exp_le by blast

       moreover have "0 \<le> exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"

 	by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)

       ultimately show ?thesis

 	using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)

     qed

     finally have "real (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (real x)" .

   } moreover

   { 

     have x_less_zero: "real x ^ get_odd n \<le> 0"

     proof (cases "real x = 0")

       case True

       have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto

       thus ?thesis unfolding True power_0_left by auto

     next

       case False hence "real x < 0" using `real x \<le> 0` by auto

       show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `real x < 0`)

     qed

     obtain t where "\<bar>t\<bar> \<le> \<bar>real x\<bar>" and "exp (real x) = (\<Sum>m = 0..<get_odd n. (real x) ^ m / real (fact m)) + exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n)"

       using Maclaurin_exp_le by blast

     moreover have "exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n) \<le> 0"

       by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)

     ultimately have "exp (real x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * real x ^ j)"

       using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)

     also have "\<dots> \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)"

       using bounds(2) by auto

     finally have "exp (real x) \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)" .

   } ultimately show ?thesis by auto

 qed

 subsection "Compute the exponential function on the entire domain"

 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where

 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))

              else let 

                 horner = (\<lambda> x. let  y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x  in if y \<le> 0 then Float 1 -2 else y)

              in if x < - 1 then (case floor_fl x of (Float m e) \<Rightarrow> (horner (float_divl prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))

                            else horner x)" |

 "ub_exp prec x = (if 0 < x    then float_divr prec 1 (lb_exp prec (-x))

              else if x < - 1  then (case floor_fl x of (Float m e) \<Rightarrow> 

                                     (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))

                               else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"

 by pat_completeness auto

 termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if 0 < x then 1 else 0))", auto simp add: less_float_def)

 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"

 proof -

   have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto

   have "1 / 4 = real (Float 1 -2)" unfolding Float_num by auto

   also have "\<dots> \<le> real (lb_exp_horner 1 (get_even 4) 1 1 (- 1))"

     unfolding get_even_def eq4 

     by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)

   also have "\<dots> \<le> exp (real (- 1 :: float))" using bnds_exp_horner[where x="- 1"] by auto

   finally show ?thesis unfolding real_of_float_minus real_of_float_1 . 

 qed

 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"

 proof -

   let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"

   let "?horner x" = "let  y = ?lb_horner x  in if y \<le> 0 then Float 1 -2 else y"

   have pos_horner: "\<And> x. 0 < ?horner x" unfolding Let_def by (cases "?lb_horner x \<le> 0", auto simp add: le_float_def less_float_def)

   moreover { fix x :: float fix num :: nat

     have "0 < real (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def real_of_float_0] by (rule zero_less_power)

     also have "\<dots> = real ((?horner x) ^ num)" using float_power by auto

     finally have "0 < real ((?horner x) ^ num)" .

   }

   ultimately show ?thesis

     unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def

     by (cases "floor_fl x", cases "x < - 1", auto simp add: float_power le_float_def less_float_def)

 qed

 lemma exp_boundaries': assumes "x \<le> 0"

   shows "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"

 proof -

   let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"

   let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"

   have "real x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto

   show ?thesis

   proof (cases "x < - 1")

     case False hence "- 1 \<le> real x" unfolding less_float_def by auto

     show ?thesis

     proof (cases "?lb_exp_horner x \<le> 0")

       from `\<not> x < - 1` have "- 1 \<le> real x" unfolding less_float_def by auto

       hence "exp (- 1) \<le> exp (real x)" unfolding exp_le_cancel_iff .

       from order_trans[OF exp_m1_ge_quarter this]

       have "real (Float 1 -2) \<le> exp (real x)" unfolding Float_num .

       moreover case True

       ultimately show ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto

     next

       case False thus ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)

     qed

   next

     case True

     obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)

     let ?num = "nat (- m) * 2 ^ nat e"

     have "real (floor_fl x) < - 1" using floor_fl `x < - 1` unfolding le_float_def less_float_def real_of_float_minus real_of_float_1 by (rule order_le_less_trans)

     hence "real (floor_fl x) < 0" unfolding Float_floor real_of_float_simp using zero_less_pow2[of xe] by auto

     hence "m < 0"

       unfolding less_float_def real_of_float_0 Float_floor real_of_float_simp

       unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded real_mult_commute] by auto

     hence "1 \<le> - m" by auto

     hence "0 < nat (- m)" by auto

     moreover

     have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto

     hence "(0::nat) < 2 ^ nat e" by auto

     ultimately have "0 < ?num"  by auto

     hence "real ?num \<noteq> 0" by auto

     have e_nat: "int (nat e) = e" using `0 \<le> e` by auto

     have num_eq: "real ?num = real (- floor_fl x)" using `0 < nat (- m)`

       unfolding Float_floor real_of_float_minus real_of_float_simp real_of_nat_mult pow2_int[of "nat e", unfolded e_nat] realpow_real_of_nat[symmetric] by auto

     have "0 < - floor_fl x" using `0 < ?num`[unfolded real_of_nat_less_iff[symmetric]] unfolding less_float_def num_eq[symmetric] real_of_float_0 real_of_nat_zero .

     hence "real (floor_fl x) < 0" unfolding less_float_def by auto

     have "exp (real x) \<le> real (ub_exp prec x)"

     proof -

       have div_less_zero: "real (float_divr prec x (- floor_fl x)) \<le> 0" 

 	using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def real_of_float_0 .

       have "exp (real x) = exp (real ?num * (real x / real ?num))" using `real ?num \<noteq> 0` by auto

       also have "\<dots> = exp (real x / real ?num) ^ ?num" unfolding exp_real_of_nat_mult ..

       also have "\<dots> \<le> exp (real (float_divr prec x (- floor_fl x))) ^ ?num" unfolding num_eq

 	by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto

       also have "\<dots> \<le> real ((?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num)" unfolding float_power

 	by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)

       finally show ?thesis unfolding ub_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def .

     qed

     moreover 

     have "real (lb_exp prec x) \<le> exp (real x)"

     proof -

       let ?divl = "float_divl prec x (- Float m e)"

       let ?horner = "?lb_exp_horner ?divl"

       show ?thesis

       proof (cases "?horner \<le> 0")

 	case False hence "0 \<le> real ?horner" unfolding le_float_def by auto

 	have div_less_zero: "real (float_divl prec x (- floor_fl x)) \<le> 0"

 	  using `real (floor_fl x) < 0` `real x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)

 	have "real ((?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num) \<le>  

           exp (real (float_divl prec x (- floor_fl x))) ^ ?num" unfolding float_power 

 	  using `0 \<le> real ?horner`[unfolded Float_floor[symmetric]] bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] by (auto intro!: power_mono)

 	also have "\<dots> \<le> exp (real x / real ?num) ^ ?num" unfolding num_eq

 	  using float_divl by (auto intro!: power_mono simp del: real_of_float_minus)

 	also have "\<dots> = exp (real ?num * (real x / real ?num))" unfolding exp_real_of_nat_mult ..

 	also have "\<dots> = exp (real x)" using `real ?num \<noteq> 0` by auto

 	finally show ?thesis

 	  unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def if_not_P[OF False] by auto

       next

 	case True

 	have "real (floor_fl x) \<noteq> 0" and "real (floor_fl x) \<le> 0" using `real (floor_fl x) < 0` by auto

 	from divide_right_mono_neg[OF floor_fl[of x] `real (floor_fl x) \<le> 0`, unfolded divide_self[OF `real (floor_fl x) \<noteq> 0`]]

 	have "- 1 \<le> real x / real (- floor_fl x)" unfolding real_of_float_minus by auto

 	from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]

 	have "real (Float 1 -2) \<le> exp (real x / real (- floor_fl x))" unfolding Float_num .

 	hence "real (Float 1 -2) ^ ?num \<le> exp (real x / real (- floor_fl x)) ^ ?num"

 	  by (auto intro!: power_mono simp add: Float_num)

 	also have "\<dots> = exp (real x)" unfolding num_eq exp_real_of_nat_mult[symmetric] using `real (floor_fl x) \<noteq> 0` by auto

 	finally show ?thesis

 	  unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] float.cases Float_floor Let_def if_P[OF True] float_power .

       qed

     qed

     ultimately show ?thesis by auto

   qed

 qed

 lemma exp_boundaries: "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"

 proof -

   show ?thesis

   proof (cases "0 < x")

     case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto 

     from exp_boundaries'[OF this] show ?thesis .

   next

     case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto

     have "real (lb_exp prec x) \<le> exp (real x)"

     proof -

       from exp_boundaries'[OF `-x \<le> 0`]

       have ub_exp: "exp (- real x) \<le> real (ub_exp prec (-x))" unfolding atLeastAtMost_iff real_of_float_minus by auto

       have "real (float_divl prec 1 (ub_exp prec (-x))) \<le> 1 / real (ub_exp prec (-x))" using float_divl[where x=1] by auto

       also have "\<dots> \<le> exp (real x)"

 	using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]

 	unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto

       finally show ?thesis unfolding lb_exp.simps if_P[OF True] .

     qed

     moreover

     have "exp (real x) \<le> real (ub_exp prec x)"

     proof -

       have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto

       from exp_boundaries'[OF `-x \<le> 0`]

       have lb_exp: "real (lb_exp prec (-x)) \<le> exp (- real x)" unfolding atLeastAtMost_iff real_of_float_minus by auto

       have "exp (real x) \<le> real (1 :: float) / real (lb_exp prec (-x))"

 	using lb_exp[unfolded inverse_le_iff_le[OF exp_gt_zero lb_exp_pos[OF `\<not> 0 < -x`, unfolded less_float_def real_of_float_0], 

 	                                        symmetric]]

 	unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide real_of_float_1 by auto

       also have "\<dots> \<le> real (float_divr prec 1 (lb_exp prec (-x)))" using float_divr .

       finally show ?thesis unfolding ub_exp.simps if_P[OF True] .

     qed

     ultimately show ?thesis by auto

   qed

 qed

 lemma bnds_exp: "\<forall> x lx ux. (l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> exp x \<and> exp x \<le> real u"

 proof (rule allI, rule allI, rule allI, rule impI)

   fix x lx ux

   assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {real lx .. real ux}"

   hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto

   { from exp_boundaries[of lx prec, unfolded l]

     have "real l \<le> exp (real lx)" by (auto simp del: lb_exp.simps)

     also have "\<dots> \<le> exp x" using x by auto

     finally have "real l \<le> exp x" .

   } moreover

   { have "exp x \<le> exp (real ux)" using x by auto

     also have "\<dots> \<le> real u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)

     finally have "exp x \<le> real u" .

   } ultimately show "real l \<le> exp x \<and> exp x \<le> real u" ..

 qed

 section "Logarithm"

 subsection "Compute the logarithm series"

 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" 

 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

 "ub_ln_horner prec 0 i x       = 0" |

 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |

 "lb_ln_horner prec 0 i x       = 0" |

 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"

 lemma ln_bounds:

   assumes "0 \<le> x" and "x < 1"

   shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x ^ (Suc i)) \<le> ln (x + 1)" (is "?lb")

   and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub")

 proof -

   let "?a n" = "(1/real (n +1)) * x ^ (Suc n)"

   have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"

     using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto

   have "norm x < 1" using assms by auto

   have "?a ----> 0" unfolding Suc_plus1[symmetric] inverse_eq_divide[symmetric] 

     using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto

   { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }

   { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]

     proof (rule mult_mono)

       show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)

       have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 real_mult_assoc[symmetric] 

 	by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)

       thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto

     qed auto }

   from summable_Leibniz'(2,4)[OF `?a ----> 0` `\<And>n. 0 \<le> ?a n`, OF `\<And>n. ?a (Suc n) \<le> ?a n`, unfolded ln_eq]

   show "?lb" and "?ub" by auto

 qed

 lemma ln_float_bounds: 

   assumes "0 \<le> real x" and "real x < 1"

   shows "real (x * lb_ln_horner prec (get_even n) 1 x) \<le> ln (real x + 1)" (is "?lb \<le> ?ln")

   and "ln (real x + 1) \<le> real (x * ub_ln_horner prec (get_odd n) 1 x)" (is "?ln \<le> ?ub")

 proof -

   obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..

   obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..

   let "?s n" = "-1^n * (1 / real (1 + n)) * (real x)^(Suc n)"

   have "?lb \<le> setsum ?s {0 ..< 2 * ev}" unfolding power_Suc2 real_mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "real x"] ev

     using horner_bounds(1)[where G="\<lambda> i k. Suc k" and F="\<lambda>x. x" and f="\<lambda>x. x" and lb="\<lambda>n i k x. lb_ln_horner prec n k x" and ub="\<lambda>n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*ev",

       OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`

     by (rule mult_right_mono)

   also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> real x` `real x < 1`] by auto

   finally show "?lb \<le> ?ln" . 

   have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> real x` `real x < 1`] by auto

   also have "\<dots> \<le> ?ub" unfolding power_Suc2 real_mult_assoc[symmetric] real_of_float_mult setsum_left_distrib[symmetric] unfolding real_mult_commute[of "real x"] od

     using horner_bounds(2)[where G="\<lambda> i k. Suc k" and F="\<lambda>x. x" and f="\<lambda>x. x" and lb="\<lambda>n i k x. lb_ln_horner prec n k x" and ub="\<lambda>n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*od+1",

       OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`

     by (rule mult_right_mono)

   finally show "?ln \<le> ?ub" . 

 qed

 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"

 proof -

   have "x \<noteq> 0" using assms by auto

   have "x + y = x * (1 + y / x)" unfolding right_distrib times_divide_eq_right nonzero_mult_divide_cancel_left[OF `x \<noteq> 0`] by auto

   moreover 

   have "0 < y / x" using assms divide_pos_pos by auto

   hence "0 < 1 + y / x" by auto

   ultimately show ?thesis using ln_mult assms by auto

 qed

 subsection "Compute the logarithm of 2"

 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3 

                                         in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) + 

                                            (third * ub_ln_horner prec (get_odd prec) 1 third))"

 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3 

                                         in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) + 

                                            (third * lb_ln_horner prec (get_even prec) 1 third))"

 lemma ub_ln2: "ln 2 \<le> real (ub_ln2 prec)" (is "?ub_ln2")

   and lb_ln2: "real (lb_ln2 prec) \<le> ln 2" (is "?lb_ln2")

 proof -

   let ?uthird = "rapprox_rat (max prec 1) 1 3"

   let ?lthird = "lapprox_rat prec 1 3"

   have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"

     using ln_add[of "3 / 2" "1 / 2"] by auto

   have lb3: "real ?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto

   hence lb3_ub: "real ?lthird < 1" by auto

   have lb3_lb: "0 \<le> real ?lthird" using lapprox_rat_bottom[of 1 3] by auto

   have ub3: "1 / 3 \<le> real ?uthird" using rapprox_rat[of 1 3] by auto

   hence ub3_lb: "0 \<le> real ?uthird" by auto

   have lb2: "0 \<le> real (Float 1 -1)" and ub2: "real (Float 1 -1) < 1" unfolding Float_num by auto

   have "0 \<le> (1::int)" and "0 < (3::int)" by auto

   have ub3_ub: "real ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]

     by (rule rapprox_posrat_less1, auto)

   have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto

   have uthird_gt0: "0 < real ?uthird + 1" using ub3_lb by auto

   have lthird_gt0: "0 < real ?lthird + 1" using lb3_lb by auto

   show ?ub_ln2 unfolding ub_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]

   proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])

     have "ln (1 / 3 + 1) \<le> ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto

     also have "\<dots> \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)"

       using ln_float_bounds(2)[OF ub3_lb ub3_ub] .

     finally show "ln (1 / 3 + 1) \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .

   qed

   show ?lb_ln2 unfolding lb_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]

   proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])

     have "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (real ?lthird + 1)"

       using ln_float_bounds(1)[OF lb3_lb lb3_ub] .

     also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto

     finally show "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .

   qed

 qed

 subsection "Compute the logarithm in the entire domain"

 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where

 "ub_ln prec x = (if x \<le> 0          then None

             else if x < 1          then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))

             else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in

                  if x \<le> Float 3 -1 then Some (horner (x - 1))

             else if x < Float 1 1  then Some (horner (Float 1 -1) + horner (x * rapprox_rat prec 2 3 - 1))

                                    else let l = bitlen (mantissa x) - 1 in

                                         Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))" |

 "lb_ln prec x = (if x \<le> 0          then None

             else if x < 1          then Some (- the (ub_ln prec (float_divr prec 1 x)))

             else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in

                  if x \<le> Float 3 -1 then Some (horner (x - 1))

             else if x < Float 1 1  then Some (horner (Float 1 -1) +

                                               horner (max (x * lapprox_rat prec 2 3 - 1) 0))

                                    else let l = bitlen (mantissa x) - 1 in

                                         Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))"

 by pat_completeness auto

 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)

   fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"

   hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto

   from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]

   show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto

 next

   fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"

   hence "0 < x" unfolding less_float_def le_float_def by auto

   from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]

   show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto

 qed

 lemma ln_shifted_float: assumes "0 < m" shows "ln (real (Float m e)) = ln 2 * real (e + (bitlen m - 1)) + ln (real (Float m (- (bitlen m - 1))))"

 proof -

   let ?B = "2^nat (bitlen m - 1)"

   have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto

   hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto

   show ?thesis

   proof (cases "0 \<le> e")

     case True

     show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]

       unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]

       unfolding real_of_float_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]

       ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto

   next

     case False hence "0 < -e" by auto

     hence pow_gt0: "(0::real) < 2^nat (-e)" by auto

     hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto

     show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]

       unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]

       unfolding real_of_float_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]

       ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto

   qed

 qed

 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"

   shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"

   (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")

 proof (cases "x < Float 1 1")

   case True

   hence "real (x - 1) < 1" and "real x < 2" unfolding less_float_def Float_num by auto

   have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto

   hence "0 \<le> real (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto

   have [simp]: "real (Float 3 -1) = 3 / 2" by (simp add: real_of_float_def pow2_def)

   show ?thesis

   proof (cases "x \<le> Float 3 -1")

     case True

     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def

       using ln_float_bounds[OF `0 \<le> real (x - 1)` `real (x - 1) < 1`, of prec] `\<not> x \<le> 0` `\<not> x < 1` True

       by auto

   next

     case False hence *: "3 / 2 < real x" by (auto simp add: le_float_def)

     with ln_add[of "3 / 2" "real x - 3 / 2"]

     have add: "ln (real x) = ln (3 / 2) + ln (real x * 2 / 3)"

       by (auto simp add: algebra_simps diff_divide_distrib)

     let "?ub_horner x" = "x * ub_ln_horner prec (get_odd prec) 1 x"

     let "?lb_horner x" = "x * lb_ln_horner prec (get_even prec) 1 x"

     { have up: "real (rapprox_rat prec 2 3) \<le> 1"

 	by (rule rapprox_rat_le1) simp_all

       have low: "2 / 3 \<le> real (rapprox_rat prec 2 3)"

 	by (rule order_trans[OF _ rapprox_rat]) simp

       from mult_less_le_imp_less[OF * low] *

       have pos: "0 < real (x * rapprox_rat prec 2 3 - 1)" by auto

       have "ln (real x * 2/3)

 	\<le> ln (real (x * rapprox_rat prec 2 3 - 1) + 1)"

       proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])

 	show "real x * 2 / 3 \<le> real (x * rapprox_rat prec 2 3 - 1) + 1"

 	  using * low by auto

 	show "0 < real x * 2 / 3" using * by simp

 	show "0 < real (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto

       qed

       also have "\<dots> \<le> real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"

       proof (rule ln_float_bounds(2))

 	from mult_less_le_imp_less[OF `real x < 2` up] low *

 	show "real (x * rapprox_rat prec 2 3 - 1) < 1" by auto

 	show "0 \<le> real (x * rapprox_rat prec 2 3 - 1)" using pos by auto

       qed

       finally have "ln (real x)

 	\<le> real (?ub_horner (Float 1 -1))

 	  + real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"

 	using ln_float_bounds(2)[of "Float 1 -1" prec prec] add by auto }

     moreover

     { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0"

       have up: "real (lapprox_rat prec 2 3) \<le> 2/3"

 	by (rule order_trans[OF lapprox_rat], simp)

       have low: "0 \<le> real (lapprox_rat prec 2 3)"

 	using lapprox_rat_bottom[of 2 3 prec] by simp

       have "real (?lb_horner ?max)

 	\<le> ln (real ?max + 1)"

       proof (rule ln_float_bounds(1))

 	from mult_less_le_imp_less[OF `real x < 2` up] * low

 	show "real ?max < 1" by (cases "real (lapprox_rat prec 2 3) = 0",

 	  auto simp add: real_of_float_max)

 	show "0 \<le> real ?max" by (auto simp add: real_of_float_max)

       qed

       also have "\<dots> \<le> ln (real x * 2/3)"

       proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])

 	show "0 < real ?max + 1" by (auto simp add: real_of_float_max)

 	show "0 < real x * 2/3" using * by auto

 	show "real ?max + 1 \<le> real x * 2/3" using * up

 	  by (cases "0 < real x * real (lapprox_posrat prec 2 3) - 1",

 	      auto simp add: real_of_float_max max_def)

       qed

       finally have "real (?lb_horner (Float 1 -1)) + real (?lb_horner ?max)

 	\<le> ln (real x)"

 	using ln_float_bounds(1)[of "Float 1 -1" prec prec] add by auto }

     ultimately

     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def

       using `\<not> x \<le> 0` `\<not> x < 1` True False by auto

   qed

 next

   case False

   hence "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" "\<not> x \<le> Float 3 -1"

     using `1 \<le> x` unfolding less_float_def le_float_def real_of_float_simp pow2_def

     by auto

   show ?thesis

   proof (cases x)

     case (Float m e)

     let ?s = "Float (e + (bitlen m - 1)) 0"

     let ?x = "Float m (- (bitlen m - 1))"

     have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto

     {

       have "real (lb_ln2 prec * ?s) \<le> ln 2 * real (e + (bitlen m - 1))" (is "?lb2 \<le> _")

 	unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right

 	using lb_ln2[of prec]

       proof (rule mult_right_mono)

 	have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto

 	from float_gt1_scale[OF this]

 	show "0 \<le> real (e + (bitlen m - 1))" by auto

       qed

       moreover

       from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]

       have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto

       from ln_float_bounds(1)[OF this]

       have "real ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln (real ?x)" (is "?lb_horner \<le> _") by auto

       ultimately have "?lb2 + ?lb_horner \<le> ln (real x)"

 	unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto

     }

     moreover

     {

       from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]

       have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto

       from ln_float_bounds(2)[OF this]

       have "ln (real ?x) \<le> real ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> ?ub_horner") by auto

       moreover

       have "ln 2 * real (e + (bitlen m - 1)) \<le> real (ub_ln2 prec * ?s)" (is "_ \<le> ?ub2")

 	unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right

 	using ub_ln2[of prec]

       proof (rule mult_right_mono)

 	have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto

 	from float_gt1_scale[OF this]

 	show "0 \<le> real (e + (bitlen m - 1))" by auto

       qed

       ultimately have "ln (real x) \<le> ?ub2 + ?ub_horner"

 	unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto

     }

     ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps

       unfolding if_not_P[OF `\<not> x \<le> 0`] if_not_P[OF `\<not> x < 1`] if_not_P[OF False] if_not_P[OF `\<not> x \<le> Float 3 -1`] Let_def

       unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] real_of_float_add

       by auto

   qed

 qed

 lemma ub_ln_lb_ln_bounds: assumes "0 < x"

   shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"

   (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")

 proof (cases "x < 1")

   case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto

   show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .

 next

   case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto

   have "0 < real x" and "real x \<noteq> 0" using `0 < x` unfolding less_float_def by auto

   hence A: "0 < 1 / real x" by auto

   {

     let ?divl = "float_divl (max prec 1) 1 x"

     have A': "1 \<le> ?divl" using float_divl_pos_less1_bound[OF `0 < x` `x < 1`] unfolding le_float_def less_float_def by auto

     hence B: "0 < real ?divl" unfolding le_float_def by auto

     have "ln (real ?divl) \<le> ln (1 / real x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto

     hence "ln (real x) \<le> - ln (real ?divl)" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto

     from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]

     have "?ln \<le> real (- the (lb_ln prec ?divl))" unfolding real_of_float_minus by (rule order_trans)

   } moreover

   {

     let ?divr = "float_divr prec 1 x"

     have A': "1 \<le> ?divr" using float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`] unfolding le_float_def less_float_def by auto

     hence B: "0 < real ?divr" unfolding le_float_def by auto

     have "ln (1 / real x) \<le> ln (real ?divr)" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto

     hence "- ln (real ?divr) \<le> ln (real x)" unfolding nonzero_inverse_eq_divide[OF `real x \<noteq> 0`, symmetric] ln_inverse[OF `0 < real x`] by auto

     from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this

     have "real (- the (ub_ln prec ?divr)) \<le> ?ln" unfolding real_of_float_minus by (rule order_trans)

   }

   ultimately show ?thesis unfolding lb_ln.simps[where x=x]  ub_ln.simps[where x=x]

     unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto

 qed

 lemma lb_ln: assumes "Some y = lb_ln prec x"

   shows "real y \<le> ln (real x)" and "0 < real x"

 proof -

   have "0 < x"

   proof (rule ccontr)

     assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto

     thus False using assms by auto

   qed

   thus "0 < real x" unfolding less_float_def by auto

   have "real (the (lb_ln prec x)) \<le> ln (real x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..

   thus "real y \<le> ln (real x)" unfolding assms[symmetric] by auto

 qed

 lemma ub_ln: assumes "Some y = ub_ln prec x"

   shows "ln (real x) \<le> real y" and "0 < real x"

 proof -

   have "0 < x"

   proof (rule ccontr)

     assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto

     thus False using assms by auto

   qed

   thus "0 < real x" unfolding less_float_def by auto

   have "ln (real x) \<le> real (the (ub_ln prec x))" using ub_ln_lb_ln_bounds[OF `0 < x`] ..

   thus "ln (real x) \<le> real y" unfolding assms[symmetric] by auto

 qed

 lemma bnds_ln: "\<forall> x lx ux. (Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {real lx .. real ux} \<longrightarrow> real l \<le> ln x \<and> ln x \<le> real u"

 proof (rule allI, rule allI, rule allI, rule impI)

   fix x lx ux

   assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {real lx .. real ux}"

   hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \<in> {real lx .. real ux}" by auto

   have "ln (real ux) \<le> real u" and "0 < real ux" using ub_ln u by auto

   have "real l \<le> ln (real lx)" and "0 < real lx" and "0 < x" using lb_ln[OF l] x by auto

   from ln_le_cancel_iff[OF `0 < real lx` `0 < x`] `real l \<le> ln (real lx)`

   have "real l \<le> ln x" using x unfolding atLeastAtMost_iff by auto

   moreover

   from ln_le_cancel_iff[OF `0 < x` `0 < real ux`] `ln (real ux) \<le> real u`

   have "ln x \<le> real u" using x unfolding atLeastAtMost_iff by auto

   ultimately show "real l \<le> ln x \<and> ln x \<le> real u" ..

 qed

 section "Implement floatarith"

 subsection "Define syntax and semantics"

 datatype floatarith

   = Add floatarith floatarith

   | Minus floatarith

   | Mult floatarith floatarith

   | Inverse floatarith

   | Cos floatarith

   | Arctan floatarith

   | Abs floatarith

   | Max floatarith floatarith

   | Min floatarith floatarith

   | Pi

   | Sqrt floatarith

   | Exp floatarith

   | Ln floatarith

   | Power floatarith nat

   | Atom nat

   | Num float

 fun interpret_floatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real"

 where

 "interpret_floatarith (Add a b) vs   = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |

 "interpret_floatarith (Minus a) vs    = - (interpret_floatarith a vs)" |

 "interpret_floatarith (Mult a b) vs   = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |

 "interpret_floatarith (Inverse a) vs  = inverse (interpret_floatarith a vs)" |

 "interpret_floatarith (Cos a) vs      = cos (interpret_floatarith a vs)" |

 "interpret_floatarith (Arctan a) vs   = arctan (interpret_floatarith a vs)" |

 "interpret_floatarith (Min a b) vs    = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |

 "interpret_floatarith (Max a b) vs    = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |

 "interpret_floatarith (Abs a) vs      = abs (interpret_floatarith a vs)" |

 "interpret_floatarith Pi vs           = pi" |

 "interpret_floatarith (Sqrt a) vs     = sqrt (interpret_floatarith a vs)" |

 "interpret_floatarith (Exp a) vs      = exp (interpret_floatarith a vs)" |

 "interpret_floatarith (Ln a) vs       = ln (interpret_floatarith a vs)" |

 "interpret_floatarith (Power a n) vs  = (interpret_floatarith a vs)^n" |

 "interpret_floatarith (Num f) vs      = real f" |

 "interpret_floatarith (Atom n) vs     = vs ! n"

 subsection "Implement approximation function"

 fun lift_bin' :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where

 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |

 "lift_bin' a b f = None"

 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where

 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)

                                              | t \<Rightarrow> None)" |

 "lift_un b f = None"

 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where

 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |

 "lift_un' b f = None"

 fun bounded_by :: "real list \<Rightarrow> (float * float) list \<Rightarrow> bool " where

 bounded_by_Cons: "bounded_by (v#vs) ((l, u)#bs) = ((real l \<le> v \<and> v \<le> real u) \<and> bounded_by vs bs)" |

 bounded_by_Nil: "bounded_by [] [] = True" |

 "bounded_by _ _ = False"

 lemma bounded_by: assumes "bounded_by vs bs" and "i < length bs"

   shows "real (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> real (snd (bs ! i))"

   using `bounded_by vs bs` and `i < length bs`

 proof (induct arbitrary: i rule: bounded_by.induct)

   fix v :: real and vs :: "real list" and l u :: float and bs :: "(float * float) list" and i :: nat

   assume hyp: "\<And>i. \<lbrakk>bounded_by vs bs; i < length bs\<rbrakk> \<Longrightarrow> real (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> real (snd (bs ! i))"

   assume bounded: "bounded_by (v # vs) ((l, u) # bs)" and length: "i < length ((l, u) # bs)"

   show "real (fst (((l, u) # bs) ! i)) \<le> (v # vs) ! i \<and> (v # vs) ! i \<le> real (snd (((l, u) # bs) ! i))"

   proof (cases i)

     case 0

     show ?thesis using bounded unfolding 0 nth_Cons_0 fst_conv snd_conv bounded_by.simps ..

   next

     case (Suc i) with length have "i < length bs" by auto

     show ?thesis unfolding Suc nth_Cons_Suc bounded_by.simps

       using hyp[OF bounded[unfolded bounded_by.simps, THEN conjunct2] `i < length bs`] .

   qed

 qed auto

 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) list \<Rightarrow> (float * float) option" where

 "approx' prec a bs          = (case (approx prec a bs) of Some (l, u) \<Rightarrow> Some (round_down prec l, round_up prec u) | None \<Rightarrow> None)" |

 "approx prec (Add a b) bs  = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (l1 + l2, u1 + u2))" | 

 "approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |

 "approx prec (Mult a b) bs  = lift_bin' (approx' prec a bs) (approx' prec b bs)

                                     (\<lambda> a1 a2 b1 b2. (float_nprt a1 * float_pprt b2 + float_nprt a2 * float_nprt b2 + float_pprt a1 * float_pprt b1 + float_pprt a2 * float_nprt b1, 

                                                      float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |

 "approx prec (Inverse a) bs = lift_un (approx' prec a bs) (\<lambda> l u. if (0 < l \<or> u < 0) then (Some (float_divl prec 1 u), Some (float_divr prec 1 l)) else (None, None))" |

 "approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |

 "approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |

 "approx prec (Min a b) bs   = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (min l1 l2, min u1 u2))" |

 "approx prec (Max a b) bs   = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (max l1 l2, max u1 u2))" |

 "approx prec (Abs a) bs     = lift_un' (approx' prec a bs) (\<lambda>l u. (if l < 0 \<and> 0 < u then 0 else min \<bar>l\<bar> \<bar>u\<bar>, max \<bar>l\<bar> \<bar>u\<bar>))" |

 "approx prec (Arctan a) bs  = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |

 "approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |

 "approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |

 "approx prec (Ln a) bs      = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |

 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |

 "approx prec (Num f) bs     = Some (f, f)" |

 "approx prec (Atom i) bs    = (if i < length bs then Some (bs ! i) else None)"

 lemma lift_bin'_ex:

   assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"

   shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"

 proof (cases a)

   case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..

   thus ?thesis using lift_bin'_Some by auto

 next

   case (Some a')

   show ?thesis

   proof (cases b)

     case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..

     thus ?thesis using lift_bin'_Some by auto

   next

     case (Some b')

     obtain la ua where a': "a' = (la, ua)" by (cases a', auto)

     obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)

     thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto

   qed

 qed

 lemma lift_bin'_f:

   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"

   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a" and Pb: "\<And>l u. Some (l, u) = g b \<Longrightarrow> P l u b"

   shows "\<exists> l1 u1 l2 u2. P l1 u1 a \<and> P l2 u2 b \<and> l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"

 proof -

   obtain l1 u1 l2 u2

     where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto

   have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto 

   have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto

   thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto 

 qed

 lemma approx_approx':

   assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"

   and approx': "Some (l, u) = approx' prec a vs"

   shows "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"

 proof -

   obtain l' u' where S: "Some (l', u') = approx prec a vs"

     using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)

   have l': "l = round_down prec l'" and u': "u = round_up prec u'"

     using approx' unfolding approx'.simps S[symmetric] by auto

   show ?thesis unfolding l' u' 

     using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]

     using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto

 qed

 lemma lift_bin':

   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"

   and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")

   and Pb: "\<And>l u. Some (l, u) = approx prec b bs \<Longrightarrow> real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u"

   shows "\<exists> l1 u1 l2 u2. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 

                         (real l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u2) \<and> 

                         l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"

 proof -

   { fix l u assume "Some (l, u) = approx' prec a bs"

     with approx_approx'[of prec a bs, OF _ this] Pa

     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this

   { fix l u assume "Some (l, u) = approx' prec b bs"

     with approx_approx'[of prec b bs, OF _ this] Pb

     have "real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u" by auto } note Pb = this

   from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]

   show ?thesis by auto

 qed

 lemma lift_un'_ex:

   assumes lift_un'_Some: "Some (l, u) = lift_un' a f"

   shows "\<exists> l u. Some (l, u) = a"

 proof (cases a)

   case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..

   thus ?thesis using lift_un'_Some by auto

 next

   case (Some a')

   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)

   thus ?thesis unfolding `a = Some a'` a' by auto

 qed

 lemma lift_un'_f:

   assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"

   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"

   shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"

 proof -

   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto

   have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto

   have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto

   thus ?thesis using Pa[OF Sa] by auto

 qed

 lemma lift_un':

   assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"

   and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")

   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 

                         l = fst (f l1 u1) \<and> u = snd (f l1 u1)"

 proof -

   { fix l u assume "Some (l, u) = approx' prec a bs"

     with approx_approx'[of prec a bs, OF _ this] Pa

     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this

   from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]

   show ?thesis by auto

 qed

 lemma lift_un'_bnds:

   assumes bnds: "\<forall> x lx ux. (l, u) = f lx ux \<and> x \<in> { real lx .. real ux } \<longrightarrow> real l \<le> f' x \<and> f' x \<le> real u"

   and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"

   and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"

   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"

 proof -

   from lift_un'[OF lift_un'_Some Pa]

   obtain l1 u1 where "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" and "l = fst (f l1 u1)" and "u = snd (f l1 u1)" by blast

   hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto

   thus ?thesis using bnds by auto

 qed

 lemma lift_un_ex:

   assumes lift_un_Some: "Some (l, u) = lift_un a f"

   shows "\<exists> l u. Some (l, u) = a"

 proof (cases a)

   case None hence "None = lift_un a f" unfolding None lift_un.simps ..

   thus ?thesis using lift_un_Some by auto

 next

   case (Some a')

   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)

   thus ?thesis unfolding `a = Some a'` a' by auto

 qed

 lemma lift_un_f:

   assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"

   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"

   shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"

 proof -

   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto

   have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"

   proof (rule ccontr)

     assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"

     hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto

     hence "lift_un (g a) f = None" 

     proof (cases "fst (f l1 u1) = None")

       case True

       then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)

       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto

     next

       case False hence "snd (f l1 u1) = None" using or by auto

       with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)

       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto

     qed

     thus False using lift_un_Some by auto

   qed

   then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)

   from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]

   have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto

   thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto

 qed

 lemma lift_un:

   assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"

   and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")

   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 

                   Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"

 proof -

   { fix l u assume "Some (l, u) = approx' prec a bs"

     with approx_approx'[of prec a bs, OF _ this] Pa

     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this

   from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]

   show ?thesis by auto

 qed

 lemma lift_un_bnds:

   assumes bnds: "\<forall> x lx ux. (Some l, Some u) = f lx ux \<and> x \<in> { real lx .. real ux } \<longrightarrow> real l \<le> f' x \<and> f' x \<le> real u"

   and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"

   and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow> real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"

   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"

 proof -

   from lift_un[OF lift_un_Some Pa]

   obtain l1 u1 where "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" and "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" by blast

   hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto

   thus ?thesis using bnds by auto

 qed

 lemma approx:

   assumes "bounded_by xs vs"

   and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")

   shows "real l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> real u" (is "?P l u arith")

   using `Some (l, u) = approx prec arith vs` 

 proof (induct arith arbitrary: l u x)

   case (Add a b)

   from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps

   obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"

     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"

     "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast

   thus ?case unfolding interpret_floatarith.simps by auto

 next

   case (Minus a)

   from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps

   obtain l1 u1 where "l = -u1" and "u = -l1"

     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" unfolding fst_conv snd_conv by blast

   thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto

 next

   case (Mult a b)

   from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps

   obtain l1 u1 l2 u2 

     where l: "l = float_nprt l1 * float_pprt u2 + float_nprt u1 * float_nprt u2 + float_pprt l1 * float_pprt l2 + float_pprt u1 * float_nprt l2"

     and u: "u = float_pprt u1 * float_pprt u2 + float_pprt l1 * float_nprt u2 + float_nprt u1 * float_pprt l2 + float_nprt l1 * float_nprt l2"

     and "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"

     and "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast

   thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt 

     using mult_le_prts mult_ge_prts by auto

 next

   case (Inverse a)

   from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps

   obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)" 

     and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"

     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1" by blast

   have either: "0 < l1 \<or> u1 < 0" proof (rule ccontr) assume P: "\<not> (0 < l1 \<or> u1 < 0)" show False using l' unfolding if_not_P[OF P] by auto qed

   moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto

   ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto

   have inv: "inverse (real u1) \<le> inverse (interpret_floatarith a xs)

            \<and> inverse (interpret_floatarith a xs) \<le> inverse (real l1)"

   proof (cases "0 < l1")

     case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs" 

       unfolding less_float_def using l1_le_u1 l1 by auto

     show ?thesis

       unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]

 	inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]

       using l1 u1 by auto

   next

     case False hence "u1 < 0" using either by blast

     hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0" 

       unfolding less_float_def using l1_le_u1 u1 by auto

     show ?thesis

       unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]

 	inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]

       using l1 u1 by auto

   qed

   from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)

   hence "real l \<le> inverse (real u1)" unfolding nonzero_inverse_eq_divide[OF `real u1 \<noteq> 0`] using float_divl[of prec 1 u1] by auto

   also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto

   finally have "real l \<le> inverse (interpret_floatarith a xs)" .

   moreover

   from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)

   hence "inverse (real l1) \<le> real u" unfolding nonzero_inverse_eq_divide[OF `real l1 \<noteq> 0`] using float_divr[of 1 l1 prec] by auto

   hence "inverse (interpret_floatarith a xs) \<le> real u" by (rule order_trans[OF inv[THEN conjunct2]])

   ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto

 next

   case (Abs x)

   from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps

   obtain l1 u1 where l': "l = (if l1 < 0 \<and> 0 < u1 then 0 else min \<bar>l1\<bar> \<bar>u1\<bar>)" and u': "u = max \<bar>l1\<bar> \<bar>u1\<bar>"

     and l1: "real l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> real u1" by blast

   thus ?case unfolding l' u' by (cases "l1 < 0 \<and> 0 < u1", auto simp add: real_of_float_min real_of_float_max real_of_float_abs less_float_def)

 next

   case (Min a b)

   from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps

   obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"

     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"

     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast

   thus ?case unfolding l' u' by (auto simp add: real_of_float_min)

 next

   case (Max a b)

   from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps

   obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"

     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"

     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast

   thus ?case unfolding l' u' by (auto simp add: real_of_float_max)

 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto

 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto

 next case Pi with pi_boundaries show ?case by auto

 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto

 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto

 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto

 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto

 next case (Num f) thus ?case by auto

 next

   case (Atom n) 

   show ?case

   proof (cases "n < length vs")

     case True

     with Atom have "vs ! n = (l, u)" by auto

     thus ?thesis using bounded_by[OF assms(1) True] by auto

   next

     case False thus ?thesis using Atom by auto

   qed

 qed

 datatype inequality = Less floatarith floatarith 

                     | LessEqual floatarith floatarith 

 fun interpret_inequality :: "inequality \<Rightarrow> real list \<Rightarrow> bool" where 

 "interpret_inequality (Less a b) vs                   = (interpret_floatarith a vs < interpret_floatarith b vs)" |

 "interpret_inequality (LessEqual a b) vs              = (interpret_floatarith a vs \<le> interpret_floatarith b vs)"

 fun approx_inequality :: "nat \<Rightarrow> inequality \<Rightarrow> (float * float) list \<Rightarrow> bool" where 

 "approx_inequality prec (Less a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u < l' | _ \<Rightarrow> False)" |

 "approx_inequality prec (LessEqual a b) bs = (case (approx prec a bs, approx prec b bs) of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l' | _ \<Rightarrow> False)"

 lemma approx_inequality: fixes m :: nat assumes "bounded_by vs bs" and "approx_inequality prec eq bs"

   shows "interpret_inequality eq vs"

 proof (cases eq)

   case (Less a b)

   show ?thesis

   proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and> 

                              approx prec b bs = Some (l', u')")

     case True

     then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"

       and b_approx: "approx prec b bs = Some (l', u') " by auto

     with `approx_inequality prec eq bs` have "real u < real l'"

       unfolding Less approx_inequality.simps less_float_def by auto

     moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`

     have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"

       using approx by auto

     ultimately show ?thesis unfolding interpret_inequality.simps Less by auto

   next

     case False

     hence "approx prec a bs = None \<or> approx prec b bs = None"

       unfolding not_Some_eq[symmetric] by auto

     hence "\<not> approx_inequality prec eq bs" unfolding Less approx_inequality.simps 

       by (cases "approx prec a bs = None", auto)

     thus ?thesis using assms by auto

   qed

 next

   case (LessEqual a b)

   show ?thesis

   proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and> 

                              approx prec b bs = Some (l', u')")

     case True

     then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"

       and b_approx: "approx prec b bs = Some (l', u') " by auto

     with `approx_inequality prec eq bs` have "real u \<le> real l'"

       unfolding LessEqual approx_inequality.simps le_float_def by auto

     moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`

     have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"

       using approx by auto

     ultimately show ?thesis unfolding interpret_inequality.simps LessEqual by auto

   next

     case False

     hence "approx prec a bs = None \<or> approx prec b bs = None"

       unfolding not_Some_eq[symmetric] by auto

     hence "\<not> approx_inequality prec eq bs" unfolding LessEqual approx_inequality.simps 

       by (cases "approx prec a bs = None", auto)

     thus ?thesis using assms by auto

   qed

 qed

 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"

   unfolding real_divide_def interpret_floatarith.simps ..

 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"

   unfolding real_diff_def interpret_floatarith.simps ..

 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs = 

   sin (interpret_floatarith a vs)"

   unfolding sin_cos_eq interpret_floatarith.simps

             interpret_floatarith_divide interpret_floatarith_diff real_diff_def

   by auto

 lemma interpret_floatarith_tan:

   "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =

    tan (interpret_floatarith a vs)"

   unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def real_divide_def

   by auto

 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"

   unfolding powr_def interpret_floatarith.simps ..

 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"

   unfolding log_def interpret_floatarith.simps real_divide_def ..

 lemma interpret_floatarith_num:

   shows "interpret_floatarith (Num (Float 0 0)) vs = 0"

   and "interpret_floatarith (Num (Float 1 0)) vs = 1"

   and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto

 subsection {* Implement proof method \texttt{approximation} *}

 lemma bounded_divl: assumes "real a / real b \<le> x" shows "real (float_divl p a b) \<le> x" by (rule order_trans[OF _ assms], rule float_divl)

 lemma bounded_divr: assumes "x \<le> real a / real b" shows "x \<le> real (float_divr p a b)" by (rule order_trans[OF assms _], rule float_divr)

 lemma bounded_num: shows "real (Float 5 1) = 10" and "real (Float 0 0) = 0" and "real (Float 1 0) = 1" and "real (Float (number_of n) 0) = (number_of n)"

                      and "0 * pow2 e = real (Float 0 e)" and "1 * pow2 e = real (Float 1 e)" and "number_of m * pow2 e = real (Float (number_of m) e)"

                      and "real (Float (number_of A) (int B)) = (number_of A) * 2^B"

                      and "real (Float 1 (int B)) = 2^B"

                      and "real (Float (number_of A) (- int B)) = (number_of A) / 2^B"

                      and "real (Float 1 (- int B)) = 1 / 2^B"

   by (auto simp add: real_of_float_simp pow2_def real_divide_def)

 lemmas bounded_by_equations = bounded_by_Cons bounded_by_Nil float_power bounded_divl bounded_divr bounded_num HOL.simp_thms

 lemmas interpret_inequality_equations = interpret_inequality.simps interpret_floatarith.simps interpret_floatarith_num

   interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log

   interpret_floatarith_sin

 ML {*

 structure Float_Arith =

 struct

 @{code_datatype float = Float}

 @{code_datatype floatarith = Add | Minus | Mult | Inverse | Cos | Arctan

                            | Abs | Max | Min | Pi | Sqrt | Exp | Ln | Power | Atom | Num }

 @{code_datatype inequality = Less | LessEqual }

 val approx_inequality = @{code approx_inequality}

 end

 *}

 code_reserved Eval Float_Arith

 code_type float (Eval "Float'_Arith.float")

 code_const Float (Eval "Float'_Arith.Float/ (_,/ _)")

 code_type floatarith (Eval "Float'_Arith.floatarith")

 code_const Add and Minus and Mult and Inverse and Cos and Arctan and Abs and Max and Min and

            Pi and Sqrt  and Exp and Ln and Power and Atom and Num

   (Eval "Float'_Arith.Add/ (_,/ _)" and "Float'_Arith.Minus" and "Float'_Arith.Mult/ (_,/ _)" and

         "Float'_Arith.Inverse" and "Float'_Arith.Cos" and

         "Float'_Arith.Arctan" and "Float'_Arith.Abs" and "Float'_Arith.Max/ (_,/ _)" and

         "Float'_Arith.Min/ (_,/ _)" and "Float'_Arith.Pi" and "Float'_Arith.Sqrt" and

         "Float'_Arith.Exp" and "Float'_Arith.Ln" and "Float'_Arith.Power/ (_,/ _)" and

         "Float'_Arith.Atom" and "Float'_Arith.Num")

 code_type inequality (Eval "Float'_Arith.inequality")

 code_const Less and LessEqual (Eval "Float'_Arith.Less/ (_,/ _)" and "Float'_Arith.LessEqual/ (_,/ _)")

 code_const approx_inequality (Eval "Float'_Arith.approx'_inequality")

 ML {*

   val ineq_equations = PureThy.get_thms @{theory} "interpret_inequality_equations";

   val bounded_by_equations = PureThy.get_thms @{theory} "bounded_by_equations";

   val bounded_by_simpset = (HOL_basic_ss addsimps bounded_by_equations)

   fun reify_ineq ctxt i = (fn st =>

     let

       val to = HOLogic.dest_Trueprop (Logic.strip_imp_concl (List.nth (prems_of st, i - 1)))

     in (Reflection.genreify_tac ctxt ineq_equations (SOME to) i) st

     end)

   fun rule_ineq ctxt prec i thm = let

     fun conv_num typ = HOLogic.dest_number #> snd #> HOLogic.mk_number typ

     val to_natc = conv_num @{typ "nat"} #> Thm.cterm_of (ProofContext.theory_of ctxt)

     val to_nat = conv_num @{typ "nat"}

     val to_int = conv_num @{typ "int"}

     fun int_to_float A = @{term "Float"} $ to_int A $ @{term "0 :: int"}

     val prec' = to_nat prec

     fun bot_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))

                    = @{term "Float"} $ to_int mantisse $ to_int exp

       | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))

       | bot_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)

       | bot_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "float_divl"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)

       | bot_float (Const (@{const_name "divide"}, _) $ A $ B)

                    = @{term "float_divl"} $ prec' $ int_to_float A $ int_to_float B

       | bot_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)

                    = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)

       | bot_float A = int_to_float A

     fun top_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))

                    = @{term "Float"} $ to_int mantisse $ to_int exp

       | top_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))

       | top_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)

       | top_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))

                    = @{term "float_divr"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)

       | top_float (Const (@{const_name "divide"}, _) $ A $ B)

                    = @{term "float_divr"} $ prec' $ int_to_float A $ int_to_float B

       | top_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)

                    = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)

       | top_float A = int_to_float A

     val goal' : term = List.nth (prems_of thm, i - 1)

     fun lift_bnd (t as (Const (@{const_name "op &"}, _) $

                         (Const (@{const_name "less_eq"}, _) $

                          bottom $ (Free (name, _))) $

                         (Const (@{const_name "less_eq"}, _) $ _ $ top)))

          = ((name, HOLogic.mk_prod (bot_float bottom, top_float top))

             handle TERM (txt, ts) => raise TERM ("Invalid bound number format: " ^

                                   (Syntax.string_of_term ctxt t), [t]))

       | lift_bnd t = raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^

                                  (Syntax.string_of_term ctxt t), [t])

     val bound_eqs = map (HOLogic.dest_Trueprop #> lift_bnd)  (Logic.strip_imp_prems goal')

     fun lift_var (Free (varname, _)) = (case AList.lookup (op =) bound_eqs varname of

                                           SOME bound => bound

                                         | NONE => raise TERM ("No bound equations found for " ^ varname, []))

       | lift_var t = raise TERM ("Can not convert expression " ^

                                  (Syntax.string_of_term ctxt t), [t])

     val _ $ vs = HOLogic.dest_Trueprop (Logic.strip_imp_concl goal')

     val bs = (HOLogic.dest_list #> map lift_var #> HOLogic.mk_list @{typ "float * float"}) vs

     val map = [(@{cpat "?prec::nat"}, to_natc prec),

                (@{cpat "?bs::(float * float) list"}, Thm.cterm_of (ProofContext.theory_of ctxt) bs)]

   in rtac (Thm.instantiate ([], map) @{thm "approx_inequality"}) i thm end

   val eval_tac = CSUBGOAL (fn (ct, i) => rtac (eval_oracle ct) i)

   fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)

                                THEN' rtac TrueI

 *}

 method_setup approximation = {*

   Args.term >>

   (fn prec => fn ctxt =>

     SIMPLE_METHOD' (fn i =>

      (DETERM (reify_ineq ctxt i)

       THEN rule_ineq ctxt prec i

       THEN Simplifier.asm_full_simp_tac bounded_by_simpset i

       THEN (TRY (filter_prems_tac (fn t => false) i))

       THEN (gen_eval_tac eval_oracle ctxt) i)))

 *} "real number approximation"

 end