src/HOL/Decision_Procs/Approximation.thy
author nipkow
Wed Jun 24 09:41:14 2009 +0200 (2009-06-24)
changeset 31790 05c92381363c
parent 31508 1ea1c70aba00
child 31809 06372924e86c
permissions -rw-r--r--
corrected and unified thm names
     1 (* Author:     Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 *)
     2 
     3 header {* Prove Real Valued Inequalities by Computation *}
     4 
     5 theory Approximation
     6 imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat
     7 begin
     8 
     9 section "Horner Scheme"
    10 
    11 subsection {* Define auxiliary helper @{text horner} function *}
    12 
    13 primrec horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where
    14 "horner F G 0 i k x       = 0" |
    15 "horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"
    16 
    17 lemma horner_schema': fixes x :: real  and a :: "nat \<Rightarrow> real"
    18   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)"
    19 proof -
    20   have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto
    21   show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]
    22     setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n  *a n * x^n"] by auto
    23 qed
    24 
    25 lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"
    26   assumes f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
    27   shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / real (f (j' + j))) * x ^ j)"
    28 proof (induct n arbitrary: i k j')
    29   case (Suc n)
    30 
    31   show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
    32     using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto
    33 qed auto
    34 
    35 lemma horner_bounds':
    36   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
    37   and lb_0: "\<And> i k x. lb 0 i k x = 0"
    38   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)"
    39   and ub_0: "\<And> i k x. ub 0 i k x = 0"
    40   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)"
    41   shows "real (lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') (real x) \<and> 
    42          horner F G n ((F ^^ j') s) (f j') (real x) \<le> real (ub n ((F ^^ j') s) (f j') x)"
    43   (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
    44 proof (induct n arbitrary: j')
    45   case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto
    46 next
    47   case (Suc n)
    48   have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps real_of_float_sub diff_def
    49   proof (rule add_mono)
    50     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
    51     from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> real x`
    52     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))"
    53       unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)
    54   qed
    55   moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps real_of_float_sub diff_def
    56   proof (rule add_mono)
    57     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
    58     from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> real x`
    59     show "- (real x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) (real x)) \<le> 
    60           - real (x * lb n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x)"
    61       unfolding real_of_float_mult neg_le_iff_le by (rule mult_left_mono)
    62   qed
    63   ultimately show ?case by blast
    64 qed
    65 
    66 subsection "Theorems for floating point functions implementing the horner scheme"
    67 
    68 text {*
    69 
    70 Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are
    71 all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.
    72 
    73 *}
    74 
    75 lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
    76   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
    77   and lb_0: "\<And> i k x. lb 0 i k x = 0"
    78   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)"
    79   and ub_0: "\<And> i k x. ub 0 i k x = 0"
    80   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)"
    81   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 
    82     "(\<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")
    83 proof -
    84   have "?lb  \<and> ?ub" 
    85     using horner_bounds'[where lb=lb, OF `0 \<le> real x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
    86     unfolding horner_schema[where f=f, OF f_Suc] .
    87   thus "?lb" and "?ub" by auto
    88 qed
    89 
    90 lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
    91   assumes "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
    92   and lb_0: "\<And> i k x. lb 0 i k x = 0"
    93   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)"
    94   and ub_0: "\<And> i k x. ub 0 i k x = 0"
    95   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)"
    96   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 
    97     "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (real x ^ j)) \<le> real (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
    98 proof -
    99   { fix x y z :: float have "x - y * z = x + - y * z"
   100       by (cases x, cases y, cases z, simp add: plus_float.simps minus_float_def uminus_float.simps times_float.simps algebra_simps)
   101   } note diff_mult_minus = this
   102 
   103   { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this
   104 
   105   have move_minus: "real (-x) = -1 * real x" by auto
   106 
   107   have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * real x ^ j) = 
   108     (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j)"
   109   proof (rule setsum_cong, simp)
   110     fix j assume "j \<in> {0 ..< n}"
   111     show "1 / real (f (j' + j)) * real x ^ j = -1 ^ j * (1 / real (f (j' + j))) * real (- x) ^ j"
   112       unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]
   113       unfolding real_mult_commute unfolding real_mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]
   114       by auto
   115   qed
   116 
   117   have "0 \<le> real (-x)" using assms by auto
   118   from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
   119     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,
   120     OF this f_Suc lb_0 refl ub_0 refl]
   121   show "?lb" and "?ub" unfolding minus_minus sum_eq
   122     by auto
   123 qed
   124 
   125 subsection {* Selectors for next even or odd number *}
   126 
   127 text {*
   128 
   129 The horner scheme computes alternating series. To get the upper and lower bounds we need to
   130 guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}.
   131 
   132 *}
   133 
   134 definition get_odd :: "nat \<Rightarrow> nat" where
   135   "get_odd n = (if odd n then n else (Suc n))"
   136 
   137 definition get_even :: "nat \<Rightarrow> nat" where
   138   "get_even n = (if even n then n else (Suc n))"
   139 
   140 lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto)
   141 lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto)
   142 lemma get_odd_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"
   143 proof (cases "odd n")
   144   case True hence "0 < n" by (rule odd_pos)
   145   from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto
   146   thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast
   147 next
   148   case False hence "odd (Suc n)" by auto
   149   thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast
   150 qed
   151 
   152 lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .
   153 lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto
   154 
   155 section "Power function"
   156 
   157 definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
   158 "float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
   159                       else if u < 0         then (u ^ n, l ^ n)
   160                                             else (0, (max (-l) u) ^ n))"
   161 
   162 lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {real l .. real u}"
   163   shows "x ^ n \<in> {real l1..real u1}"
   164 proof (cases "even n")
   165   case True
   166   show ?thesis
   167   proof (cases "0 < l")
   168     case True hence "odd n \<or> 0 < l" and "0 \<le> real l" unfolding less_float_def by auto
   169     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
   170     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
   171     thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
   172   next
   173     case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto
   174     show ?thesis
   175     proof (cases "u < 0")
   176       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
   177       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] 
   178 	unfolding power_minus_even[OF `even n`] by auto
   179       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
   180       ultimately show ?thesis using float_power by auto
   181     next
   182       case False
   183       have "\<bar>x\<bar> \<le> real (max (-l) u)"
   184       proof (cases "-l \<le> u")
   185 	case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto
   186       next
   187 	case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto
   188       qed
   189       hence x_abs: "\<bar>x\<bar> \<le> \<bar>real (max (-l) u)\<bar>" by auto
   190       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
   191       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
   192     qed
   193   qed
   194 next
   195   case False hence "odd n \<or> 0 < l" by auto
   196   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
   197   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
   198   thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto
   199 qed
   200 
   201 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"
   202   using float_power_bnds by auto
   203 
   204 section "Square root"
   205 
   206 text {*
   207 
   208 The square root computation is implemented as newton iteration. As first first step we use the
   209 nearest power of two greater than the square root.
   210 
   211 *}
   212 
   213 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
   214 "sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |
   215 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
   216                                   in Float 1 -1 * (y + float_divr prec x y))"
   217 
   218 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where
   219 "ub_sqrt prec x = (if 0 < x then (sqrt_iteration prec prec x)
   220               else if x < 0 then - lb_sqrt prec (- x)
   221                             else 0)" |
   222 "lb_sqrt prec x = (if 0 < x then (float_divl prec x (sqrt_iteration prec prec x))
   223               else if x < 0 then - ub_sqrt prec (- x)
   224                             else 0)"
   225 by pat_completeness auto
   226 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)
   227 
   228 declare lb_sqrt.simps[simp del]
   229 declare ub_sqrt.simps[simp del]
   230 
   231 lemma sqrt_ub_pos_pos_1:
   232   assumes "sqrt x < b" and "0 < b" and "0 < x"
   233   shows "sqrt x < (b + x / b)/2"
   234 proof -
   235   from assms have "0 < (b - sqrt x) ^ 2 " by simp
   236   also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra
   237   also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)
   238   finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption
   239   hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms
   240     by (simp add: field_simps power2_eq_square)
   241   thus ?thesis by (simp add: field_simps)
   242 qed
   243 
   244 lemma sqrt_iteration_bound: assumes "0 < real x"
   245   shows "sqrt (real x) < real (sqrt_iteration prec n x)"
   246 proof (induct n)
   247   case 0
   248   show ?case
   249   proof (cases x)
   250     case (Float m e)
   251     hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto
   252     hence "0 < sqrt (real m)" by auto
   253 
   254     have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto
   255 
   256     have "real x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"
   257       unfolding pow2_add pow2_int Float real_of_float_simp by auto
   258     also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"
   259     proof (rule mult_strict_right_mono, auto)
   260       show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
   261 	unfolding real_of_int_less_iff[of m, symmetric] by auto
   262     qed
   263     finally have "sqrt (real x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto
   264     also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"
   265     proof -
   266       let ?E = "e + bitlen m"
   267       have E_mod_pow: "pow2 (?E mod 2) < 4"
   268       proof (cases "?E mod 2 = 1")
   269 	case True thus ?thesis by auto
   270       next
   271 	case False
   272 	have "0 \<le> ?E mod 2" by auto
   273 	have "?E mod 2 < 2" by auto
   274 	from this[THEN zless_imp_add1_zle]
   275 	have "?E mod 2 \<le> 0" using False by auto
   276 	from xt1(5)[OF `0 \<le> ?E mod 2` this]
   277 	show ?thesis by auto
   278       qed
   279       hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto
   280       hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto
   281 
   282       have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto
   283       have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"
   284 	unfolding E_eq unfolding pow2_add ..
   285       also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"
   286 	unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto
   287       also have "\<dots> < pow2 (?E div 2) * 2"
   288 	by (rule mult_strict_left_mono, auto intro: E_mod_pow)
   289       also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto
   290       finally show ?thesis by auto
   291     qed
   292     finally show ?thesis
   293       unfolding Float sqrt_iteration.simps real_of_float_simp by auto
   294   qed
   295 next
   296   case (Suc n)
   297   let ?b = "sqrt_iteration prec n x"
   298   have "0 < sqrt (real x)" using `0 < real x` by auto
   299   also have "\<dots> < real ?b" using Suc .
   300   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
   301   also have "\<dots> \<le> (real ?b + real (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
   302   also have "\<dots> = real (Float 1 -1) * (real ?b + real (float_divr prec x ?b))" by auto
   303   finally show ?case unfolding sqrt_iteration.simps Let_def real_of_float_mult real_of_float_add right_distrib .
   304 qed
   305 
   306 lemma sqrt_iteration_lower_bound: assumes "0 < real x"
   307   shows "0 < real (sqrt_iteration prec n x)" (is "0 < ?sqrt")
   308 proof -
   309   have "0 < sqrt (real x)" using assms by auto
   310   also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .
   311   finally show ?thesis .
   312 qed
   313 
   314 lemma lb_sqrt_lower_bound: assumes "0 \<le> real x"
   315   shows "0 \<le> real (lb_sqrt prec x)"
   316 proof (cases "0 < x")
   317   case True hence "0 < real x" and "0 \<le> x" using `0 \<le> real x` unfolding less_float_def le_float_def by auto
   318   hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto 
   319   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
   320   thus ?thesis unfolding lb_sqrt.simps using True by auto
   321 next
   322   case False with `0 \<le> real x` have "real x = 0" unfolding less_float_def by auto
   323   thus ?thesis unfolding lb_sqrt.simps less_float_def by auto
   324 qed
   325 
   326 lemma bnds_sqrt':
   327   shows "sqrt (real x) \<in> { real (lb_sqrt prec x) .. real (ub_sqrt prec x) }"
   328 proof -
   329   { fix x :: float assume "0 < x"
   330     hence "0 < real x" and "0 \<le> real x" unfolding less_float_def by auto
   331     hence sqrt_gt0: "0 < sqrt (real x)" by auto
   332     hence sqrt_ub: "sqrt (real x) < real (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto
   333 
   334     have "real (float_divl prec x (sqrt_iteration prec prec x)) \<le>
   335           real x / real (sqrt_iteration prec prec x)" by (rule float_divl)
   336     also have "\<dots> < real x / sqrt (real x)"
   337       by (rule divide_strict_left_mono[OF sqrt_ub `0 < real x`
   338                mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])
   339     also have "\<dots> = sqrt (real x)" 
   340       unfolding inverse_eq_iff_eq[of _ "sqrt (real x)", symmetric]
   341 	        sqrt_divide_self_eq[OF `0 \<le> real x`, symmetric] by auto
   342     finally have "real (lb_sqrt prec x) \<le> sqrt (real x)"
   343       unfolding lb_sqrt.simps if_P[OF `0 < x`] by auto }
   344   note lb = this
   345 
   346   { fix x :: float assume "0 < x"
   347     hence "0 < real x" unfolding less_float_def by auto
   348     hence "0 < sqrt (real x)" by auto
   349     hence "sqrt (real x) < real (sqrt_iteration prec prec x)"
   350       using sqrt_iteration_bound by auto
   351     hence "sqrt (real x) \<le> real (ub_sqrt prec x)"
   352       unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }
   353   note ub = this
   354 
   355   show ?thesis
   356   proof (cases "0 < x")
   357     case True with lb ub show ?thesis by auto
   358   next case False show ?thesis
   359   proof (cases "real x = 0")
   360     case True thus ?thesis 
   361       by (auto simp add: less_float_def lb_sqrt.simps ub_sqrt.simps)
   362   next
   363     case False with `\<not> 0 < x` have "x < 0" and "0 < -x"
   364       by (auto simp add: less_float_def)
   365 
   366     with `\<not> 0 < x`
   367     show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]
   368       by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)
   369   qed qed
   370 qed
   371 
   372 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"
   373 proof ((rule allI) +, rule impI, erule conjE, rule conjI)
   374   fix x lx ux
   375   assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)"
   376     and x: "x \<in> {real lx .. real ux}"
   377   hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto
   378 
   379   have "sqrt (real lx) \<le> sqrt x" using x by auto
   380   from order_trans[OF _ this]
   381   show "real l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto
   382 
   383   have "sqrt x \<le> sqrt (real ux)" using x by auto
   384   from order_trans[OF this]
   385   show "sqrt x \<le> real u" unfolding u using bnds_sqrt'[of ux prec] by auto
   386 qed
   387 
   388 section "Arcus tangens and \<pi>"
   389 
   390 subsection "Compute arcus tangens series"
   391 
   392 text {*
   393 
   394 As first step we implement the computation of the arcus tangens series. This is only valid in the range
   395 @{term "{-1 :: real .. 1}"}. This is used to compute \<pi> and then the entire arcus tangens.
   396 
   397 *}
   398 
   399 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
   400 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
   401   "ub_arctan_horner prec 0 k x = 0"
   402 | "ub_arctan_horner prec (Suc n) k x = 
   403     (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"
   404 | "lb_arctan_horner prec 0 k x = 0"
   405 | "lb_arctan_horner prec (Suc n) k x = 
   406     (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"
   407 
   408 lemma arctan_0_1_bounds': assumes "0 \<le> real x" "real x \<le> 1" and "even n"
   409   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))}"
   410 proof -
   411   let "?c i" = "-1^i * (1 / real (i * 2 + 1) * real x ^ (i * 2 + 1))"
   412   let "?S n" = "\<Sum> i=0..<n. ?c i"
   413 
   414   have "0 \<le> real (x * x)" by auto
   415   from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto
   416   
   417   have "arctan (real x) \<in> { ?S n .. ?S (Suc n) }"
   418   proof (cases "real x = 0")
   419     case False
   420     hence "0 < real x" using `0 \<le> real x` by auto
   421     hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * real x ^ (0 * 2 + 1)" by auto 
   422 
   423     have "\<bar> real x \<bar> \<le> 1"  using `0 \<le> real x` `real x \<le> 1` by auto
   424     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`]
   425     show ?thesis unfolding arctan_series[OF `\<bar> real x \<bar> \<le> 1`] Suc_eq_plus1  .
   426   qed auto
   427   note arctan_bounds = this[unfolded atLeastAtMost_iff]
   428 
   429   have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto
   430 
   431   note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 
   432     and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
   433     and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", 
   434     OF `0 \<le> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
   435 
   436   { have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
   437       using bounds(1) `0 \<le> real x`
   438       unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
   439       unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
   440       by (auto intro!: mult_left_mono)
   441     also have "\<dots> \<le> arctan (real x)" using arctan_bounds ..
   442     finally have "real (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (real x)" . }
   443   moreover
   444   { have "arctan (real x) \<le> ?S (Suc n)" using arctan_bounds ..
   445     also have "\<dots> \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
   446       using bounds(2)[of "Suc n"] `0 \<le> real x`
   447       unfolding real_of_float_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
   448       unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
   449       by (auto intro!: mult_left_mono)
   450     finally have "arctan (real x) \<le> real (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
   451   ultimately show ?thesis by auto
   452 qed
   453 
   454 lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
   455   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))}"
   456 proof (cases "even n")
   457   case True
   458   obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto
   459   hence "even n'" unfolding even_Suc by auto
   460   have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
   461     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
   462   moreover
   463   have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"
   464     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
   465   ultimately show ?thesis by auto
   466 next
   467   case False hence "0 < n" by (rule odd_pos)
   468   from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..
   469   from False[unfolded this even_Suc]
   470   have "even n'" and "even (Suc (Suc n'))" by auto
   471   have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .
   472 
   473   have "arctan (real x) \<le> real (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"
   474     unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real x \<le> 1` `even n'`] by auto
   475   moreover
   476   have "real (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (real x)"
   477     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
   478   ultimately show ?thesis by auto
   479 qed
   480 
   481 subsection "Compute \<pi>"
   482 
   483 definition ub_pi :: "nat \<Rightarrow> float" where
   484   "ub_pi prec = (let A = rapprox_rat prec 1 5 ; 
   485                      B = lapprox_rat prec 1 239
   486                  in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) - 
   487                                                   B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
   488 
   489 definition lb_pi :: "nat \<Rightarrow> float" where
   490   "lb_pi prec = (let A = lapprox_rat prec 1 5 ; 
   491                      B = rapprox_rat prec 1 239
   492                  in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) - 
   493                                                   B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
   494 
   495 lemma pi_boundaries: "pi \<in> {real (lb_pi n) .. real (ub_pi n)}"
   496 proof -
   497   have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
   498 
   499   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
   500     let ?k = "rapprox_rat prec 1 k"
   501     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
   502       
   503     have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
   504     have "real ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]
   505       by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)
   506 
   507     have "1 / real k \<le> real ?k" using rapprox_rat[where x=1 and y=k] by auto
   508     hence "arctan (1 / real k) \<le> arctan (real ?k)" by (rule arctan_monotone')
   509     also have "\<dots> \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
   510       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
   511     finally have "arctan (1 / (real k)) \<le> real (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .
   512   } note ub_arctan = this
   513 
   514   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
   515     let ?k = "lapprox_rat prec 1 k"
   516     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
   517     have "1 / real k \<le> 1" using `1 < k` by auto
   518 
   519     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`)
   520     have "\<And>n. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)
   521 
   522     have "real ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto
   523 
   524     have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (real ?k)"
   525       using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
   526     also have "\<dots> \<le> arctan (1 / real k)" using `real ?k \<le> 1 / real k` by (rule arctan_monotone')
   527     finally have "real (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .
   528   } note lb_arctan = this
   529 
   530   have "pi \<le> real (ub_pi n)"
   531     unfolding ub_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub unfolding Float_num
   532     using lb_arctan[of 239] ub_arctan[of 5]
   533     by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
   534   moreover
   535   have "real (lb_pi n) \<le> pi"
   536     unfolding lb_pi_def machin_pi Let_def real_of_float_mult real_of_float_sub Float_num
   537     using lb_arctan[of 5] ub_arctan[of 239]
   538     by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)
   539   ultimately show ?thesis by auto
   540 qed
   541 
   542 subsection "Compute arcus tangens in the entire domain"
   543 
   544 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
   545   "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
   546                            lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
   547     in (if x < 0          then - ub_arctan prec (-x) else
   548         if x \<le> Float 1 -1 then lb_horner x else
   549         if x \<le> Float 1 1  then Float 1 1 * lb_horner (float_divl prec x (1 + ub_sqrt prec (1 + x * x)))
   550                           else (let inv = float_divr prec 1 x
   551                                 in if inv > 1 then 0
   552                                               else lb_pi prec * Float 1 -1 - ub_horner inv)))"
   553 
   554 | "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
   555                            ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
   556     in (if x < 0          then - lb_arctan prec (-x) else
   557         if x \<le> Float 1 -1 then ub_horner x else
   558         if x \<le> Float 1 1  then let y = float_divr prec x (1 + lb_sqrt prec (1 + x * x))
   559                                in if y > 1 then ub_pi prec * Float 1 -1
   560                                            else Float 1 1 * ub_horner y
   561                           else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"
   562 by pat_completeness auto
   563 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)
   564 
   565 declare ub_arctan_horner.simps[simp del]
   566 declare lb_arctan_horner.simps[simp del]
   567 
   568 lemma lb_arctan_bound': assumes "0 \<le> real x"
   569   shows "real (lb_arctan prec x) \<le> arctan (real x)"
   570 proof -
   571   have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
   572   let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
   573     and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
   574 
   575   show ?thesis
   576   proof (cases "x \<le> Float 1 -1")
   577     case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
   578     show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
   579       using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
   580   next
   581     case False hence "0 < real x" unfolding le_float_def Float_num by auto
   582     let ?R = "1 + sqrt (1 + real x * real x)"
   583     let ?fR = "1 + ub_sqrt prec (1 + x * x)"
   584     let ?DIV = "float_divl prec x ?fR"
   585 
   586     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
   587     hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
   588 
   589     have "sqrt (real (1 + x * x)) \<le> real (ub_sqrt prec (1 + x * x))"
   590       using bnds_sqrt'[of "1 + x * x"] by auto
   591 
   592     hence "?R \<le> real ?fR" by auto
   593     hence "0 < ?fR" and "0 < real ?fR" unfolding less_float_def using `0 < ?R` by auto
   594 
   595     have monotone: "real (float_divl prec x ?fR) \<le> real x / ?R"
   596     proof -
   597       have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)
   598       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]])
   599       finally show ?thesis .
   600     qed
   601 
   602     show ?thesis
   603     proof (cases "x \<le> Float 1 1")
   604       case True
   605 
   606       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
   607       also have "\<dots> \<le> real (ub_sqrt prec (1 + x * x))"
   608 	using bnds_sqrt'[of "1 + x * x"] by auto
   609       finally have "real x \<le> real ?fR" by auto
   610       moreover have "real ?DIV \<le> real x / real ?fR" by (rule float_divl)
   611       ultimately have "real ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < real ?fR`, symmetric] by auto
   612 
   613       have "0 \<le> real ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto
   614 
   615       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
   616 	using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
   617       also have "\<dots> \<le> 2 * arctan (real x / ?R)"
   618 	using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
   619       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 . 
   620       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] .
   621     next
   622       case False
   623       hence "2 < real x" unfolding le_float_def Float_num by auto
   624       hence "1 \<le> real x" by auto
   625 
   626       let "?invx" = "float_divr prec 1 x"
   627       have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
   628 
   629       show ?thesis
   630       proof (cases "1 < ?invx")
   631 	case True
   632 	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] 
   633 	  using `0 \<le> arctan (real x)` by auto
   634       next
   635 	case False
   636 	hence "real ?invx \<le> 1" unfolding less_float_def by auto
   637 	have "0 \<le> real ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> real x`)
   638 
   639 	have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto
   640 
   641 	have "arctan (1 / real x) \<le> arctan (real ?invx)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divr)
   642 	also have "\<dots> \<le> real (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
   643 	finally have "pi / 2 - real (?ub_horner ?invx) \<le> arctan (real x)"
   644 	  using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`]
   645 	  unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
   646 	moreover
   647 	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
   648 	ultimately
   649 	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]
   650 	  by auto
   651       qed
   652     qed
   653   qed
   654 qed
   655 
   656 lemma ub_arctan_bound': assumes "0 \<le> real x"
   657   shows "arctan (real x) \<le> real (ub_arctan prec x)"
   658 proof -
   659   have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> real x` by auto
   660 
   661   let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
   662     and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
   663 
   664   show ?thesis
   665   proof (cases "x \<le> Float 1 -1")
   666     case True hence "real x \<le> 1" unfolding le_float_def Float_num by auto
   667     show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
   668       using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
   669   next
   670     case False hence "0 < real x" unfolding le_float_def Float_num by auto
   671     let ?R = "1 + sqrt (1 + real x * real x)"
   672     let ?fR = "1 + lb_sqrt prec (1 + x * x)"
   673     let ?DIV = "float_divr prec x ?fR"
   674 
   675     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
   676     hence "0 \<le> real (1 + x*x)" by auto
   677 
   678     hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
   679 
   680     have "real (lb_sqrt prec (1 + x * x)) \<le> sqrt (real (1 + x * x))"
   681       using bnds_sqrt'[of "1 + x * x"] by auto
   682     hence "real ?fR \<le> ?R" by auto
   683     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)`])
   684 
   685     have monotone: "real x / ?R \<le> real (float_divr prec x ?fR)"
   686     proof -
   687       from divide_left_mono[OF `real ?fR \<le> ?R` `0 \<le> real x` mult_pos_pos[OF divisor_gt0 `0 < real ?fR`]]
   688       have "real x / ?R \<le> real x / real ?fR" .
   689       also have "\<dots> \<le> real ?DIV" by (rule float_divr)
   690       finally show ?thesis .
   691     qed
   692 
   693     show ?thesis
   694     proof (cases "x \<le> Float 1 1")
   695       case True
   696       show ?thesis
   697       proof (cases "?DIV > 1")
   698 	case True
   699 	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
   700 	from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le]
   701 	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] .
   702       next
   703 	case False
   704 	hence "real ?DIV \<le> 1" unfolding less_float_def by auto
   705 
   706 	have "0 \<le> real x / ?R" using `0 \<le> real x` `0 < ?R` unfolding real_0_le_divide_iff by auto
   707 	hence "0 \<le> real ?DIV" using monotone by (rule order_trans)
   708 
   709 	have "arctan (real x) = 2 * arctan (real x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .
   710 	also have "\<dots> \<le> 2 * arctan (real ?DIV)"
   711 	  using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
   712 	also have "\<dots> \<le> real (Float 1 1 * ?ub_horner ?DIV)" unfolding real_of_float_mult[of "Float 1 1"] Float_num
   713 	  using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
   714 	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] .
   715       qed
   716     next
   717       case False
   718       hence "2 < real x" unfolding le_float_def Float_num by auto
   719       hence "1 \<le> real x" by auto
   720       hence "0 < real x" by auto
   721       hence "0 < x" unfolding less_float_def by auto
   722 
   723       let "?invx" = "float_divl prec 1 x"
   724       have "0 \<le> arctan (real x)" using arctan_monotone'[OF `0 \<le> real x`] using arctan_tan[of 0, unfolded tan_zero] by auto
   725 
   726       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`])
   727       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`)
   728 
   729       have "1 / real x \<noteq> 0" and "0 < 1 / real x" using `0 < real x` by auto
   730 
   731       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
   732       also have "\<dots> \<le> arctan (1 / real x)" unfolding real_of_float_1[symmetric] by (rule arctan_monotone', rule float_divl)
   733       finally have "arctan (real x) \<le> pi / 2 - real (?lb_horner ?invx)"
   734 	using `0 \<le> arctan (real x)` arctan_inverse[OF `1 / real x \<noteq> 0`] 
   735 	unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
   736       moreover
   737       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
   738       ultimately
   739       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]
   740 	by auto
   741     qed
   742   qed
   743 qed
   744 
   745 lemma arctan_boundaries:
   746   "arctan (real x) \<in> {real (lb_arctan prec x) .. real (ub_arctan prec x)}"
   747 proof (cases "0 \<le> x")
   748   case True hence "0 \<le> real x" unfolding le_float_def by auto
   749   show ?thesis using ub_arctan_bound'[OF `0 \<le> real x`] lb_arctan_bound'[OF `0 \<le> real x`] unfolding atLeastAtMost_iff by auto
   750 next
   751   let ?mx = "-x"
   752   case False hence "x < 0" and "0 \<le> real ?mx" unfolding le_float_def less_float_def by auto
   753   hence bounds: "real (lb_arctan prec ?mx) \<le> arctan (real ?mx) \<and> arctan (real ?mx) \<le> real (ub_arctan prec ?mx)"
   754     using ub_arctan_bound'[OF `0 \<le> real ?mx`] lb_arctan_bound'[OF `0 \<le> real ?mx`] by auto
   755   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`]
   756     unfolding atLeastAtMost_iff using bounds[unfolded real_of_float_minus arctan_minus] by auto
   757 qed
   758 
   759 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"
   760 proof (rule allI, rule allI, rule allI, rule impI)
   761   fix x lx ux
   762   assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {real lx .. real ux}"
   763   hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto
   764 
   765   { from arctan_boundaries[of lx prec, unfolded l]
   766     have "real l \<le> arctan (real lx)" by (auto simp del: lb_arctan.simps)
   767     also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')
   768     finally have "real l \<le> arctan x" .
   769   } moreover
   770   { have "arctan x \<le> arctan (real ux)" using x by (auto intro: arctan_monotone')
   771     also have "\<dots> \<le> real u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
   772     finally have "arctan x \<le> real u" .
   773   } ultimately show "real l \<le> arctan x \<and> arctan x \<le> real u" ..
   774 qed
   775 
   776 section "Sinus and Cosinus"
   777 
   778 subsection "Compute the cosinus and sinus series"
   779 
   780 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
   781 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
   782   "ub_sin_cos_aux prec 0 i k x = 0"
   783 | "ub_sin_cos_aux prec (Suc n) i k x = 
   784     (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
   785 | "lb_sin_cos_aux prec 0 i k x = 0"
   786 | "lb_sin_cos_aux prec (Suc n) i k x = 
   787     (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
   788 
   789 lemma cos_aux:
   790   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")
   791   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")
   792 proof -
   793   have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
   794   let "?f n" = "fact (2 * n)"
   795 
   796   { fix n 
   797     have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
   798     have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ n) 1 + 1)"
   799       unfolding F by auto } note f_eq = this
   800     
   801   from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 
   802     OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
   803   show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"])
   804 qed
   805 
   806 lemma cos_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"
   807   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))}"
   808 proof (cases "real x = 0")
   809   case False hence "real x \<noteq> 0" by auto
   810   hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
   811   have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
   812     using mult_pos_pos[where a="real x" and b="real x"] by auto
   813 
   814   { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^ (2 * i))
   815     = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")
   816   proof -
   817     have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto
   818     also have "\<dots> = 
   819       (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
   820     also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"
   821       unfolding sum_split_even_odd ..
   822     also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"
   823       by (rule setsum_cong2) auto
   824     finally show ?thesis by assumption
   825   qed } note morph_to_if_power = this
   826 
   827 
   828   { fix n :: nat assume "0 < n"
   829     hence "0 < 2 * n" by auto
   830     obtain t where "0 < t" and "t < real x" and
   831       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) 
   832       + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (real x)^(2*n)" 
   833       (is "_ = ?SUM + ?rest / ?fact * ?pow")
   834       using Maclaurin_cos_expansion2[OF `0 < real x` `0 < 2 * n`] by auto
   835 
   836     have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto
   837     also have "\<dots> = cos (t + real n * pi)"  using cos_add by auto
   838     also have "\<dots> = ?rest" by auto
   839     finally have "cos t * -1^n = ?rest" .
   840     moreover
   841     have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto
   842     hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
   843     ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
   844 
   845     have "0 < ?fact" by auto
   846     have "0 < ?pow" using `0 < real x` by auto
   847 
   848     {
   849       assume "even n"
   850       have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"
   851 	unfolding morph_to_if_power[symmetric] using cos_aux by auto 
   852       also have "\<dots> \<le> cos (real x)"
   853       proof -
   854 	from even[OF `even n`] `0 < ?fact` `0 < ?pow`
   855 	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
   856 	thus ?thesis unfolding cos_eq by auto
   857       qed
   858       finally have "real (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (real x)" .
   859     } note lb = this
   860 
   861     {
   862       assume "odd n"
   863       have "cos (real x) \<le> ?SUM"
   864       proof -
   865 	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
   866 	have "0 \<le> (- ?rest) / ?fact * ?pow"
   867 	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
   868 	thus ?thesis unfolding cos_eq by auto
   869       qed
   870       also have "\<dots> \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))"
   871 	unfolding morph_to_if_power[symmetric] using cos_aux by auto
   872       finally have "cos (real x) \<le> real (ub_sin_cos_aux prec n 1 1 (x * x))" .
   873     } note ub = this and lb
   874   } note ub = this(1) and lb = this(2)
   875 
   876   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] .
   877   moreover have "real (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (real x)" 
   878   proof (cases "0 < get_even n")
   879     case True show ?thesis using lb[OF True get_even] .
   880   next
   881     case False
   882     hence "get_even n = 0" by auto
   883     have "- (pi / 2) \<le> real x" by (rule order_trans[OF _ `0 < real x`[THEN less_imp_le]], auto)
   884     with `real x \<le> pi / 2`
   885     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
   886   qed
   887   ultimately show ?thesis by auto
   888 next
   889   case True
   890   show ?thesis
   891   proof (cases "n = 0")
   892     case True 
   893     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
   894   next
   895     case False with not0_implies_Suc obtain m where "n = Suc m" by blast
   896     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)
   897   qed
   898 qed
   899 
   900 lemma sin_aux: assumes "0 \<le> real x"
   901   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")
   902   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")
   903 proof -
   904   have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto
   905   let "?f n" = "fact (2 * n + 1)"
   906 
   907   { fix n 
   908     have F: "\<And>m. ((\<lambda>i. i + 2) ^^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)
   909     have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ n) 2 + 1)"
   910       unfolding F by auto } note f_eq = this
   911     
   912   from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
   913     OF `0 \<le> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
   914   show "?lb" and "?ub" using `0 \<le> real x` unfolding real_of_float_mult
   915     unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]
   916     unfolding real_mult_commute
   917     by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real x"])
   918 qed
   919 
   920 lemma sin_boundaries: assumes "0 \<le> real x" and "real x \<le> pi / 2"
   921   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))}"
   922 proof (cases "real x = 0")
   923   case False hence "real x \<noteq> 0" by auto
   924   hence "0 < x" and "0 < real x" using `0 \<le> real x` unfolding less_float_def by auto
   925   have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0
   926     using mult_pos_pos[where a="real x" and b="real x"] by auto
   927 
   928   { 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))
   929     = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")
   930     proof -
   931       have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto
   932       have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto
   933       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)"
   934 	unfolding sum_split_even_odd ..
   935       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)"
   936 	by (rule setsum_cong2) auto
   937       finally show ?thesis by assumption
   938     qed } note setsum_morph = this
   939 
   940   { fix n :: nat assume "0 < n"
   941     hence "0 < 2 * n + 1" by auto
   942     obtain t where "0 < t" and "t < real x" and
   943       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) 
   944       + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (real x)^(2*n + 1)" 
   945       (is "_ = ?SUM + ?rest / ?fact * ?pow")
   946       using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < real x`] by auto
   947 
   948     have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto
   949     moreover
   950     have "t \<le> pi / 2" using `t < real x` and `real x \<le> pi / 2` by auto
   951     hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto
   952     ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto
   953 
   954     have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)
   955     have "0 < ?pow" using `0 < real x` by (rule zero_less_power)
   956 
   957     {
   958       assume "even n"
   959       have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> 
   960             (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (real x) ^ i)"
   961 	using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
   962       also have "\<dots> \<le> ?SUM" by auto
   963       also have "\<dots> \<le> sin (real x)"
   964       proof -
   965 	from even[OF `even n`] `0 < ?fact` `0 < ?pow`
   966 	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
   967 	thus ?thesis unfolding sin_eq by auto
   968       qed
   969       finally have "real (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (real x)" .
   970     } note lb = this
   971 
   972     {
   973       assume "odd n"
   974       have "sin (real x) \<le> ?SUM"
   975       proof -
   976 	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]
   977 	have "0 \<le> (- ?rest) / ?fact * ?pow"
   978 	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
   979 	thus ?thesis unfolding sin_eq by auto
   980       qed
   981       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)"
   982 	 by auto
   983       also have "\<dots> \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" 
   984 	using sin_aux[OF `0 \<le> real x`] unfolding setsum_morph[symmetric] by auto
   985       finally have "sin (real x) \<le> real (x * ub_sin_cos_aux prec n 2 1 (x * x))" .
   986     } note ub = this and lb
   987   } note ub = this(1) and lb = this(2)
   988 
   989   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] .
   990   moreover have "real (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (real x)" 
   991   proof (cases "0 < get_even n")
   992     case True show ?thesis using lb[OF True get_even] .
   993   next
   994     case False
   995     hence "get_even n = 0" by auto
   996     with `real x \<le> pi / 2` `0 \<le> real x`
   997     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
   998   qed
   999   ultimately show ?thesis by auto
  1000 next
  1001   case True
  1002   show ?thesis
  1003   proof (cases "n = 0")
  1004     case True 
  1005     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
  1006   next
  1007     case False with not0_implies_Suc obtain m where "n = Suc m" by blast
  1008     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)
  1009   qed
  1010 qed
  1011 
  1012 subsection "Compute the cosinus in the entire domain"
  1013 
  1014 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1015 "lb_cos prec x = (let
  1016     horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
  1017     half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
  1018   in if x < Float 1 -1 then horner x
  1019 else if x < 1          then half (horner (x * Float 1 -1))
  1020                        else half (half (horner (x * Float 1 -2))))"
  1021 
  1022 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1023 "ub_cos prec x = (let
  1024     horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
  1025     half = \<lambda> x. Float 1 1 * x * x - 1
  1026   in if x < Float 1 -1 then horner x
  1027 else if x < 1          then half (horner (x * Float 1 -1))
  1028                        else half (half (horner (x * Float 1 -2))))"
  1029 
  1030 lemma lb_cos: assumes "0 \<le> real x" and "real x \<le> pi"
  1031   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) }")
  1032 proof -
  1033   { fix x :: real
  1034     have "cos x = cos (x / 2 + x / 2)" by auto
  1035     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"
  1036       unfolding cos_add by auto
  1037     also have "\<dots> = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra
  1038     finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" .
  1039   } note x_half = this[symmetric]
  1040 
  1041   have "\<not> x < 0" using `0 \<le> real x` unfolding less_float_def by auto
  1042   let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
  1043   let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
  1044   let "?ub_half x" = "Float 1 1 * x * x - 1"
  1045   let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
  1046 
  1047   show ?thesis
  1048   proof (cases "x < Float 1 -1")
  1049     case True hence "real x \<le> pi / 2" unfolding less_float_def using pi_ge_two by auto
  1050     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
  1051       using cos_boundaries[OF `0 \<le> real x` `real x \<le> pi / 2`] .
  1052   next
  1053     case False
  1054     { fix y x :: float let ?x2 = "real (x * Float 1 -1)"
  1055       assume "real y \<le> cos ?x2" and "-pi \<le> real x" and "real x \<le> pi"
  1056       hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
  1057       hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
  1058 
  1059       have "real (?lb_half y) \<le> cos (real x)"
  1060       proof (cases "y < 0")
  1061 	case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto
  1062       next
  1063 	case False
  1064 	hence "0 \<le> real y" unfolding less_float_def by auto
  1065 	from mult_mono[OF `real y \<le> cos ?x2` `real y \<le> cos ?x2` `0 \<le> cos ?x2` this]
  1066 	have "real y * real y \<le> cos ?x2 * cos ?x2" .
  1067 	hence "2 * real y * real y \<le> 2 * cos ?x2 * cos ?x2" by auto
  1068 	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
  1069 	thus ?thesis unfolding if_not_P[OF False] x_half Float_num real_of_float_mult real_of_float_sub by auto
  1070       qed
  1071     } note lb_half = this
  1072 
  1073     { fix y x :: float let ?x2 = "real (x * Float 1 -1)"
  1074       assume ub: "cos ?x2 \<le> real y" and "- pi \<le> real x" and "real x \<le> pi"
  1075       hence "- (pi / 2) \<le> ?x2" and "?x2 \<le> pi / 2" using pi_ge_two unfolding real_of_float_mult Float_num by auto
  1076       hence "0 \<le> cos ?x2" by (rule cos_ge_zero)
  1077 
  1078       have "cos (real x) \<le> real (?ub_half y)"
  1079       proof -
  1080 	have "0 \<le> real y" using `0 \<le> cos ?x2` ub by (rule order_trans)
  1081 	from mult_mono[OF ub ub this `0 \<le> cos ?x2`]
  1082 	have "cos ?x2 * cos ?x2 \<le> real y * real y" .
  1083 	hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real y * real y" by auto
  1084 	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
  1085 	thus ?thesis unfolding x_half real_of_float_mult Float_num real_of_float_sub by auto
  1086       qed
  1087     } note ub_half = this
  1088 
  1089     let ?x2 = "x * Float 1 -1"
  1090     let ?x4 = "x * Float 1 -1 * Float 1 -1"
  1091 
  1092     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)
  1093 
  1094     show ?thesis
  1095     proof (cases "x < 1")
  1096       case True hence "real x \<le> 1" unfolding less_float_def by auto
  1097       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
  1098       from cos_boundaries[OF this]
  1099       have lb: "real (?lb_horner ?x2) \<le> ?cos ?x2" and ub: "?cos ?x2 \<le> real (?ub_horner ?x2)" by auto
  1100 
  1101       have "real (?lb x) \<le> ?cos x"
  1102       proof -
  1103 	from lb_half[OF lb `-pi \<le> real x` `real x \<le> pi`]
  1104 	show ?thesis unfolding lb_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto
  1105       qed
  1106       moreover have "?cos x \<le> real (?ub x)"
  1107       proof -
  1108 	from ub_half[OF ub `-pi \<le> real x` `real x \<le> pi`]
  1109 	show ?thesis unfolding ub_cos_def[where x=x] Let_def using `\<not> x < 0` `\<not> x < Float 1 -1` `x < 1` by auto 
  1110       qed
  1111       ultimately show ?thesis by auto
  1112     next
  1113       case False
  1114       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
  1115       from cos_boundaries[OF this]
  1116       have lb: "real (?lb_horner ?x4) \<le> ?cos ?x4" and ub: "?cos ?x4 \<le> real (?ub_horner ?x4)" by auto
  1117 
  1118       have eq_4: "?x2 * Float 1 -1 = x * Float 1 -2" by (cases x, auto simp add: times_float.simps)
  1119 
  1120       have "real (?lb x) \<le> ?cos x"
  1121       proof -
  1122 	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
  1123 	from lb_half[OF lb_half[OF lb this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]
  1124 	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 .
  1125       qed
  1126       moreover have "?cos x \<le> real (?ub x)"
  1127       proof -
  1128 	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
  1129 	from ub_half[OF ub_half[OF ub this] `-pi \<le> real x` `real x \<le> pi`, unfolded eq_4]
  1130 	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 .
  1131       qed
  1132       ultimately show ?thesis by auto
  1133     qed
  1134   qed
  1135 qed
  1136 
  1137 lemma lb_cos_minus: assumes "-pi \<le> real x" and "real x \<le> 0"
  1138   shows "cos (real (-x)) \<in> {real (lb_cos prec (-x)) .. real (ub_cos prec (-x))}"
  1139 proof -
  1140   have "0 \<le> real (-x)" and "real (-x) \<le> pi" using `-pi \<le> real x` `real x \<le> 0` by auto
  1141   from lb_cos[OF this] show ?thesis .
  1142 qed
  1143 
  1144 definition bnds_cos :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
  1145 "bnds_cos prec lx ux = (let
  1146     lpi = round_down prec (lb_pi prec) ;
  1147     upi = round_up prec (ub_pi prec) ;
  1148     k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ;
  1149     lx = lx - k * 2 * (if k < 0 then lpi else upi) ;
  1150     ux = ux - k * 2 * (if k < 0 then upi else lpi)
  1151   in   if - lpi \<le> lx \<and> ux \<le> 0    then (lb_cos prec (-lx), ub_cos prec (-ux))
  1152   else if 0 \<le> lx \<and> ux \<le> lpi      then (lb_cos prec ux, ub_cos prec lx)
  1153   else if - lpi \<le> lx \<and> ux \<le> lpi  then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0)
  1154   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))))
  1155   else if -2 * lpi \<le> lx \<and> ux \<le> 0 then (Float -1 0, max (ub_cos prec (lx + 2 * lpi)) (ub_cos prec (-ux)))
  1156                                  else (Float -1 0, Float 1 0))"
  1157 
  1158 lemma floor_int:
  1159   obtains k :: int where "real k = real (floor_fl f)"
  1160 proof -
  1161   assume *: "\<And> k :: int. real k = real (floor_fl f) \<Longrightarrow> thesis"
  1162   obtain m e where fl: "Float m e = floor_fl f" by (cases "floor_fl f", auto)
  1163   from floor_pos_exp[OF this]
  1164   have "real (m* 2^(nat e)) = real (floor_fl f)"
  1165     by (auto simp add: fl[symmetric] real_of_float_def pow2_def)
  1166   from *[OF this] show thesis by blast
  1167 qed
  1168 
  1169 lemma float_remove_real_numeral[simp]: "real (number_of k :: float) = number_of k"
  1170 proof -
  1171   have "real (number_of k :: float) = real k"
  1172     unfolding number_of_float_def real_of_float_def pow2_def by auto
  1173   also have "\<dots> = real (number_of k :: int)"
  1174     by (simp add: number_of_is_id)
  1175   finally show ?thesis by auto
  1176 qed
  1177 
  1178 lemma cos_periodic_nat[simp]: fixes n :: nat shows "cos (x + real n * 2 * pi) = cos x"
  1179 proof (induct n arbitrary: x)
  1180   case (Suc n)
  1181   have split_pi_off: "x + real (Suc n) * 2 * pi = (x + real n * 2 * pi) + 2 * pi"
  1182     unfolding Suc_eq_plus1 real_of_nat_add real_of_one real_add_mult_distrib by auto
  1183   show ?case unfolding split_pi_off using Suc by auto
  1184 qed auto
  1185 
  1186 lemma cos_periodic_int[simp]: fixes i :: int shows "cos (x + real i * 2 * pi) = cos x"
  1187 proof (cases "0 \<le> i")
  1188   case True hence i_nat: "real i = real (nat i)" by auto
  1189   show ?thesis unfolding i_nat by auto
  1190 next
  1191   case False hence i_nat: "real i = - real (nat (-i))" by auto
  1192   have "cos x = cos (x + real i * 2 * pi - real i * 2 * pi)" by auto
  1193   also have "\<dots> = cos (x + real i * 2 * pi)"
  1194     unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat)
  1195   finally show ?thesis by auto
  1196 qed
  1197 
  1198 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"
  1199 proof ((rule allI | rule impI | erule conjE) +)
  1200   fix x lx ux
  1201   assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \<in> {real lx .. real ux}"
  1202 
  1203   let ?lpi = "round_down prec (lb_pi prec)"
  1204   let ?upi = "round_up prec (ub_pi prec)"
  1205   let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))"
  1206   let ?lx = "lx - ?k * 2 * (if ?k < 0 then ?lpi else ?upi)"
  1207   let ?ux = "ux - ?k * 2 * (if ?k < 0 then ?upi else ?lpi)"
  1208 
  1209   obtain k :: int where k: "real k = real ?k" using floor_int .
  1210 
  1211   have upi: "pi \<le> real ?upi" and lpi: "real ?lpi \<le> pi"
  1212     using round_up[of "ub_pi prec" prec] pi_boundaries[of prec]
  1213           round_down[of prec "lb_pi prec"] by auto
  1214   hence "real ?lx \<le> x - real k * 2 * pi \<and> x - real k * 2 * pi \<le> real ?ux"
  1215     using x by (cases "k = 0") (auto intro!: add_mono
  1216                 simp add: real_diff_def k[symmetric] less_float_def)
  1217   note lx = this[THEN conjunct1] and ux = this[THEN conjunct2]
  1218   hence lx_less_ux: "real ?lx \<le> real ?ux" by (rule order_trans)
  1219 
  1220   { assume "- ?lpi \<le> ?lx" and x_le_0: "x - real k * 2 * pi \<le> 0"
  1221     with lpi[THEN le_imp_neg_le] lx
  1222     have pi_lx: "- pi \<le> real ?lx" and lx_0: "real ?lx \<le> 0"
  1223       by (simp_all add: le_float_def)
  1224 
  1225     have "real (lb_cos prec (- ?lx)) \<le> cos (real (- ?lx))"
  1226       using lb_cos_minus[OF pi_lx lx_0] by simp
  1227     also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"
  1228       using cos_monotone_minus_pi_0'[OF pi_lx lx x_le_0]
  1229       by (simp only: real_of_float_minus real_of_int_minus
  1230 	cos_minus real_diff_def mult_minus_left)
  1231     finally have "real (lb_cos prec (- ?lx)) \<le> cos x"
  1232       unfolding cos_periodic_int . }
  1233   note negative_lx = this
  1234 
  1235   { assume "0 \<le> ?lx" and pi_x: "x - real k * 2 * pi \<le> pi"
  1236     with lx
  1237     have pi_lx: "real ?lx \<le> pi" and lx_0: "0 \<le> real ?lx"
  1238       by (auto simp add: le_float_def)
  1239 
  1240     have "cos (x + real (-k) * 2 * pi) \<le> cos (real ?lx)"
  1241       using cos_monotone_0_pi'[OF lx_0 lx pi_x]
  1242       by (simp only: real_of_float_minus real_of_int_minus
  1243 	cos_minus real_diff_def mult_minus_left)
  1244     also have "\<dots> \<le> real (ub_cos prec ?lx)"
  1245       using lb_cos[OF lx_0 pi_lx] by simp
  1246     finally have "cos x \<le> real (ub_cos prec ?lx)"
  1247       unfolding cos_periodic_int . }
  1248   note positive_lx = this
  1249 
  1250   { assume pi_x: "- pi \<le> x - real k * 2 * pi" and "?ux \<le> 0"
  1251     with ux
  1252     have pi_ux: "- pi \<le> real ?ux" and ux_0: "real ?ux \<le> 0"
  1253       by (simp_all add: le_float_def)
  1254 
  1255     have "cos (x + real (-k) * 2 * pi) \<le> cos (real (- ?ux))"
  1256       using cos_monotone_minus_pi_0'[OF pi_x ux ux_0]
  1257       by (simp only: real_of_float_minus real_of_int_minus
  1258 	  cos_minus real_diff_def mult_minus_left)
  1259     also have "\<dots> \<le> real (ub_cos prec (- ?ux))"
  1260       using lb_cos_minus[OF pi_ux ux_0, of prec] by simp
  1261     finally have "cos x \<le> real (ub_cos prec (- ?ux))"
  1262       unfolding cos_periodic_int . }
  1263   note negative_ux = this
  1264 
  1265   { assume "?ux \<le> ?lpi" and x_ge_0: "0 \<le> x - real k * 2 * pi"
  1266     with lpi ux
  1267     have pi_ux: "real ?ux \<le> pi" and ux_0: "0 \<le> real ?ux"
  1268       by (simp_all add: le_float_def)
  1269 
  1270     have "real (lb_cos prec ?ux) \<le> cos (real ?ux)"
  1271       using lb_cos[OF ux_0 pi_ux] by simp
  1272     also have "\<dots> \<le> cos (x + real (-k) * 2 * pi)"
  1273       using cos_monotone_0_pi'[OF x_ge_0 ux pi_ux]
  1274       by (simp only: real_of_float_minus real_of_int_minus
  1275 	cos_minus real_diff_def mult_minus_left)
  1276     finally have "real (lb_cos prec ?ux) \<le> cos x"
  1277       unfolding cos_periodic_int . }
  1278   note positive_ux = this
  1279 
  1280   show "real l \<le> cos x \<and> cos x \<le> real u"
  1281   proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> 0")
  1282     case True with bnds
  1283     have l: "l = lb_cos prec (-?lx)"
  1284       and u: "u = ub_cos prec (-?ux)"
  1285       by (auto simp add: bnds_cos_def Let_def)
  1286 
  1287     from True lpi[THEN le_imp_neg_le] lx ux
  1288     have "- pi \<le> x - real k * 2 * pi"
  1289       and "x - real k * 2 * pi \<le> 0"
  1290       by (auto simp add: le_float_def)
  1291     with True negative_ux negative_lx
  1292     show ?thesis unfolding l u by simp
  1293   next case False note 1 = this show ?thesis
  1294   proof (cases "0 \<le> ?lx \<and> ?ux \<le> ?lpi")
  1295     case True with bnds 1
  1296     have l: "l = lb_cos prec ?ux"
  1297       and u: "u = ub_cos prec ?lx"
  1298       by (auto simp add: bnds_cos_def Let_def)
  1299 
  1300     from True lpi lx ux
  1301     have "0 \<le> x - real k * 2 * pi"
  1302       and "x - real k * 2 * pi \<le> pi"
  1303       by (auto simp add: le_float_def)
  1304     with True positive_ux positive_lx
  1305     show ?thesis unfolding l u by simp
  1306   next case False note 2 = this show ?thesis
  1307   proof (cases "- ?lpi \<le> ?lx \<and> ?ux \<le> ?lpi")
  1308     case True note Cond = this with bnds 1 2
  1309     have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)"
  1310       and u: "u = Float 1 0"
  1311       by (auto simp add: bnds_cos_def Let_def)
  1312 
  1313     show ?thesis unfolding u l using negative_lx positive_ux Cond
  1314       by (cases "x - real k * 2 * pi < 0", simp_all add: real_of_float_min)
  1315   next case False note 3 = this show ?thesis
  1316   proof (cases "0 \<le> ?lx \<and> ?ux \<le> 2 * ?lpi")
  1317     case True note Cond = this with bnds 1 2 3
  1318     have l: "l = Float -1 0"
  1319       and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))"
  1320       by (auto simp add: bnds_cos_def Let_def)
  1321 
  1322     have "cos x \<le> real u"
  1323     proof (cases "x - real k * 2 * pi < pi")
  1324       case True hence "x - real k * 2 * pi \<le> pi" by simp
  1325       from positive_lx[OF Cond[THEN conjunct1] this]
  1326       show ?thesis unfolding u by (simp add: real_of_float_max)
  1327     next
  1328       case False hence "pi \<le> x - real k * 2 * pi" by simp
  1329       hence pi_x: "- pi \<le> x - real k * 2 * pi - 2 * pi" by simp
  1330 
  1331       have "real ?ux \<le> 2 * pi" using Cond lpi by (auto simp add: le_float_def)
  1332       hence "x - real k * 2 * pi - 2 * pi \<le> 0" using ux by simp
  1333 
  1334       have ux_0: "real (?ux - 2 * ?lpi) \<le> 0"
  1335 	using Cond by (auto simp add: le_float_def)
  1336 
  1337       from 2 and Cond have "\<not> ?ux \<le> ?lpi" by auto
  1338       hence "- ?lpi \<le> ?ux - 2 * ?lpi" by (auto simp add: le_float_def)
  1339       hence pi_ux: "- pi \<le> real (?ux - 2 * ?lpi)"
  1340 	using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
  1341 
  1342       have x_le_ux: "x - real k * 2 * pi - 2 * pi \<le> real (?ux - 2 * ?lpi)"
  1343 	using ux lpi by auto
  1344 
  1345       have "cos x = cos (x + real (-k) * 2 * pi + real (-1 :: int) * 2 * pi)"
  1346 	unfolding cos_periodic_int ..
  1347       also have "\<dots> \<le> cos (real (?ux - 2 * ?lpi))"
  1348 	using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0]
  1349 	by (simp only: real_of_float_minus real_of_int_minus real_of_one
  1350 	    number_of_Min real_diff_def mult_minus_left real_mult_1)
  1351       also have "\<dots> = cos (real (- (?ux - 2 * ?lpi)))"
  1352 	unfolding real_of_float_minus cos_minus ..
  1353       also have "\<dots> \<le> real (ub_cos prec (- (?ux - 2 * ?lpi)))"
  1354 	using lb_cos_minus[OF pi_ux ux_0] by simp
  1355       finally show ?thesis unfolding u by (simp add: real_of_float_max)
  1356     qed
  1357     thus ?thesis unfolding l by auto
  1358   next case False note 4 = this show ?thesis
  1359   proof (cases "-2 * ?lpi \<le> ?lx \<and> ?ux \<le> 0")
  1360     case True note Cond = this with bnds 1 2 3 4
  1361     have l: "l = Float -1 0"
  1362       and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))"
  1363       by (auto simp add: bnds_cos_def Let_def)
  1364 
  1365     have "cos x \<le> real u"
  1366     proof (cases "-pi < x - real k * 2 * pi")
  1367       case True hence "-pi \<le> x - real k * 2 * pi" by simp
  1368       from negative_ux[OF this Cond[THEN conjunct2]]
  1369       show ?thesis unfolding u by (simp add: real_of_float_max)
  1370     next
  1371       case False hence "x - real k * 2 * pi \<le> -pi" by simp
  1372       hence pi_x: "x - real k * 2 * pi + 2 * pi \<le> pi" by simp
  1373 
  1374       have "-2 * pi \<le> real ?lx" using Cond lpi by (auto simp add: le_float_def)
  1375 
  1376       hence "0 \<le> x - real k * 2 * pi + 2 * pi" using lx by simp
  1377 
  1378       have lx_0: "0 \<le> real (?lx + 2 * ?lpi)"
  1379 	using Cond lpi by (auto simp add: le_float_def)
  1380 
  1381       from 1 and Cond have "\<not> -?lpi \<le> ?lx" by auto
  1382       hence "?lx + 2 * ?lpi \<le> ?lpi" by (auto simp add: le_float_def)
  1383       hence pi_lx: "real (?lx + 2 * ?lpi) \<le> pi"
  1384 	using lpi[THEN le_imp_neg_le] by (auto simp add: le_float_def)
  1385 
  1386       have lx_le_x: "real (?lx + 2 * ?lpi) \<le> x - real k * 2 * pi + 2 * pi"
  1387 	using lx lpi by auto
  1388 
  1389       have "cos x = cos (x + real (-k) * 2 * pi + real (1 :: int) * 2 * pi)"
  1390 	unfolding cos_periodic_int ..
  1391       also have "\<dots> \<le> cos (real (?lx + 2 * ?lpi))"
  1392 	using cos_monotone_0_pi'[OF lx_0 lx_le_x pi_x]
  1393 	by (simp only: real_of_float_minus real_of_int_minus real_of_one
  1394 	  number_of_Min real_diff_def mult_minus_left real_mult_1)
  1395       also have "\<dots> \<le> real (ub_cos prec (?lx + 2 * ?lpi))"
  1396 	using lb_cos[OF lx_0 pi_lx] by simp
  1397       finally show ?thesis unfolding u by (simp add: real_of_float_max)
  1398     qed
  1399     thus ?thesis unfolding l by auto
  1400   next
  1401     case False with bnds 1 2 3 4 show ?thesis by (auto simp add: bnds_cos_def Let_def)
  1402   qed qed qed qed qed
  1403 qed
  1404 
  1405 section "Exponential function"
  1406 
  1407 subsection "Compute the series of the exponential function"
  1408 
  1409 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
  1410 "ub_exp_horner prec 0 i k x       = 0" |
  1411 "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" |
  1412 "lb_exp_horner prec 0 i k x       = 0" |
  1413 "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"
  1414 
  1415 lemma bnds_exp_horner: assumes "real x \<le> 0"
  1416   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) }"
  1417 proof -
  1418   { fix n
  1419     have F: "\<And> m. ((\<lambda>i. i + 1) ^^ n) m = n + m" by (induct n, auto)
  1420     have "fact (Suc n) = fact n * ((\<lambda>i. i + 1) ^^ n) 1" unfolding F by auto } note f_eq = this
  1421 
  1422   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,
  1423     OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps]
  1424 
  1425   { 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)"
  1426       using bounds(1) by auto
  1427     also have "\<dots> \<le> exp (real x)"
  1428     proof -
  1429       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)"
  1430 	using Maclaurin_exp_le by blast
  1431       moreover have "0 \<le> exp t / real (fact (get_even n)) * (real x) ^ (get_even n)"
  1432 	by (auto intro!: mult_nonneg_nonneg divide_nonneg_pos simp add: get_even zero_le_even_power exp_gt_zero)
  1433       ultimately show ?thesis
  1434 	using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
  1435     qed
  1436     finally have "real (lb_exp_horner prec (get_even n) 1 1 x) \<le> exp (real x)" .
  1437   } moreover
  1438   { 
  1439     have x_less_zero: "real x ^ get_odd n \<le> 0"
  1440     proof (cases "real x = 0")
  1441       case True
  1442       have "(get_odd n) \<noteq> 0" using get_odd[THEN odd_pos] by auto
  1443       thus ?thesis unfolding True power_0_left by auto
  1444     next
  1445       case False hence "real x < 0" using `real x \<le> 0` by auto
  1446       show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq get_odd `real x < 0`)
  1447     qed
  1448 
  1449     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)"
  1450       using Maclaurin_exp_le by blast
  1451     moreover have "exp t / real (fact (get_odd n)) * (real x) ^ (get_odd n) \<le> 0"
  1452       by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero exp_gt_zero)
  1453     ultimately have "exp (real x) \<le> (\<Sum>j = 0..<get_odd n. 1 / real (fact j) * real x ^ j)"
  1454       using get_odd exp_gt_zero by (auto intro!: pordered_cancel_semiring_class.mult_nonneg_nonneg)
  1455     also have "\<dots> \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)"
  1456       using bounds(2) by auto
  1457     finally have "exp (real x) \<le> real (ub_exp_horner prec (get_odd n) 1 1 x)" .
  1458   } ultimately show ?thesis by auto
  1459 qed
  1460 
  1461 subsection "Compute the exponential function on the entire domain"
  1462 
  1463 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1464 "lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
  1465              else let 
  1466                 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)
  1467              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))
  1468                            else horner x)" |
  1469 "ub_exp prec x = (if 0 < x    then float_divr prec 1 (lb_exp prec (-x))
  1470              else if x < - 1  then (case floor_fl x of (Float m e) \<Rightarrow> 
  1471                                     (ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- Float m e))) ^ (nat (-m) * 2 ^ nat e))
  1472                               else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
  1473 by pat_completeness auto
  1474 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)
  1475 
  1476 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
  1477 proof -
  1478   have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
  1479 
  1480   have "1 / 4 = real (Float 1 -2)" unfolding Float_num by auto
  1481   also have "\<dots> \<le> real (lb_exp_horner 1 (get_even 4) 1 1 (- 1))"
  1482     unfolding get_even_def eq4 
  1483     by (auto simp add: lapprox_posrat_def rapprox_posrat_def normfloat.simps)
  1484   also have "\<dots> \<le> exp (real (- 1 :: float))" using bnds_exp_horner[where x="- 1"] by auto
  1485   finally show ?thesis unfolding real_of_float_minus real_of_float_1 . 
  1486 qed
  1487 
  1488 lemma lb_exp_pos: assumes "\<not> 0 < x" shows "0 < lb_exp prec x"
  1489 proof -
  1490   let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
  1491   let "?horner x" = "let  y = ?lb_horner x  in if y \<le> 0 then Float 1 -2 else y"
  1492   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)
  1493   moreover { fix x :: float fix num :: nat
  1494     have "0 < real (?horner x) ^ num" using `0 < ?horner x`[unfolded less_float_def real_of_float_0] by (rule zero_less_power)
  1495     also have "\<dots> = real ((?horner x) ^ num)" using float_power by auto
  1496     finally have "0 < real ((?horner x) ^ num)" .
  1497   }
  1498   ultimately show ?thesis
  1499     unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def
  1500     by (cases "floor_fl x", cases "x < - 1", auto simp add: float_power le_float_def less_float_def)
  1501 qed
  1502 
  1503 lemma exp_boundaries': assumes "x \<le> 0"
  1504   shows "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"
  1505 proof -
  1506   let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x"
  1507   let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x"
  1508 
  1509   have "real x \<le> 0" and "\<not> x > 0" using `x \<le> 0` unfolding le_float_def less_float_def by auto
  1510   show ?thesis
  1511   proof (cases "x < - 1")
  1512     case False hence "- 1 \<le> real x" unfolding less_float_def by auto
  1513     show ?thesis
  1514     proof (cases "?lb_exp_horner x \<le> 0")
  1515       from `\<not> x < - 1` have "- 1 \<le> real x" unfolding less_float_def by auto
  1516       hence "exp (- 1) \<le> exp (real x)" unfolding exp_le_cancel_iff .
  1517       from order_trans[OF exp_m1_ge_quarter this]
  1518       have "real (Float 1 -2) \<le> exp (real x)" unfolding Float_num .
  1519       moreover case True
  1520       ultimately show ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by auto
  1521     next
  1522       case False thus ?thesis using bnds_exp_horner `real x \<le> 0` `\<not> x > 0` `\<not> x < - 1` by (auto simp add: Let_def)
  1523     qed
  1524   next
  1525     case True
  1526     
  1527     obtain m e where Float_floor: "floor_fl x = Float m e" by (cases "floor_fl x", auto)
  1528     let ?num = "nat (- m) * 2 ^ nat e"
  1529     
  1530     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)
  1531     hence "real (floor_fl x) < 0" unfolding Float_floor real_of_float_simp using zero_less_pow2[of xe] by auto
  1532     hence "m < 0"
  1533       unfolding less_float_def real_of_float_0 Float_floor real_of_float_simp
  1534       unfolding pos_prod_lt[OF zero_less_pow2[of e], unfolded real_mult_commute] by auto
  1535     hence "1 \<le> - m" by auto
  1536     hence "0 < nat (- m)" by auto
  1537     moreover
  1538     have "0 \<le> e" using floor_pos_exp Float_floor[symmetric] by auto
  1539     hence "(0::nat) < 2 ^ nat e" by auto
  1540     ultimately have "0 < ?num"  by auto
  1541     hence "real ?num \<noteq> 0" by auto
  1542     have e_nat: "int (nat e) = e" using `0 \<le> e` by auto
  1543     have num_eq: "real ?num = real (- floor_fl x)" using `0 < nat (- m)`
  1544       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
  1545     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 .
  1546     hence "real (floor_fl x) < 0" unfolding less_float_def by auto
  1547     
  1548     have "exp (real x) \<le> real (ub_exp prec x)"
  1549     proof -
  1550       have div_less_zero: "real (float_divr prec x (- floor_fl x)) \<le> 0" 
  1551 	using float_divr_nonpos_pos_upper_bound[OF `x \<le> 0` `0 < - floor_fl x`] unfolding le_float_def real_of_float_0 .
  1552       
  1553       have "exp (real x) = exp (real ?num * (real x / real ?num))" using `real ?num \<noteq> 0` by auto
  1554       also have "\<dots> = exp (real x / real ?num) ^ ?num" unfolding exp_real_of_nat_mult ..
  1555       also have "\<dots> \<le> exp (real (float_divr prec x (- floor_fl x))) ^ ?num" unfolding num_eq
  1556 	by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto
  1557       also have "\<dots> \<le> real ((?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num)" unfolding float_power
  1558 	by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
  1559       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 .
  1560     qed
  1561     moreover 
  1562     have "real (lb_exp prec x) \<le> exp (real x)"
  1563     proof -
  1564       let ?divl = "float_divl prec x (- Float m e)"
  1565       let ?horner = "?lb_exp_horner ?divl"
  1566       
  1567       show ?thesis
  1568       proof (cases "?horner \<le> 0")
  1569 	case False hence "0 \<le> real ?horner" unfolding le_float_def by auto
  1570 	
  1571 	have div_less_zero: "real (float_divl prec x (- floor_fl x)) \<le> 0"
  1572 	  using `real (floor_fl x) < 0` `real x \<le> 0` by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg)
  1573 	
  1574 	have "real ((?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num) \<le>  
  1575           exp (real (float_divl prec x (- floor_fl x))) ^ ?num" unfolding float_power 
  1576 	  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)
  1577 	also have "\<dots> \<le> exp (real x / real ?num) ^ ?num" unfolding num_eq
  1578 	  using float_divl by (auto intro!: power_mono simp del: real_of_float_minus)
  1579 	also have "\<dots> = exp (real ?num * (real x / real ?num))" unfolding exp_real_of_nat_mult ..
  1580 	also have "\<dots> = exp (real x)" using `real ?num \<noteq> 0` by auto
  1581 	finally show ?thesis
  1582 	  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
  1583       next
  1584 	case True
  1585 	have "real (floor_fl x) \<noteq> 0" and "real (floor_fl x) \<le> 0" using `real (floor_fl x) < 0` by auto
  1586 	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`]]
  1587 	have "- 1 \<le> real x / real (- floor_fl x)" unfolding real_of_float_minus by auto
  1588 	from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]]
  1589 	have "real (Float 1 -2) \<le> exp (real x / real (- floor_fl x))" unfolding Float_num .
  1590 	hence "real (Float 1 -2) ^ ?num \<le> exp (real x / real (- floor_fl x)) ^ ?num"
  1591 	  by (auto intro!: power_mono simp add: Float_num)
  1592 	also have "\<dots> = exp (real x)" unfolding num_eq exp_real_of_nat_mult[symmetric] using `real (floor_fl x) \<noteq> 0` by auto
  1593 	finally show ?thesis
  1594 	  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 .
  1595       qed
  1596     qed
  1597     ultimately show ?thesis by auto
  1598   qed
  1599 qed
  1600 
  1601 lemma exp_boundaries: "exp (real x) \<in> { real (lb_exp prec x) .. real (ub_exp prec x)}"
  1602 proof -
  1603   show ?thesis
  1604   proof (cases "0 < x")
  1605     case False hence "x \<le> 0" unfolding less_float_def le_float_def by auto 
  1606     from exp_boundaries'[OF this] show ?thesis .
  1607   next
  1608     case True hence "-x \<le> 0" unfolding less_float_def le_float_def by auto
  1609     
  1610     have "real (lb_exp prec x) \<le> exp (real x)"
  1611     proof -
  1612       from exp_boundaries'[OF `-x \<le> 0`]
  1613       have ub_exp: "exp (- real x) \<le> real (ub_exp prec (-x))" unfolding atLeastAtMost_iff real_of_float_minus by auto
  1614       
  1615       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
  1616       also have "\<dots> \<le> exp (real x)"
  1617 	using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]]
  1618 	unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto
  1619       finally show ?thesis unfolding lb_exp.simps if_P[OF True] .
  1620     qed
  1621     moreover
  1622     have "exp (real x) \<le> real (ub_exp prec x)"
  1623     proof -
  1624       have "\<not> 0 < -x" using `0 < x` unfolding less_float_def by auto
  1625       
  1626       from exp_boundaries'[OF `-x \<le> 0`]
  1627       have lb_exp: "real (lb_exp prec (-x)) \<le> exp (- real x)" unfolding atLeastAtMost_iff real_of_float_minus by auto
  1628       
  1629       have "exp (real x) \<le> real (1 :: float) / real (lb_exp prec (-x))"
  1630 	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], 
  1631 	                                        symmetric]]
  1632 	unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide real_of_float_1 by auto
  1633       also have "\<dots> \<le> real (float_divr prec 1 (lb_exp prec (-x)))" using float_divr .
  1634       finally show ?thesis unfolding ub_exp.simps if_P[OF True] .
  1635     qed
  1636     ultimately show ?thesis by auto
  1637   qed
  1638 qed
  1639 
  1640 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"
  1641 proof (rule allI, rule allI, rule allI, rule impI)
  1642   fix x lx ux
  1643   assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \<and> x \<in> {real lx .. real ux}"
  1644   hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \<in> {real lx .. real ux}" by auto
  1645 
  1646   { from exp_boundaries[of lx prec, unfolded l]
  1647     have "real l \<le> exp (real lx)" by (auto simp del: lb_exp.simps)
  1648     also have "\<dots> \<le> exp x" using x by auto
  1649     finally have "real l \<le> exp x" .
  1650   } moreover
  1651   { have "exp x \<le> exp (real ux)" using x by auto
  1652     also have "\<dots> \<le> real u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps)
  1653     finally have "exp x \<le> real u" .
  1654   } ultimately show "real l \<le> exp x \<and> exp x \<le> real u" ..
  1655 qed
  1656 
  1657 section "Logarithm"
  1658 
  1659 subsection "Compute the logarithm series"
  1660 
  1661 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" 
  1662 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
  1663 "ub_ln_horner prec 0 i x       = 0" |
  1664 "ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
  1665 "lb_ln_horner prec 0 i x       = 0" |
  1666 "lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
  1667 
  1668 lemma ln_bounds:
  1669   assumes "0 \<le> x" and "x < 1"
  1670   shows "(\<Sum>i=0..<2*n. -1^i * (1 / real (i + 1)) * x ^ (Suc i)) \<le> ln (x + 1)" (is "?lb")
  1671   and "ln (x + 1) \<le> (\<Sum>i=0..<2*n + 1. -1^i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub")
  1672 proof -
  1673   let "?a n" = "(1/real (n +1)) * x ^ (Suc n)"
  1674 
  1675   have ln_eq: "(\<Sum> i. -1^i * ?a i) = ln (x + 1)"
  1676     using ln_series[of "x + 1"] `0 \<le> x` `x < 1` by auto
  1677 
  1678   have "norm x < 1" using assms by auto
  1679   have "?a ----> 0" unfolding Suc_eq_plus1[symmetric] inverse_eq_divide[symmetric] 
  1680     using LIMSEQ_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF `norm x < 1`]]] by auto
  1681   { fix n have "0 \<le> ?a n" by (rule mult_nonneg_nonneg, auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`) }
  1682   { fix n have "?a (Suc n) \<le> ?a n" unfolding inverse_eq_divide[symmetric]
  1683     proof (rule mult_mono)
  1684       show "0 \<le> x ^ Suc (Suc n)" by (auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
  1685       have "x ^ Suc (Suc n) \<le> x ^ Suc n * 1" unfolding power_Suc2 real_mult_assoc[symmetric] 
  1686 	by (rule mult_left_mono, fact less_imp_le[OF `x < 1`], auto intro!: mult_nonneg_nonneg simp add: `0 \<le> x`)
  1687       thus "x ^ Suc (Suc n) \<le> x ^ Suc n" by auto
  1688     qed auto }
  1689   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]
  1690   show "?lb" and "?ub" by auto
  1691 qed
  1692 
  1693 lemma ln_float_bounds: 
  1694   assumes "0 \<le> real x" and "real x < 1"
  1695   shows "real (x * lb_ln_horner prec (get_even n) 1 x) \<le> ln (real x + 1)" (is "?lb \<le> ?ln")
  1696   and "ln (real x + 1) \<le> real (x * ub_ln_horner prec (get_odd n) 1 x)" (is "?ln \<le> ?ub")
  1697 proof -
  1698   obtain ev where ev: "get_even n = 2 * ev" using get_even_double ..
  1699   obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double ..
  1700 
  1701   let "?s n" = "-1^n * (1 / real (1 + n)) * (real x)^(Suc n)"
  1702 
  1703   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
  1704     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",
  1705       OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`
  1706     by (rule mult_right_mono)
  1707   also have "\<dots> \<le> ?ln" using ln_bounds(1)[OF `0 \<le> real x` `real x < 1`] by auto
  1708   finally show "?lb \<le> ?ln" . 
  1709 
  1710   have "?ln \<le> setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF `0 \<le> real x` `real x < 1`] by auto
  1711   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
  1712     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",
  1713       OF `0 \<le> real x` refl lb_ln_horner.simps ub_ln_horner.simps] `0 \<le> real x`
  1714     by (rule mult_right_mono)
  1715   finally show "?ln \<le> ?ub" . 
  1716 qed
  1717 
  1718 lemma ln_add: assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)"
  1719 proof -
  1720   have "x \<noteq> 0" using assms by auto
  1721   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
  1722   moreover 
  1723   have "0 < y / x" using assms divide_pos_pos by auto
  1724   hence "0 < 1 + y / x" by auto
  1725   ultimately show ?thesis using ln_mult assms by auto
  1726 qed
  1727 
  1728 subsection "Compute the logarithm of 2"
  1729 
  1730 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3 
  1731                                         in (Float 1 -1 * ub_ln_horner prec (get_odd prec) 1 (Float 1 -1)) + 
  1732                                            (third * ub_ln_horner prec (get_odd prec) 1 third))"
  1733 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3 
  1734                                         in (Float 1 -1 * lb_ln_horner prec (get_even prec) 1 (Float 1 -1)) + 
  1735                                            (third * lb_ln_horner prec (get_even prec) 1 third))"
  1736 
  1737 lemma ub_ln2: "ln 2 \<le> real (ub_ln2 prec)" (is "?ub_ln2")
  1738   and lb_ln2: "real (lb_ln2 prec) \<le> ln 2" (is "?lb_ln2")
  1739 proof -
  1740   let ?uthird = "rapprox_rat (max prec 1) 1 3"
  1741   let ?lthird = "lapprox_rat prec 1 3"
  1742 
  1743   have ln2_sum: "ln 2 = ln (1/2 + 1) + ln (1 / 3 + 1)"
  1744     using ln_add[of "3 / 2" "1 / 2"] by auto
  1745   have lb3: "real ?lthird \<le> 1 / 3" using lapprox_rat[of prec 1 3] by auto
  1746   hence lb3_ub: "real ?lthird < 1" by auto
  1747   have lb3_lb: "0 \<le> real ?lthird" using lapprox_rat_bottom[of 1 3] by auto
  1748   have ub3: "1 / 3 \<le> real ?uthird" using rapprox_rat[of 1 3] by auto
  1749   hence ub3_lb: "0 \<le> real ?uthird" by auto
  1750 
  1751   have lb2: "0 \<le> real (Float 1 -1)" and ub2: "real (Float 1 -1) < 1" unfolding Float_num by auto
  1752 
  1753   have "0 \<le> (1::int)" and "0 < (3::int)" by auto
  1754   have ub3_ub: "real ?uthird < 1" unfolding rapprox_rat.simps(2)[OF `0 \<le> 1` `0 < 3`]
  1755     by (rule rapprox_posrat_less1, auto)
  1756 
  1757   have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
  1758   have uthird_gt0: "0 < real ?uthird + 1" using ub3_lb by auto
  1759   have lthird_gt0: "0 < real ?lthird + 1" using lb3_lb by auto
  1760 
  1761   show ?ub_ln2 unfolding ub_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
  1762   proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
  1763     have "ln (1 / 3 + 1) \<le> ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
  1764     also have "\<dots> \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)"
  1765       using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
  1766     finally show "ln (1 / 3 + 1) \<le> real (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .
  1767   qed
  1768   show ?lb_ln2 unfolding lb_ln2_def Let_def real_of_float_add ln2_sum Float_num(4)[symmetric]
  1769   proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
  1770     have "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (real ?lthird + 1)"
  1771       using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
  1772     also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
  1773     finally show "real (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .
  1774   qed
  1775 qed
  1776 
  1777 subsection "Compute the logarithm in the entire domain"
  1778 
  1779 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
  1780 "ub_ln prec x = (if x \<le> 0          then None
  1781             else if x < 1          then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
  1782             else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in
  1783                  if x \<le> Float 3 -1 then Some (horner (x - 1))
  1784             else if x < Float 1 1  then Some (horner (Float 1 -1) + horner (x * rapprox_rat prec 2 3 - 1))
  1785                                    else let l = bitlen (mantissa x) - 1 in
  1786                                         Some (ub_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))" |
  1787 "lb_ln prec x = (if x \<le> 0          then None
  1788             else if x < 1          then Some (- the (ub_ln prec (float_divr prec 1 x)))
  1789             else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in
  1790                  if x \<le> Float 3 -1 then Some (horner (x - 1))
  1791             else if x < Float 1 1  then Some (horner (Float 1 -1) +
  1792                                               horner (max (x * lapprox_rat prec 2 3 - 1) 0))
  1793                                    else let l = bitlen (mantissa x) - 1 in
  1794                                         Some (lb_ln2 prec * (Float (scale x + l) 0) + horner (Float (mantissa x) (- l) - 1)))"
  1795 by pat_completeness auto
  1796 
  1797 termination proof (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 1 then 1 else 0))", auto)
  1798   fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divl (max prec (Suc 0)) 1 x < 1"
  1799   hence "0 < x" and "0 < max prec (Suc 0)" unfolding less_float_def le_float_def by auto
  1800   from float_divl_pos_less1_bound[OF `0 < x` `x < 1` `0 < max prec (Suc 0)`]
  1801   show False using `float_divl (max prec (Suc 0)) 1 x < 1` unfolding less_float_def le_float_def by auto
  1802 next
  1803   fix prec x assume "\<not> x \<le> 0" and "x < 1" and "float_divr prec 1 x < 1"
  1804   hence "0 < x" unfolding less_float_def le_float_def by auto
  1805   from float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`, of prec]
  1806   show False using `float_divr prec 1 x < 1` unfolding less_float_def le_float_def by auto
  1807 qed
  1808 
  1809 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))))"
  1810 proof -
  1811   let ?B = "2^nat (bitlen m - 1)"
  1812   have "0 < real m" and "\<And>X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \<noteq> 0" using assms by auto
  1813   hence "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
  1814   show ?thesis
  1815   proof (cases "0 \<le> e")
  1816     case True
  1817     show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
  1818       unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
  1819       unfolding real_of_float_ge0_exp[OF True] ln_mult[OF `0 < real m` `0 < 2^nat e`]
  1820       ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` True by auto
  1821   next
  1822     case False hence "0 < -e" by auto
  1823     hence pow_gt0: "(0::real) < 2^nat (-e)" by auto
  1824     hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto
  1825     show ?thesis unfolding normalized_float[OF `m \<noteq> 0`]
  1826       unfolding ln_div[OF `0 < real m` `0 < ?B`] real_of_int_add ln_realpow[OF `0 < 2`]
  1827       unfolding real_of_float_nge0_exp[OF False] ln_mult[OF `0 < real m` inv_gt0] ln_inverse[OF pow_gt0]
  1828       ln_realpow[OF `0 < 2`] algebra_simps using `0 \<le> bitlen m - 1` False by auto
  1829   qed
  1830 qed
  1831 
  1832 lemma ub_ln_lb_ln_bounds': assumes "1 \<le> x"
  1833   shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"
  1834   (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
  1835 proof (cases "x < Float 1 1")
  1836   case True
  1837   hence "real (x - 1) < 1" and "real x < 2" unfolding less_float_def Float_num by auto
  1838   have "\<not> x \<le> 0" and "\<not> x < 1" using `1 \<le> x` unfolding less_float_def le_float_def by auto
  1839   hence "0 \<le> real (x - 1)" using `1 \<le> x` unfolding less_float_def Float_num by auto
  1840 
  1841   have [simp]: "real (Float 3 -1) = 3 / 2" by (simp add: real_of_float_def pow2_def)
  1842 
  1843   show ?thesis
  1844   proof (cases "x \<le> Float 3 -1")
  1845     case True
  1846     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
  1847       using ln_float_bounds[OF `0 \<le> real (x - 1)` `real (x - 1) < 1`, of prec] `\<not> x \<le> 0` `\<not> x < 1` True
  1848       by auto
  1849   next
  1850     case False hence *: "3 / 2 < real x" by (auto simp add: le_float_def)
  1851 
  1852     with ln_add[of "3 / 2" "real x - 3 / 2"]
  1853     have add: "ln (real x) = ln (3 / 2) + ln (real x * 2 / 3)"
  1854       by (auto simp add: algebra_simps diff_divide_distrib)
  1855 
  1856     let "?ub_horner x" = "x * ub_ln_horner prec (get_odd prec) 1 x"
  1857     let "?lb_horner x" = "x * lb_ln_horner prec (get_even prec) 1 x"
  1858 
  1859     { have up: "real (rapprox_rat prec 2 3) \<le> 1"
  1860 	by (rule rapprox_rat_le1) simp_all
  1861       have low: "2 / 3 \<le> real (rapprox_rat prec 2 3)"
  1862 	by (rule order_trans[OF _ rapprox_rat]) simp
  1863       from mult_less_le_imp_less[OF * low] *
  1864       have pos: "0 < real (x * rapprox_rat prec 2 3 - 1)" by auto
  1865 
  1866       have "ln (real x * 2/3)
  1867 	\<le> ln (real (x * rapprox_rat prec 2 3 - 1) + 1)"
  1868       proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
  1869 	show "real x * 2 / 3 \<le> real (x * rapprox_rat prec 2 3 - 1) + 1"
  1870 	  using * low by auto
  1871 	show "0 < real x * 2 / 3" using * by simp
  1872 	show "0 < real (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto
  1873       qed
  1874       also have "\<dots> \<le> real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"
  1875       proof (rule ln_float_bounds(2))
  1876 	from mult_less_le_imp_less[OF `real x < 2` up] low *
  1877 	show "real (x * rapprox_rat prec 2 3 - 1) < 1" by auto
  1878 	show "0 \<le> real (x * rapprox_rat prec 2 3 - 1)" using pos by auto
  1879       qed
  1880       finally have "ln (real x)
  1881 	\<le> real (?ub_horner (Float 1 -1))
  1882 	  + real (?ub_horner (x * rapprox_rat prec 2 3 - 1))"
  1883 	using ln_float_bounds(2)[of "Float 1 -1" prec prec] add by auto }
  1884     moreover
  1885     { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0"
  1886 
  1887       have up: "real (lapprox_rat prec 2 3) \<le> 2/3"
  1888 	by (rule order_trans[OF lapprox_rat], simp)
  1889 
  1890       have low: "0 \<le> real (lapprox_rat prec 2 3)"
  1891 	using lapprox_rat_bottom[of 2 3 prec] by simp
  1892 
  1893       have "real (?lb_horner ?max)
  1894 	\<le> ln (real ?max + 1)"
  1895       proof (rule ln_float_bounds(1))
  1896 	from mult_less_le_imp_less[OF `real x < 2` up] * low
  1897 	show "real ?max < 1" by (cases "real (lapprox_rat prec 2 3) = 0",
  1898 	  auto simp add: real_of_float_max)
  1899 	show "0 \<le> real ?max" by (auto simp add: real_of_float_max)
  1900       qed
  1901       also have "\<dots> \<le> ln (real x * 2/3)"
  1902       proof (rule ln_le_cancel_iff[symmetric, THEN iffD1])
  1903 	show "0 < real ?max + 1" by (auto simp add: real_of_float_max)
  1904 	show "0 < real x * 2/3" using * by auto
  1905 	show "real ?max + 1 \<le> real x * 2/3" using * up
  1906 	  by (cases "0 < real x * real (lapprox_posrat prec 2 3) - 1",
  1907 	      auto simp add: real_of_float_max max_def)
  1908       qed
  1909       finally have "real (?lb_horner (Float 1 -1)) + real (?lb_horner ?max)
  1910 	\<le> ln (real x)"
  1911 	using ln_float_bounds(1)[of "Float 1 -1" prec prec] add by auto }
  1912     ultimately
  1913     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
  1914       using `\<not> x \<le> 0` `\<not> x < 1` True False by auto
  1915   qed
  1916 next
  1917   case False
  1918   hence "\<not> x \<le> 0" and "\<not> x < 1" "0 < x" "\<not> x \<le> Float 3 -1"
  1919     using `1 \<le> x` unfolding less_float_def le_float_def real_of_float_simp pow2_def
  1920     by auto
  1921   show ?thesis
  1922   proof (cases x)
  1923     case (Float m e)
  1924     let ?s = "Float (e + (bitlen m - 1)) 0"
  1925     let ?x = "Float m (- (bitlen m - 1))"
  1926 
  1927     have "0 < m" and "m \<noteq> 0" using float_pos_m_pos `0 < x` Float by auto
  1928 
  1929     {
  1930       have "real (lb_ln2 prec * ?s) \<le> ln 2 * real (e + (bitlen m - 1))" (is "?lb2 \<le> _")
  1931 	unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
  1932 	using lb_ln2[of prec]
  1933       proof (rule mult_right_mono)
  1934 	have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
  1935 	from float_gt1_scale[OF this]
  1936 	show "0 \<le> real (e + (bitlen m - 1))" by auto
  1937       qed
  1938       moreover
  1939       from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
  1940       have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
  1941       from ln_float_bounds(1)[OF this]
  1942       have "real ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln (real ?x)" (is "?lb_horner \<le> _") by auto
  1943       ultimately have "?lb2 + ?lb_horner \<le> ln (real x)"
  1944 	unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
  1945     }
  1946     moreover
  1947     {
  1948       from bitlen_div[OF `0 < m`, unfolded normalized_float[OF `m \<noteq> 0`, symmetric]]
  1949       have "0 \<le> real (?x - 1)" and "real (?x - 1) < 1" by auto
  1950       from ln_float_bounds(2)[OF this]
  1951       have "ln (real ?x) \<le> real ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> ?ub_horner") by auto
  1952       moreover
  1953       have "ln 2 * real (e + (bitlen m - 1)) \<le> real (ub_ln2 prec * ?s)" (is "_ \<le> ?ub2")
  1954 	unfolding real_of_float_mult real_of_float_ge0_exp[OF order_refl] nat_0 power_0 mult_1_right
  1955 	using ub_ln2[of prec]
  1956       proof (rule mult_right_mono)
  1957 	have "1 \<le> Float m e" using `1 \<le> x` Float unfolding le_float_def by auto
  1958 	from float_gt1_scale[OF this]
  1959 	show "0 \<le> real (e + (bitlen m - 1))" by auto
  1960       qed
  1961       ultimately have "ln (real x) \<le> ?ub2 + ?ub_horner"
  1962 	unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
  1963     }
  1964     ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
  1965       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
  1966       unfolding scale.simps[of m e, unfolded Float[symmetric]] mantissa.simps[of m e, unfolded Float[symmetric]] real_of_float_add
  1967       by auto
  1968   qed
  1969 qed
  1970 
  1971 lemma ub_ln_lb_ln_bounds: assumes "0 < x"
  1972   shows "real (the (lb_ln prec x)) \<le> ln (real x) \<and> ln (real x) \<le> real (the (ub_ln prec x))"
  1973   (is "?lb \<le> ?ln \<and> ?ln \<le> ?ub")
  1974 proof (cases "x < 1")
  1975   case False hence "1 \<le> x" unfolding less_float_def le_float_def by auto
  1976   show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
  1977 next
  1978   case True have "\<not> x \<le> 0" using `0 < x` unfolding less_float_def le_float_def by auto
  1979 
  1980   have "0 < real x" and "real x \<noteq> 0" using `0 < x` unfolding less_float_def by auto
  1981   hence A: "0 < 1 / real x" by auto
  1982 
  1983   {
  1984     let ?divl = "float_divl (max prec 1) 1 x"
  1985     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
  1986     hence B: "0 < real ?divl" unfolding le_float_def by auto
  1987 
  1988     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
  1989     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
  1990     from this ub_ln_lb_ln_bounds'[OF A', THEN conjunct1, THEN le_imp_neg_le]
  1991     have "?ln \<le> real (- the (lb_ln prec ?divl))" unfolding real_of_float_minus by (rule order_trans)
  1992   } moreover
  1993   {
  1994     let ?divr = "float_divr prec 1 x"
  1995     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
  1996     hence B: "0 < real ?divr" unfolding le_float_def by auto
  1997 
  1998     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
  1999     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
  2000     from ub_ln_lb_ln_bounds'[OF A', THEN conjunct2, THEN le_imp_neg_le] this
  2001     have "real (- the (ub_ln prec ?divr)) \<le> ?ln" unfolding real_of_float_minus by (rule order_trans)
  2002   }
  2003   ultimately show ?thesis unfolding lb_ln.simps[where x=x]  ub_ln.simps[where x=x]
  2004     unfolding if_not_P[OF `\<not> x \<le> 0`] if_P[OF True] by auto
  2005 qed
  2006 
  2007 lemma lb_ln: assumes "Some y = lb_ln prec x"
  2008   shows "real y \<le> ln (real x)" and "0 < real x"
  2009 proof -
  2010   have "0 < x"
  2011   proof (rule ccontr)
  2012     assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
  2013     thus False using assms by auto
  2014   qed
  2015   thus "0 < real x" unfolding less_float_def by auto
  2016   have "real (the (lb_ln prec x)) \<le> ln (real x)" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
  2017   thus "real y \<le> ln (real x)" unfolding assms[symmetric] by auto
  2018 qed
  2019 
  2020 lemma ub_ln: assumes "Some y = ub_ln prec x"
  2021   shows "ln (real x) \<le> real y" and "0 < real x"
  2022 proof -
  2023   have "0 < x"
  2024   proof (rule ccontr)
  2025     assume "\<not> 0 < x" hence "x \<le> 0" unfolding le_float_def less_float_def by auto
  2026     thus False using assms by auto
  2027   qed
  2028   thus "0 < real x" unfolding less_float_def by auto
  2029   have "ln (real x) \<le> real (the (ub_ln prec x))" using ub_ln_lb_ln_bounds[OF `0 < x`] ..
  2030   thus "ln (real x) \<le> real y" unfolding assms[symmetric] by auto
  2031 qed
  2032 
  2033 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"
  2034 proof (rule allI, rule allI, rule allI, rule impI)
  2035   fix x lx ux
  2036   assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \<and> x \<in> {real lx .. real ux}"
  2037   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
  2038 
  2039   have "ln (real ux) \<le> real u" and "0 < real ux" using ub_ln u by auto
  2040   have "real l \<le> ln (real lx)" and "0 < real lx" and "0 < x" using lb_ln[OF l] x by auto
  2041 
  2042   from ln_le_cancel_iff[OF `0 < real lx` `0 < x`] `real l \<le> ln (real lx)`
  2043   have "real l \<le> ln x" using x unfolding atLeastAtMost_iff by auto
  2044   moreover
  2045   from ln_le_cancel_iff[OF `0 < x` `0 < real ux`] `ln (real ux) \<le> real u`
  2046   have "ln x \<le> real u" using x unfolding atLeastAtMost_iff by auto
  2047   ultimately show "real l \<le> ln x \<and> ln x \<le> real u" ..
  2048 qed
  2049 
  2050 section "Implement floatarith"
  2051 
  2052 subsection "Define syntax and semantics"
  2053 
  2054 datatype floatarith
  2055   = Add floatarith floatarith
  2056   | Minus floatarith
  2057   | Mult floatarith floatarith
  2058   | Inverse floatarith
  2059   | Cos floatarith
  2060   | Arctan floatarith
  2061   | Abs floatarith
  2062   | Max floatarith floatarith
  2063   | Min floatarith floatarith
  2064   | Pi
  2065   | Sqrt floatarith
  2066   | Exp floatarith
  2067   | Ln floatarith
  2068   | Power floatarith nat
  2069   | Atom nat
  2070   | Num float
  2071 
  2072 fun interpret_floatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real"
  2073 where
  2074 "interpret_floatarith (Add a b) vs   = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
  2075 "interpret_floatarith (Minus a) vs    = - (interpret_floatarith a vs)" |
  2076 "interpret_floatarith (Mult a b) vs   = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
  2077 "interpret_floatarith (Inverse a) vs  = inverse (interpret_floatarith a vs)" |
  2078 "interpret_floatarith (Cos a) vs      = cos (interpret_floatarith a vs)" |
  2079 "interpret_floatarith (Arctan a) vs   = arctan (interpret_floatarith a vs)" |
  2080 "interpret_floatarith (Min a b) vs    = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
  2081 "interpret_floatarith (Max a b) vs    = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
  2082 "interpret_floatarith (Abs a) vs      = abs (interpret_floatarith a vs)" |
  2083 "interpret_floatarith Pi vs           = pi" |
  2084 "interpret_floatarith (Sqrt a) vs     = sqrt (interpret_floatarith a vs)" |
  2085 "interpret_floatarith (Exp a) vs      = exp (interpret_floatarith a vs)" |
  2086 "interpret_floatarith (Ln a) vs       = ln (interpret_floatarith a vs)" |
  2087 "interpret_floatarith (Power a n) vs  = (interpret_floatarith a vs)^n" |
  2088 "interpret_floatarith (Num f) vs      = real f" |
  2089 "interpret_floatarith (Atom n) vs     = vs ! n"
  2090 
  2091 subsection "Implement approximation function"
  2092 
  2093 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
  2094 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
  2095 "lift_bin' a b f = None"
  2096 
  2097 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
  2098 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
  2099                                              | t \<Rightarrow> None)" |
  2100 "lift_un b f = None"
  2101 
  2102 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
  2103 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
  2104 "lift_un' b f = None"
  2105 
  2106 fun bounded_by :: "real list \<Rightarrow> (float * float) list \<Rightarrow> bool " where
  2107 bounded_by_Cons: "bounded_by (v#vs) ((l, u)#bs) = ((real l \<le> v \<and> v \<le> real u) \<and> bounded_by vs bs)" |
  2108 bounded_by_Nil: "bounded_by [] [] = True" |
  2109 "bounded_by _ _ = False"
  2110 
  2111 lemma bounded_by: assumes "bounded_by vs bs" and "i < length bs"
  2112   shows "real (fst (bs ! i)) \<le> vs ! i \<and> vs ! i \<le> real (snd (bs ! i))"
  2113   using `bounded_by vs bs` and `i < length bs`
  2114 proof (induct arbitrary: i rule: bounded_by.induct)
  2115   fix v :: real and vs :: "real list" and l u :: float and bs :: "(float * float) list" and i :: nat
  2116   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))"
  2117   assume bounded: "bounded_by (v # vs) ((l, u) # bs)" and length: "i < length ((l, u) # bs)"
  2118   show "real (fst (((l, u) # bs) ! i)) \<le> (v # vs) ! i \<and> (v # vs) ! i \<le> real (snd (((l, u) # bs) ! i))"
  2119   proof (cases i)
  2120     case 0
  2121     show ?thesis using bounded unfolding 0 nth_Cons_0 fst_conv snd_conv bounded_by.simps ..
  2122   next
  2123     case (Suc i) with length have "i < length bs" by auto
  2124     show ?thesis unfolding Suc nth_Cons_Suc bounded_by.simps
  2125       using hyp[OF bounded[unfolded bounded_by.simps, THEN conjunct2] `i < length bs`] .
  2126   qed
  2127 qed auto
  2128 
  2129 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) list \<Rightarrow> (float * float) option" where
  2130 "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)" |
  2131 "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))" | 
  2132 "approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
  2133 "approx prec (Mult a b) bs  = lift_bin' (approx' prec a bs) (approx' prec b bs)
  2134                                     (\<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, 
  2135                                                      float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
  2136 "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))" |
  2137 "approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |
  2138 "approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |
  2139 "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))" |
  2140 "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))" |
  2141 "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>))" |
  2142 "approx prec (Arctan a) bs  = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
  2143 "approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
  2144 "approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
  2145 "approx prec (Ln a) bs      = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
  2146 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
  2147 "approx prec (Num f) bs     = Some (f, f)" |
  2148 "approx prec (Atom i) bs    = (if i < length bs then Some (bs ! i) else None)"
  2149 
  2150 lemma lift_bin'_ex:
  2151   assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
  2152   shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
  2153 proof (cases a)
  2154   case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
  2155   thus ?thesis using lift_bin'_Some by auto
  2156 next
  2157   case (Some a')
  2158   show ?thesis
  2159   proof (cases b)
  2160     case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
  2161     thus ?thesis using lift_bin'_Some by auto
  2162   next
  2163     case (Some b')
  2164     obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2165     obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
  2166     thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
  2167   qed
  2168 qed
  2169 
  2170 lemma lift_bin'_f:
  2171   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
  2172   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"
  2173   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)"
  2174 proof -
  2175   obtain l1 u1 l2 u2
  2176     where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
  2177   have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto 
  2178   have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
  2179   thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto 
  2180 qed
  2181 
  2182 lemma approx_approx':
  2183   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"
  2184   and approx': "Some (l, u) = approx' prec a vs"
  2185   shows "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
  2186 proof -
  2187   obtain l' u' where S: "Some (l', u') = approx prec a vs"
  2188     using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
  2189   have l': "l = round_down prec l'" and u': "u = round_up prec u'"
  2190     using approx' unfolding approx'.simps S[symmetric] by auto
  2191   show ?thesis unfolding l' u' 
  2192     using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
  2193     using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
  2194 qed
  2195 
  2196 lemma lift_bin':
  2197   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
  2198   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")
  2199   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"
  2200   shows "\<exists> l1 u1 l2 u2. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 
  2201                         (real l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u2) \<and> 
  2202                         l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
  2203 proof -
  2204   { fix l u assume "Some (l, u) = approx' prec a bs"
  2205     with approx_approx'[of prec a bs, OF _ this] Pa
  2206     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2207   { fix l u assume "Some (l, u) = approx' prec b bs"
  2208     with approx_approx'[of prec b bs, OF _ this] Pb
  2209     have "real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u" by auto } note Pb = this
  2210 
  2211   from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
  2212   show ?thesis by auto
  2213 qed
  2214 
  2215 lemma lift_un'_ex:
  2216   assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
  2217   shows "\<exists> l u. Some (l, u) = a"
  2218 proof (cases a)
  2219   case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
  2220   thus ?thesis using lift_un'_Some by auto
  2221 next
  2222   case (Some a')
  2223   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2224   thus ?thesis unfolding `a = Some a'` a' by auto
  2225 qed
  2226 
  2227 lemma lift_un'_f:
  2228   assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
  2229   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
  2230   shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
  2231 proof -
  2232   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
  2233   have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
  2234   have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
  2235   thus ?thesis using Pa[OF Sa] by auto
  2236 qed
  2237 
  2238 lemma lift_un':
  2239   assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
  2240   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")
  2241   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 
  2242                         l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
  2243 proof -
  2244   { fix l u assume "Some (l, u) = approx' prec a bs"
  2245     with approx_approx'[of prec a bs, OF _ this] Pa
  2246     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2247   from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
  2248   show ?thesis by auto
  2249 qed
  2250 
  2251 lemma lift_un'_bnds:
  2252   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"
  2253   and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
  2254   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"
  2255   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
  2256 proof -
  2257   from lift_un'[OF lift_un'_Some Pa]
  2258   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
  2259   hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
  2260   thus ?thesis using bnds by auto
  2261 qed
  2262 
  2263 lemma lift_un_ex:
  2264   assumes lift_un_Some: "Some (l, u) = lift_un a f"
  2265   shows "\<exists> l u. Some (l, u) = a"
  2266 proof (cases a)
  2267   case None hence "None = lift_un a f" unfolding None lift_un.simps ..
  2268   thus ?thesis using lift_un_Some by auto
  2269 next
  2270   case (Some a')
  2271   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2272   thus ?thesis unfolding `a = Some a'` a' by auto
  2273 qed
  2274 
  2275 lemma lift_un_f:
  2276   assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
  2277   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
  2278   shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
  2279 proof -
  2280   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
  2281   have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
  2282   proof (rule ccontr)
  2283     assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
  2284     hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
  2285     hence "lift_un (g a) f = None" 
  2286     proof (cases "fst (f l1 u1) = None")
  2287       case True
  2288       then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
  2289       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
  2290     next
  2291       case False hence "snd (f l1 u1) = None" using or by auto
  2292       with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
  2293       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
  2294     qed
  2295     thus False using lift_un_Some by auto
  2296   qed
  2297   then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
  2298   from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
  2299   have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
  2300   thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
  2301 qed
  2302 
  2303 lemma lift_un:
  2304   assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
  2305   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")
  2306   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and> 
  2307                   Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
  2308 proof -
  2309   { fix l u assume "Some (l, u) = approx' prec a bs"
  2310     with approx_approx'[of prec a bs, OF _ this] Pa
  2311     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2312   from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
  2313   show ?thesis by auto
  2314 qed
  2315 
  2316 lemma lift_un_bnds:
  2317   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"
  2318   and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
  2319   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"
  2320   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
  2321 proof -
  2322   from lift_un[OF lift_un_Some Pa]
  2323   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
  2324   hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
  2325   thus ?thesis using bnds by auto
  2326 qed
  2327 
  2328 lemma approx:
  2329   assumes "bounded_by xs vs"
  2330   and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
  2331   shows "real l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> real u" (is "?P l u arith")
  2332   using `Some (l, u) = approx prec arith vs` 
  2333 proof (induct arith arbitrary: l u x)
  2334   case (Add a b)
  2335   from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
  2336   obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
  2337     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
  2338     "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
  2339   thus ?case unfolding interpret_floatarith.simps by auto
  2340 next
  2341   case (Minus a)
  2342   from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
  2343   obtain l1 u1 where "l = -u1" and "u = -l1"
  2344     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" unfolding fst_conv snd_conv by blast
  2345   thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto
  2346 next
  2347   case (Mult a b)
  2348   from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
  2349   obtain l1 u1 l2 u2 
  2350     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"
  2351     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"
  2352     and "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
  2353     and "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
  2354   thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt 
  2355     using mult_le_prts mult_ge_prts by auto
  2356 next
  2357   case (Inverse a)
  2358   from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
  2359   obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)" 
  2360     and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
  2361     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1" by blast
  2362   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
  2363   moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto
  2364   ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto
  2365 
  2366   have inv: "inverse (real u1) \<le> inverse (interpret_floatarith a xs)
  2367            \<and> inverse (interpret_floatarith a xs) \<le> inverse (real l1)"
  2368   proof (cases "0 < l1")
  2369     case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs" 
  2370       unfolding less_float_def using l1_le_u1 l1 by auto
  2371     show ?thesis
  2372       unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]
  2373 	inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]
  2374       using l1 u1 by auto
  2375   next
  2376     case False hence "u1 < 0" using either by blast
  2377     hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0" 
  2378       unfolding less_float_def using l1_le_u1 u1 by auto
  2379     show ?thesis
  2380       unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]
  2381 	inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]
  2382       using l1 u1 by auto
  2383   qed
  2384 
  2385   from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
  2386   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
  2387   also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto
  2388   finally have "real l \<le> inverse (interpret_floatarith a xs)" .
  2389   moreover
  2390   from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
  2391   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
  2392   hence "inverse (interpret_floatarith a xs) \<le> real u" by (rule order_trans[OF inv[THEN conjunct2]])
  2393   ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto
  2394 next
  2395   case (Abs x)
  2396   from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
  2397   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>"
  2398     and l1: "real l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> real u1" by blast
  2399   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)
  2400 next
  2401   case (Min a b)
  2402   from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
  2403   obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
  2404     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
  2405     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
  2406   thus ?case unfolding l' u' by (auto simp add: real_of_float_min)
  2407 next
  2408   case (Max a b)
  2409   from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
  2410   obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
  2411     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
  2412     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
  2413   thus ?case unfolding l' u' by (auto simp add: real_of_float_max)
  2414 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
  2415 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
  2416 next case Pi with pi_boundaries show ?case by auto
  2417 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
  2418 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
  2419 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
  2420 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
  2421 next case (Num f) thus ?case by auto
  2422 next
  2423   case (Atom n) 
  2424   show ?case
  2425   proof (cases "n < length vs")
  2426     case True
  2427     with Atom have "vs ! n = (l, u)" by auto
  2428     thus ?thesis using bounded_by[OF assms(1) True] by auto
  2429   next
  2430     case False thus ?thesis using Atom by auto
  2431   qed
  2432 qed
  2433 
  2434 datatype inequality = Less floatarith floatarith 
  2435                     | LessEqual floatarith floatarith 
  2436 
  2437 fun interpret_inequality :: "inequality \<Rightarrow> real list \<Rightarrow> bool" where 
  2438 "interpret_inequality (Less a b) vs                   = (interpret_floatarith a vs < interpret_floatarith b vs)" |
  2439 "interpret_inequality (LessEqual a b) vs              = (interpret_floatarith a vs \<le> interpret_floatarith b vs)"
  2440 
  2441 fun approx_inequality :: "nat \<Rightarrow> inequality \<Rightarrow> (float * float) list \<Rightarrow> bool" where 
  2442 "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)" |
  2443 "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)"
  2444 
  2445 lemma approx_inequality: fixes m :: nat assumes "bounded_by vs bs" and "approx_inequality prec eq bs"
  2446   shows "interpret_inequality eq vs"
  2447 proof (cases eq)
  2448   case (Less a b)
  2449   show ?thesis
  2450   proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and> 
  2451                              approx prec b bs = Some (l', u')")
  2452     case True
  2453     then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
  2454       and b_approx: "approx prec b bs = Some (l', u') " by auto
  2455     with `approx_inequality prec eq bs` have "real u < real l'"
  2456       unfolding Less approx_inequality.simps less_float_def by auto
  2457     moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
  2458     have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"
  2459       using approx by auto
  2460     ultimately show ?thesis unfolding interpret_inequality.simps Less by auto
  2461   next
  2462     case False
  2463     hence "approx prec a bs = None \<or> approx prec b bs = None"
  2464       unfolding not_Some_eq[symmetric] by auto
  2465     hence "\<not> approx_inequality prec eq bs" unfolding Less approx_inequality.simps 
  2466       by (cases "approx prec a bs = None", auto)
  2467     thus ?thesis using assms by auto
  2468   qed
  2469 next
  2470   case (LessEqual a b)
  2471   show ?thesis
  2472   proof (cases "\<exists> u l u' l'. approx prec a bs = Some (l, u) \<and> 
  2473                              approx prec b bs = Some (l', u')")
  2474     case True
  2475     then obtain l u l' u' where a_approx: "approx prec a bs = Some (l, u)"
  2476       and b_approx: "approx prec b bs = Some (l', u') " by auto
  2477     with `approx_inequality prec eq bs` have "real u \<le> real l'"
  2478       unfolding LessEqual approx_inequality.simps le_float_def by auto
  2479     moreover from a_approx[symmetric] and b_approx[symmetric] and `bounded_by vs bs`
  2480     have "interpret_floatarith a vs \<le> real u" and "real l' \<le> interpret_floatarith b vs"
  2481       using approx by auto
  2482     ultimately show ?thesis unfolding interpret_inequality.simps LessEqual by auto
  2483   next
  2484     case False
  2485     hence "approx prec a bs = None \<or> approx prec b bs = None"
  2486       unfolding not_Some_eq[symmetric] by auto
  2487     hence "\<not> approx_inequality prec eq bs" unfolding LessEqual approx_inequality.simps 
  2488       by (cases "approx prec a bs = None", auto)
  2489     thus ?thesis using assms by auto
  2490   qed
  2491 qed
  2492 
  2493 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"
  2494   unfolding real_divide_def interpret_floatarith.simps ..
  2495 
  2496 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"
  2497   unfolding real_diff_def interpret_floatarith.simps ..
  2498 
  2499 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs = 
  2500   sin (interpret_floatarith a vs)"
  2501   unfolding sin_cos_eq interpret_floatarith.simps
  2502             interpret_floatarith_divide interpret_floatarith_diff real_diff_def
  2503   by auto
  2504 
  2505 lemma interpret_floatarith_tan:
  2506   "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =
  2507    tan (interpret_floatarith a vs)"
  2508   unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def real_divide_def
  2509   by auto
  2510 
  2511 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"
  2512   unfolding powr_def interpret_floatarith.simps ..
  2513 
  2514 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"
  2515   unfolding log_def interpret_floatarith.simps real_divide_def ..
  2516 
  2517 lemma interpret_floatarith_num:
  2518   shows "interpret_floatarith (Num (Float 0 0)) vs = 0"
  2519   and "interpret_floatarith (Num (Float 1 0)) vs = 1"
  2520   and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
  2521 
  2522 subsection {* Implement proof method \texttt{approximation} *}
  2523 
  2524 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)
  2525 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)
  2526 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)"
  2527                      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)"
  2528                      and "real (Float (number_of A) (int B)) = (number_of A) * 2^B"
  2529                      and "real (Float 1 (int B)) = 2^B"
  2530                      and "real (Float (number_of A) (- int B)) = (number_of A) / 2^B"
  2531                      and "real (Float 1 (- int B)) = 1 / 2^B"
  2532   by (auto simp add: real_of_float_simp pow2_def real_divide_def)
  2533 
  2534 lemmas bounded_by_equations = bounded_by_Cons bounded_by_Nil float_power bounded_divl bounded_divr bounded_num HOL.simp_thms
  2535 lemmas interpret_inequality_equations = interpret_inequality.simps interpret_floatarith.simps interpret_floatarith_num
  2536   interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log
  2537   interpret_floatarith_sin
  2538 
  2539 ML {*
  2540 structure Float_Arith =
  2541 struct
  2542 
  2543 @{code_datatype float = Float}
  2544 @{code_datatype floatarith = Add | Minus | Mult | Inverse | Cos | Arctan
  2545                            | Abs | Max | Min | Pi | Sqrt | Exp | Ln | Power | Atom | Num }
  2546 @{code_datatype inequality = Less | LessEqual }
  2547 
  2548 val approx_inequality = @{code approx_inequality}
  2549 
  2550 end
  2551 *}
  2552 
  2553 code_reserved Eval Float_Arith
  2554 
  2555 code_type float (Eval "Float'_Arith.float")
  2556 code_const Float (Eval "Float'_Arith.Float/ (_,/ _)")
  2557 
  2558 code_type floatarith (Eval "Float'_Arith.floatarith")
  2559 code_const Add and Minus and Mult and Inverse and Cos and Arctan and Abs and Max and Min and
  2560            Pi and Sqrt  and Exp and Ln and Power and Atom and Num
  2561   (Eval "Float'_Arith.Add/ (_,/ _)" and "Float'_Arith.Minus" and "Float'_Arith.Mult/ (_,/ _)" and
  2562         "Float'_Arith.Inverse" and "Float'_Arith.Cos" and
  2563         "Float'_Arith.Arctan" and "Float'_Arith.Abs" and "Float'_Arith.Max/ (_,/ _)" and
  2564         "Float'_Arith.Min/ (_,/ _)" and "Float'_Arith.Pi" and "Float'_Arith.Sqrt" and
  2565         "Float'_Arith.Exp" and "Float'_Arith.Ln" and "Float'_Arith.Power/ (_,/ _)" and
  2566         "Float'_Arith.Atom" and "Float'_Arith.Num")
  2567 
  2568 code_type inequality (Eval "Float'_Arith.inequality")
  2569 code_const Less and LessEqual (Eval "Float'_Arith.Less/ (_,/ _)" and "Float'_Arith.LessEqual/ (_,/ _)")
  2570 
  2571 code_const approx_inequality (Eval "Float'_Arith.approx'_inequality")
  2572 
  2573 ML {*
  2574   val ineq_equations = PureThy.get_thms @{theory} "interpret_inequality_equations";
  2575   val bounded_by_equations = PureThy.get_thms @{theory} "bounded_by_equations";
  2576   val bounded_by_simpset = (HOL_basic_ss addsimps bounded_by_equations)
  2577 
  2578   fun reify_ineq ctxt i = (fn st =>
  2579     let
  2580       val to = HOLogic.dest_Trueprop (Logic.strip_imp_concl (List.nth (prems_of st, i - 1)))
  2581     in (Reflection.genreify_tac ctxt ineq_equations (SOME to) i) st
  2582     end)
  2583 
  2584   fun rule_ineq ctxt prec i thm = let
  2585     fun conv_num typ = HOLogic.dest_number #> snd #> HOLogic.mk_number typ
  2586     val to_natc = conv_num @{typ "nat"} #> Thm.cterm_of (ProofContext.theory_of ctxt)
  2587     val to_nat = conv_num @{typ "nat"}
  2588     val to_int = conv_num @{typ "int"}
  2589     fun int_to_float A = @{term "Float"} $ to_int A $ @{term "0 :: int"}
  2590 
  2591     val prec' = to_nat prec
  2592 
  2593     fun bot_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
  2594                    = @{term "Float"} $ to_int mantisse $ to_int exp
  2595       | bot_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
  2596                    = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
  2597       | bot_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
  2598                    = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
  2599       | bot_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
  2600                    = @{term "float_divl"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
  2601       | bot_float (Const (@{const_name "divide"}, _) $ A $ B)
  2602                    = @{term "float_divl"} $ prec' $ int_to_float A $ int_to_float B
  2603       | bot_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
  2604                    = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
  2605       | bot_float A = int_to_float A
  2606 
  2607     fun top_float (Const (@{const_name "times"}, _) $ mantisse $ (Const (@{const_name "pow2"}, _) $ exp))
  2608                    = @{term "Float"} $ to_int mantisse $ to_int exp
  2609       | top_float (Const (@{const_name "divide"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
  2610                    = @{term "Float"} $ to_int mantisse $ (@{term "uminus :: int \<Rightarrow> int"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp))
  2611       | top_float (Const (@{const_name "times"}, _) $ mantisse $ (@{term "power 2 :: nat \<Rightarrow> real"} $ exp))
  2612                    = @{term "Float"} $ to_int mantisse $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
  2613       | top_float (Const (@{const_name "divide"}, _) $ A $ (@{term "power 10 :: nat \<Rightarrow> real"} $ exp))
  2614                    = @{term "float_divr"} $ prec' $ (int_to_float A) $ (@{term "power (Float 5 1)"} $ to_nat exp)
  2615       | top_float (Const (@{const_name "divide"}, _) $ A $ B)
  2616                    = @{term "float_divr"} $ prec' $ int_to_float A $ int_to_float B
  2617       | top_float (@{term "power 2 :: nat \<Rightarrow> real"} $ exp)
  2618                    = @{term "Float 1"} $ (@{term "int :: nat \<Rightarrow> int"} $ to_nat exp)
  2619       | top_float A = int_to_float A
  2620 
  2621     val goal' : term = List.nth (prems_of thm, i - 1)
  2622 
  2623     fun lift_bnd (t as (Const (@{const_name "op &"}, _) $
  2624                         (Const (@{const_name "less_eq"}, _) $
  2625                          bottom $ (Free (name, _))) $
  2626                         (Const (@{const_name "less_eq"}, _) $ _ $ top)))
  2627          = ((name, HOLogic.mk_prod (bot_float bottom, top_float top))
  2628             handle TERM (txt, ts) => raise TERM ("Invalid bound number format: " ^
  2629                                   (Syntax.string_of_term ctxt t), [t]))
  2630       | lift_bnd t = raise TERM ("Premisse needs format '<num> <= <var> & <var> <= <num>', but found " ^
  2631                                  (Syntax.string_of_term ctxt t), [t])
  2632     val bound_eqs = map (HOLogic.dest_Trueprop #> lift_bnd)  (Logic.strip_imp_prems goal')
  2633 
  2634     fun lift_var (Free (varname, _)) = (case AList.lookup (op =) bound_eqs varname of
  2635                                           SOME bound => bound
  2636                                         | NONE => raise TERM ("No bound equations found for " ^ varname, []))
  2637       | lift_var t = raise TERM ("Can not convert expression " ^
  2638                                  (Syntax.string_of_term ctxt t), [t])
  2639 
  2640     val _ $ vs = HOLogic.dest_Trueprop (Logic.strip_imp_concl goal')
  2641 
  2642     val bs = (HOLogic.dest_list #> map lift_var #> HOLogic.mk_list @{typ "float * float"}) vs
  2643     val map = [(@{cpat "?prec::nat"}, to_natc prec),
  2644                (@{cpat "?bs::(float * float) list"}, Thm.cterm_of (ProofContext.theory_of ctxt) bs)]
  2645   in rtac (Thm.instantiate ([], map) @{thm "approx_inequality"}) i thm end
  2646 
  2647   val eval_tac = CSUBGOAL (fn (ct, i) => rtac (eval_oracle ct) i)
  2648 
  2649   fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
  2650                                THEN' rtac TrueI
  2651 
  2652 *}
  2653 
  2654 method_setup approximation = {*
  2655   Args.term >>
  2656   (fn prec => fn ctxt =>
  2657     SIMPLE_METHOD' (fn i =>
  2658      (DETERM (reify_ineq ctxt i)
  2659       THEN rule_ineq ctxt prec i
  2660       THEN Simplifier.asm_full_simp_tac bounded_by_simpset i
  2661       THEN (TRY (filter_prems_tac (fn t => false) i))
  2662       THEN (gen_eval_tac eval_oracle ctxt) i)))
  2663 *} "real number approximation"
  2664 
  2665 end