src/HOL/Decision_Procs/Approximation.thy
author haftmann
Mon Apr 26 15:37:50 2010 +0200 (2010-04-26)
changeset 36409 d323e7773aa8
parent 35845 e5980f0ad025
child 36526 353041483b9b
permissions -rw-r--r--
use new classes (linordered_)field_inverse_zero
     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!: 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!: 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] real_of_nat_power 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 min_max.sup_absorb1)
  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   | Var nat
  2070   | Num float
  2071 
  2072 fun interpret_floatarith :: "floatarith \<Rightarrow> real list \<Rightarrow> real" where
  2073 "interpret_floatarith (Add a b) vs   = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
  2074 "interpret_floatarith (Minus a) vs    = - (interpret_floatarith a vs)" |
  2075 "interpret_floatarith (Mult a b) vs   = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
  2076 "interpret_floatarith (Inverse a) vs  = inverse (interpret_floatarith a vs)" |
  2077 "interpret_floatarith (Cos a) vs      = cos (interpret_floatarith a vs)" |
  2078 "interpret_floatarith (Arctan a) vs   = arctan (interpret_floatarith a vs)" |
  2079 "interpret_floatarith (Min a b) vs    = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
  2080 "interpret_floatarith (Max a b) vs    = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
  2081 "interpret_floatarith (Abs a) vs      = abs (interpret_floatarith a vs)" |
  2082 "interpret_floatarith Pi vs           = pi" |
  2083 "interpret_floatarith (Sqrt a) vs     = sqrt (interpret_floatarith a vs)" |
  2084 "interpret_floatarith (Exp a) vs      = exp (interpret_floatarith a vs)" |
  2085 "interpret_floatarith (Ln a) vs       = ln (interpret_floatarith a vs)" |
  2086 "interpret_floatarith (Power a n) vs  = (interpret_floatarith a vs)^n" |
  2087 "interpret_floatarith (Num f) vs      = real f" |
  2088 "interpret_floatarith (Var n) vs     = vs ! n"
  2089 
  2090 lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)"
  2091   unfolding real_divide_def interpret_floatarith.simps ..
  2092 
  2093 lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)"
  2094   unfolding real_diff_def interpret_floatarith.simps ..
  2095 
  2096 lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) vs =
  2097   sin (interpret_floatarith a vs)"
  2098   unfolding sin_cos_eq interpret_floatarith.simps
  2099             interpret_floatarith_divide interpret_floatarith_diff real_diff_def
  2100   by auto
  2101 
  2102 lemma interpret_floatarith_tan:
  2103   "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (Inverse (Cos a))) vs =
  2104    tan (interpret_floatarith a vs)"
  2105   unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def real_divide_def
  2106   by auto
  2107 
  2108 lemma interpret_floatarith_powr: "interpret_floatarith (Exp (Mult b (Ln a))) vs = (interpret_floatarith a vs) powr (interpret_floatarith b vs)"
  2109   unfolding powr_def interpret_floatarith.simps ..
  2110 
  2111 lemma interpret_floatarith_log: "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = log (interpret_floatarith b vs) (interpret_floatarith x vs)"
  2112   unfolding log_def interpret_floatarith.simps real_divide_def ..
  2113 
  2114 lemma interpret_floatarith_num:
  2115   shows "interpret_floatarith (Num (Float 0 0)) vs = 0"
  2116   and "interpret_floatarith (Num (Float 1 0)) vs = 1"
  2117   and "interpret_floatarith (Num (Float (number_of a) 0)) vs = number_of a" by auto
  2118 
  2119 subsection "Implement approximation function"
  2120 
  2121 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
  2122 "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
  2123 "lift_bin' a b f = None"
  2124 
  2125 fun lift_un :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> ((float option) * (float option))) \<Rightarrow> (float * float) option" where
  2126 "lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) \<Rightarrow> Some (l, u)
  2127                                              | t \<Rightarrow> None)" |
  2128 "lift_un b f = None"
  2129 
  2130 fun lift_un' :: "(float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
  2131 "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
  2132 "lift_un' b f = None"
  2133 
  2134 definition
  2135 "bounded_by xs vs \<longleftrightarrow>
  2136   (\<forall> i < length vs. case vs ! i of None \<Rightarrow> True
  2137          | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u })"
  2138 
  2139 lemma bounded_byE:
  2140   assumes "bounded_by xs vs"
  2141   shows "\<And> i. i < length vs \<Longrightarrow> case vs ! i of None \<Rightarrow> True
  2142          | Some (l, u) \<Rightarrow> xs ! i \<in> { real l .. real u }"
  2143   using assms bounded_by_def by blast
  2144 
  2145 lemma bounded_by_update:
  2146   assumes "bounded_by xs vs"
  2147   and bnd: "xs ! i \<in> { real l .. real u }"
  2148   shows "bounded_by xs (vs[i := Some (l,u)])"
  2149 proof -
  2150 { fix j
  2151   let ?vs = "vs[i := Some (l,u)]"
  2152   assume "j < length ?vs" hence [simp]: "j < length vs" by simp
  2153   have "case ?vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> xs ! j \<in> { real l .. real u }"
  2154   proof (cases "?vs ! j")
  2155     case (Some b)
  2156     thus ?thesis
  2157     proof (cases "i = j")
  2158       case True
  2159       thus ?thesis using `?vs ! j = Some b` and bnd by auto
  2160     next
  2161       case False
  2162       thus ?thesis using `bounded_by xs vs` unfolding bounded_by_def by auto
  2163     qed
  2164   qed auto }
  2165   thus ?thesis unfolding bounded_by_def by auto
  2166 qed
  2167 
  2168 lemma bounded_by_None:
  2169   shows "bounded_by xs (replicate (length xs) None)"
  2170   unfolding bounded_by_def by auto
  2171 
  2172 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option" where
  2173 "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)" |
  2174 "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))" |
  2175 "approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
  2176 "approx prec (Mult a b) bs  = lift_bin' (approx' prec a bs) (approx' prec b bs)
  2177                                     (\<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,
  2178                                                      float_pprt a2 * float_pprt b2 + float_pprt a1 * float_nprt b2 + float_nprt a2 * float_pprt b1 + float_nprt a1 * float_nprt b1))" |
  2179 "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))" |
  2180 "approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |
  2181 "approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |
  2182 "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))" |
  2183 "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))" |
  2184 "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>))" |
  2185 "approx prec (Arctan a) bs  = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
  2186 "approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
  2187 "approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
  2188 "approx prec (Ln a) bs      = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
  2189 "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
  2190 "approx prec (Num f) bs     = Some (f, f)" |
  2191 "approx prec (Var i) bs    = (if i < length bs then bs ! i else None)"
  2192 
  2193 lemma lift_bin'_ex:
  2194   assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
  2195   shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
  2196 proof (cases a)
  2197   case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
  2198   thus ?thesis using lift_bin'_Some by auto
  2199 next
  2200   case (Some a')
  2201   show ?thesis
  2202   proof (cases b)
  2203     case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps ..
  2204     thus ?thesis using lift_bin'_Some by auto
  2205   next
  2206     case (Some b')
  2207     obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2208     obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto)
  2209     thus ?thesis unfolding `a = Some a'` `b = Some b'` a' b' by auto
  2210   qed
  2211 qed
  2212 
  2213 lemma lift_bin'_f:
  2214   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
  2215   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"
  2216   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)"
  2217 proof -
  2218   obtain l1 u1 l2 u2
  2219     where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto
  2220   have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
  2221   have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto
  2222   thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto
  2223 qed
  2224 
  2225 lemma approx_approx':
  2226   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"
  2227   and approx': "Some (l, u) = approx' prec a vs"
  2228   shows "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u"
  2229 proof -
  2230   obtain l' u' where S: "Some (l', u') = approx prec a vs"
  2231     using approx' unfolding approx'.simps by (cases "approx prec a vs", auto)
  2232   have l': "l = round_down prec l'" and u': "u = round_up prec u'"
  2233     using approx' unfolding approx'.simps S[symmetric] by auto
  2234   show ?thesis unfolding l' u'
  2235     using order_trans[OF Pa[OF S, THEN conjunct2] round_up[of u']]
  2236     using order_trans[OF round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
  2237 qed
  2238 
  2239 lemma lift_bin':
  2240   assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
  2241   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")
  2242   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"
  2243   shows "\<exists> l1 u1 l2 u2. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
  2244                         (real l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u2) \<and>
  2245                         l = fst (f l1 u1 l2 u2) \<and> u = snd (f l1 u1 l2 u2)"
  2246 proof -
  2247   { fix l u assume "Some (l, u) = approx' prec a bs"
  2248     with approx_approx'[of prec a bs, OF _ this] Pa
  2249     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2250   { fix l u assume "Some (l, u) = approx' prec b bs"
  2251     with approx_approx'[of prec b bs, OF _ this] Pb
  2252     have "real l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real u" by auto } note Pb = this
  2253 
  2254   from lift_bin'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
  2255   show ?thesis by auto
  2256 qed
  2257 
  2258 lemma lift_un'_ex:
  2259   assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
  2260   shows "\<exists> l u. Some (l, u) = a"
  2261 proof (cases a)
  2262   case None hence "None = lift_un' a f" unfolding None lift_un'.simps ..
  2263   thus ?thesis using lift_un'_Some by auto
  2264 next
  2265   case (Some a')
  2266   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2267   thus ?thesis unfolding `a = Some a'` a' by auto
  2268 qed
  2269 
  2270 lemma lift_un'_f:
  2271   assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
  2272   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
  2273   shows "\<exists> l1 u1. P l1 u1 a \<and> l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
  2274 proof -
  2275   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto
  2276   have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
  2277   have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto
  2278   thus ?thesis using Pa[OF Sa] by auto
  2279 qed
  2280 
  2281 lemma lift_un':
  2282   assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
  2283   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")
  2284   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
  2285                         l = fst (f l1 u1) \<and> u = snd (f l1 u1)"
  2286 proof -
  2287   { fix l u assume "Some (l, u) = approx' prec a bs"
  2288     with approx_approx'[of prec a bs, OF _ this] Pa
  2289     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2290   from lift_un'_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
  2291   show ?thesis by auto
  2292 qed
  2293 
  2294 lemma lift_un'_bnds:
  2295   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"
  2296   and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
  2297   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"
  2298   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
  2299 proof -
  2300   from lift_un'[OF lift_un'_Some Pa]
  2301   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
  2302   hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
  2303   thus ?thesis using bnds by auto
  2304 qed
  2305 
  2306 lemma lift_un_ex:
  2307   assumes lift_un_Some: "Some (l, u) = lift_un a f"
  2308   shows "\<exists> l u. Some (l, u) = a"
  2309 proof (cases a)
  2310   case None hence "None = lift_un a f" unfolding None lift_un.simps ..
  2311   thus ?thesis using lift_un_Some by auto
  2312 next
  2313   case (Some a')
  2314   obtain la ua where a': "a' = (la, ua)" by (cases a', auto)
  2315   thus ?thesis unfolding `a = Some a'` a' by auto
  2316 qed
  2317 
  2318 lemma lift_un_f:
  2319   assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
  2320   and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
  2321   shows "\<exists> l1 u1. P l1 u1 a \<and> Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
  2322 proof -
  2323   obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto
  2324   have "fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None"
  2325   proof (rule ccontr)
  2326     assume "\<not> (fst (f l1 u1) \<noteq> None \<and> snd (f l1 u1) \<noteq> None)"
  2327     hence or: "fst (f l1 u1) = None \<or> snd (f l1 u1) = None" by auto
  2328     hence "lift_un (g a) f = None"
  2329     proof (cases "fst (f l1 u1) = None")
  2330       case True
  2331       then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto)
  2332       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
  2333     next
  2334       case False hence "snd (f l1 u1) = None" using or by auto
  2335       with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto)
  2336       thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto
  2337     qed
  2338     thus False using lift_un_Some by auto
  2339   qed
  2340   then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto)
  2341   from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
  2342   have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto
  2343   thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
  2344 qed
  2345 
  2346 lemma lift_un:
  2347   assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
  2348   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")
  2349   shows "\<exists> l1 u1. (real l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u1) \<and>
  2350                   Some l = fst (f l1 u1) \<and> Some u = snd (f l1 u1)"
  2351 proof -
  2352   { fix l u assume "Some (l, u) = approx' prec a bs"
  2353     with approx_approx'[of prec a bs, OF _ this] Pa
  2354     have "real l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real u" by auto } note Pa = this
  2355   from lift_un_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
  2356   show ?thesis by auto
  2357 qed
  2358 
  2359 lemma lift_un_bnds:
  2360   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"
  2361   and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
  2362   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"
  2363   shows "real l \<le> f' (interpret_floatarith a xs) \<and> f' (interpret_floatarith a xs) \<le> real u"
  2364 proof -
  2365   from lift_un[OF lift_un_Some Pa]
  2366   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
  2367   hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \<in> {real l1 .. real u1}" by auto
  2368   thus ?thesis using bnds by auto
  2369 qed
  2370 
  2371 lemma approx:
  2372   assumes "bounded_by xs vs"
  2373   and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
  2374   shows "real l \<le> interpret_floatarith arith xs \<and> interpret_floatarith arith xs \<le> real u" (is "?P l u arith")
  2375   using `Some (l, u) = approx prec arith vs`
  2376 proof (induct arith arbitrary: l u x)
  2377   case (Add a b)
  2378   from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
  2379   obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
  2380     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
  2381     "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
  2382   thus ?case unfolding interpret_floatarith.simps by auto
  2383 next
  2384   case (Minus a)
  2385   from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
  2386   obtain l1 u1 where "l = -u1" and "u = -l1"
  2387     "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1" unfolding fst_conv snd_conv by blast
  2388   thus ?case unfolding interpret_floatarith.simps using real_of_float_minus by auto
  2389 next
  2390   case (Mult a b)
  2391   from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
  2392   obtain l1 u1 l2 u2
  2393     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"
  2394     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"
  2395     and "real l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> real u1"
  2396     and "real l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> real u2" unfolding fst_conv snd_conv by blast
  2397   thus ?case unfolding interpret_floatarith.simps l u real_of_float_add real_of_float_mult real_of_float_nprt real_of_float_pprt
  2398     using mult_le_prts mult_ge_prts by auto
  2399 next
  2400   case (Inverse a)
  2401   from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
  2402   obtain l1 u1 where l': "Some l = (if 0 < l1 \<or> u1 < 0 then Some (float_divl prec 1 u1) else None)"
  2403     and u': "Some u = (if 0 < l1 \<or> u1 < 0 then Some (float_divr prec 1 l1) else None)"
  2404     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1" by blast
  2405   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
  2406   moreover have l1_le_u1: "real l1 \<le> real u1" using l1 u1 by auto
  2407   ultimately have "real l1 \<noteq> 0" and "real u1 \<noteq> 0" unfolding less_float_def by auto
  2408 
  2409   have inv: "inverse (real u1) \<le> inverse (interpret_floatarith a xs)
  2410            \<and> inverse (interpret_floatarith a xs) \<le> inverse (real l1)"
  2411   proof (cases "0 < l1")
  2412     case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs"
  2413       unfolding less_float_def using l1_le_u1 l1 by auto
  2414     show ?thesis
  2415       unfolding inverse_le_iff_le[OF `0 < real u1` `0 < interpret_floatarith a xs`]
  2416         inverse_le_iff_le[OF `0 < interpret_floatarith a xs` `0 < real l1`]
  2417       using l1 u1 by auto
  2418   next
  2419     case False hence "u1 < 0" using either by blast
  2420     hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0"
  2421       unfolding less_float_def using l1_le_u1 u1 by auto
  2422     show ?thesis
  2423       unfolding inverse_le_iff_le_neg[OF `real u1 < 0` `interpret_floatarith a xs < 0`]
  2424         inverse_le_iff_le_neg[OF `interpret_floatarith a xs < 0` `real l1 < 0`]
  2425       using l1 u1 by auto
  2426   qed
  2427 
  2428   from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \<or> u1 < 0", auto)
  2429   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
  2430   also have "\<dots> \<le> inverse (interpret_floatarith a xs)" using inv by auto
  2431   finally have "real l \<le> inverse (interpret_floatarith a xs)" .
  2432   moreover
  2433   from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \<or> u1 < 0", auto)
  2434   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
  2435   hence "inverse (interpret_floatarith a xs) \<le> real u" by (rule order_trans[OF inv[THEN conjunct2]])
  2436   ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto
  2437 next
  2438   case (Abs x)
  2439   from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
  2440   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>"
  2441     and l1: "real l1 \<le> interpret_floatarith x xs" and u1: "interpret_floatarith x xs \<le> real u1" by blast
  2442   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)
  2443 next
  2444   case (Min a b)
  2445   from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
  2446   obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
  2447     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
  2448     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
  2449   thus ?case unfolding l' u' by (auto simp add: real_of_float_min)
  2450 next
  2451   case (Max a b)
  2452   from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
  2453   obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
  2454     and l1: "real l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> real u1"
  2455     and l1: "real l2 \<le> interpret_floatarith b xs" and u1: "interpret_floatarith b xs \<le> real u2" by blast
  2456   thus ?case unfolding l' u' by (auto simp add: real_of_float_max)
  2457 next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto
  2458 next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto
  2459 next case Pi with pi_boundaries show ?case by auto
  2460 next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
  2461 next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto
  2462 next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto
  2463 next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto
  2464 next case (Num f) thus ?case by auto
  2465 next
  2466   case (Var n)
  2467   from this[symmetric] `bounded_by xs vs`[THEN bounded_byE, of n]
  2468   show ?case by (cases "n < length vs", auto)
  2469 qed
  2470 
  2471 datatype form = Bound floatarith floatarith floatarith form
  2472               | Assign floatarith floatarith form
  2473               | Less floatarith floatarith
  2474               | LessEqual floatarith floatarith
  2475               | AtLeastAtMost floatarith floatarith floatarith
  2476 
  2477 fun interpret_form :: "form \<Rightarrow> real list \<Rightarrow> bool" where
  2478 "interpret_form (Bound x a b f) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs } \<longrightarrow> interpret_form f vs)" |
  2479 "interpret_form (Assign x a f) vs  = (interpret_floatarith x vs = interpret_floatarith a vs \<longrightarrow> interpret_form f vs)" |
  2480 "interpret_form (Less a b) vs      = (interpret_floatarith a vs < interpret_floatarith b vs)" |
  2481 "interpret_form (LessEqual a b) vs = (interpret_floatarith a vs \<le> interpret_floatarith b vs)" |
  2482 "interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs })"
  2483 
  2484 fun approx_form' and approx_form :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> nat list \<Rightarrow> bool" where
  2485 "approx_form' prec f 0 n l u bs ss = approx_form prec f (bs[n := Some (l, u)]) ss" |
  2486 "approx_form' prec f (Suc s) n l u bs ss =
  2487   (let m = (l + u) * Float 1 -1
  2488    in (if approx_form' prec f s n l m bs ss then approx_form' prec f s n m u bs ss else False))" |
  2489 "approx_form prec (Bound (Var n) a b f) bs ss =
  2490    (case (approx prec a bs, approx prec b bs)
  2491    of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
  2492     | _ \<Rightarrow> False)" |
  2493 "approx_form prec (Assign (Var n) a f) bs ss =
  2494    (case (approx prec a bs)
  2495    of (Some (l, u)) \<Rightarrow> approx_form' prec f (ss ! n) n l u bs ss
  2496     | _ \<Rightarrow> False)" |
  2497 "approx_form prec (Less a b) bs ss =
  2498    (case (approx prec a bs, approx prec b bs)
  2499    of (Some (l, u), Some (l', u')) \<Rightarrow> u < l'
  2500     | _ \<Rightarrow> False)" |
  2501 "approx_form prec (LessEqual a b) bs ss =
  2502    (case (approx prec a bs, approx prec b bs)
  2503    of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l'
  2504     | _ \<Rightarrow> False)" |
  2505 "approx_form prec (AtLeastAtMost x a b) bs ss =
  2506    (case (approx prec x bs, approx prec a bs, approx prec b bs)
  2507    of (Some (lx, ux), Some (l, u), Some (l', u')) \<Rightarrow> u \<le> lx \<and> ux \<le> l'
  2508     | _ \<Rightarrow> False)" |
  2509 "approx_form _ _ _ _ = False"
  2510 
  2511 lemma lazy_conj: "(if A then B else False) = (A \<and> B)" by simp
  2512 
  2513 lemma approx_form_approx_form':
  2514   assumes "approx_form' prec f s n l u bs ss" and "x \<in> { real l .. real u }"
  2515   obtains l' u' where "x \<in> { real l' .. real u' }"
  2516   and "approx_form prec f (bs[n := Some (l', u')]) ss"
  2517 using assms proof (induct s arbitrary: l u)
  2518   case 0
  2519   from this(1)[of l u] this(2,3)
  2520   show thesis by auto
  2521 next
  2522   case (Suc s)
  2523 
  2524   let ?m = "(l + u) * Float 1 -1"
  2525   have "real l \<le> real ?m" and "real ?m \<le> real u"
  2526     unfolding le_float_def using Suc.prems by auto
  2527 
  2528   with `x \<in> { real l .. real u }`
  2529   have "x \<in> { real l .. real ?m} \<or> x \<in> { real ?m .. real u }" by auto
  2530   thus thesis
  2531   proof (rule disjE)
  2532     assume *: "x \<in> { real l .. real ?m }"
  2533     with Suc.hyps[OF _ _ *] Suc.prems
  2534     show thesis by (simp add: Let_def lazy_conj)
  2535   next
  2536     assume *: "x \<in> { real ?m .. real u }"
  2537     with Suc.hyps[OF _ _ *] Suc.prems
  2538     show thesis by (simp add: Let_def lazy_conj)
  2539   qed
  2540 qed
  2541 
  2542 lemma approx_form_aux:
  2543   assumes "approx_form prec f vs ss"
  2544   and "bounded_by xs vs"
  2545   shows "interpret_form f xs"
  2546 using assms proof (induct f arbitrary: vs)
  2547   case (Bound x a b f)
  2548   then obtain n
  2549     where x_eq: "x = Var n" by (cases x) auto
  2550 
  2551   with Bound.prems obtain l u' l' u
  2552     where l_eq: "Some (l, u') = approx prec a vs"
  2553     and u_eq: "Some (l', u) = approx prec b vs"
  2554     and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
  2555     by (cases "approx prec a vs", simp,
  2556         cases "approx prec b vs", auto) blast
  2557 
  2558   { assume "xs ! n \<in> { interpret_floatarith a xs .. interpret_floatarith b xs }"
  2559     with approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq]
  2560     have "xs ! n \<in> { real l .. real u}" by auto
  2561 
  2562     from approx_form_approx_form'[OF approx_form' this]
  2563     obtain lx ux where bnds: "xs ! n \<in> { real lx .. real ux }"
  2564       and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
  2565 
  2566     from `bounded_by xs vs` bnds
  2567     have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
  2568     with Bound.hyps[OF approx_form]
  2569     have "interpret_form f xs" by blast }
  2570   thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
  2571 next
  2572   case (Assign x a f)
  2573   then obtain n
  2574     where x_eq: "x = Var n" by (cases x) auto
  2575 
  2576   with Assign.prems obtain l u' l' u
  2577     where bnd_eq: "Some (l, u) = approx prec a vs"
  2578     and x_eq: "x = Var n"
  2579     and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
  2580     by (cases "approx prec a vs") auto
  2581 
  2582   { assume bnds: "xs ! n = interpret_floatarith a xs"
  2583     with approx[OF Assign.prems(2) bnd_eq]
  2584     have "xs ! n \<in> { real l .. real u}" by auto
  2585     from approx_form_approx_form'[OF approx_form' this]
  2586     obtain lx ux where bnds: "xs ! n \<in> { real lx .. real ux }"
  2587       and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .
  2588 
  2589     from `bounded_by xs vs` bnds
  2590     have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update)
  2591     with Assign.hyps[OF approx_form]
  2592     have "interpret_form f xs" by blast }
  2593   thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp
  2594 next
  2595   case (Less a b)
  2596   then obtain l u l' u'
  2597     where l_eq: "Some (l, u) = approx prec a vs"
  2598     and u_eq: "Some (l', u') = approx prec b vs"
  2599     and inequality: "u < l'"
  2600     by (cases "approx prec a vs", auto,
  2601       cases "approx prec b vs", auto)
  2602   from inequality[unfolded less_float_def] approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
  2603   show ?case by auto
  2604 next
  2605   case (LessEqual a b)
  2606   then obtain l u l' u'
  2607     where l_eq: "Some (l, u) = approx prec a vs"
  2608     and u_eq: "Some (l', u') = approx prec b vs"
  2609     and inequality: "u \<le> l'"
  2610     by (cases "approx prec a vs", auto,
  2611       cases "approx prec b vs", auto)
  2612   from inequality[unfolded le_float_def] approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
  2613   show ?case by auto
  2614 next
  2615   case (AtLeastAtMost x a b)
  2616   then obtain lx ux l u l' u'
  2617     where x_eq: "Some (lx, ux) = approx prec x vs"
  2618     and l_eq: "Some (l, u) = approx prec a vs"
  2619     and u_eq: "Some (l', u') = approx prec b vs"
  2620     and inequality: "u \<le> lx \<and> ux \<le> l'"
  2621     by (cases "approx prec x vs", auto,
  2622       cases "approx prec a vs", auto,
  2623       cases "approx prec b vs", auto, blast)
  2624   from inequality[unfolded le_float_def] approx[OF AtLeastAtMost.prems(2) l_eq] approx[OF AtLeastAtMost.prems(2) u_eq] approx[OF AtLeastAtMost.prems(2) x_eq]
  2625   show ?case by auto
  2626 qed
  2627 
  2628 lemma approx_form:
  2629   assumes "n = length xs"
  2630   assumes "approx_form prec f (replicate n None) ss"
  2631   shows "interpret_form f xs"
  2632   using approx_form_aux[OF _ bounded_by_None] assms by auto
  2633 
  2634 subsection {* Implementing Taylor series expansion *}
  2635 
  2636 fun isDERIV :: "nat \<Rightarrow> floatarith \<Rightarrow> real list \<Rightarrow> bool" where
  2637 "isDERIV x (Add a b) vs         = (isDERIV x a vs \<and> isDERIV x b vs)" |
  2638 "isDERIV x (Mult a b) vs        = (isDERIV x a vs \<and> isDERIV x b vs)" |
  2639 "isDERIV x (Minus a) vs         = isDERIV x a vs" |
  2640 "isDERIV x (Inverse a) vs       = (isDERIV x a vs \<and> interpret_floatarith a vs \<noteq> 0)" |
  2641 "isDERIV x (Cos a) vs           = isDERIV x a vs" |
  2642 "isDERIV x (Arctan a) vs        = isDERIV x a vs" |
  2643 "isDERIV x (Min a b) vs         = False" |
  2644 "isDERIV x (Max a b) vs         = False" |
  2645 "isDERIV x (Abs a) vs           = False" |
  2646 "isDERIV x Pi vs                = True" |
  2647 "isDERIV x (Sqrt a) vs          = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
  2648 "isDERIV x (Exp a) vs           = isDERIV x a vs" |
  2649 "isDERIV x (Ln a) vs            = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
  2650 "isDERIV x (Power a 0) vs       = True" |
  2651 "isDERIV x (Power a (Suc n)) vs = isDERIV x a vs" |
  2652 "isDERIV x (Num f) vs           = True" |
  2653 "isDERIV x (Var n) vs          = True"
  2654 
  2655 fun DERIV_floatarith :: "nat \<Rightarrow> floatarith \<Rightarrow> floatarith" where
  2656 "DERIV_floatarith x (Add a b)         = Add (DERIV_floatarith x a) (DERIV_floatarith x b)" |
  2657 "DERIV_floatarith x (Mult a b)        = Add (Mult a (DERIV_floatarith x b)) (Mult (DERIV_floatarith x a) b)" |
  2658 "DERIV_floatarith x (Minus a)         = Minus (DERIV_floatarith x a)" |
  2659 "DERIV_floatarith x (Inverse a)       = Minus (Mult (DERIV_floatarith x a) (Inverse (Power a 2)))" |
  2660 "DERIV_floatarith x (Cos a)           = Minus (Mult (Cos (Add (Mult Pi (Num (Float 1 -1))) (Minus a))) (DERIV_floatarith x a))" |
  2661 "DERIV_floatarith x (Arctan a)        = Mult (Inverse (Add (Num 1) (Power a 2))) (DERIV_floatarith x a)" |
  2662 "DERIV_floatarith x (Min a b)         = Num 0" |
  2663 "DERIV_floatarith x (Max a b)         = Num 0" |
  2664 "DERIV_floatarith x (Abs a)           = Num 0" |
  2665 "DERIV_floatarith x Pi                = Num 0" |
  2666 "DERIV_floatarith x (Sqrt a)          = (Mult (Inverse (Mult (Sqrt a) (Num 2))) (DERIV_floatarith x a))" |
  2667 "DERIV_floatarith x (Exp a)           = Mult (Exp a) (DERIV_floatarith x a)" |
  2668 "DERIV_floatarith x (Ln a)            = Mult (Inverse a) (DERIV_floatarith x a)" |
  2669 "DERIV_floatarith x (Power a 0)       = Num 0" |
  2670 "DERIV_floatarith x (Power a (Suc n)) = Mult (Num (Float (int (Suc n)) 0)) (Mult (Power a n) (DERIV_floatarith x a))" |
  2671 "DERIV_floatarith x (Num f)           = Num 0" |
  2672 "DERIV_floatarith x (Var n)          = (if x = n then Num 1 else Num 0)"
  2673 
  2674 lemma DERIV_floatarith:
  2675   assumes "n < length vs"
  2676   assumes isDERIV: "isDERIV n f (vs[n := x])"
  2677   shows "DERIV (\<lambda> x'. interpret_floatarith f (vs[n := x'])) x :>
  2678                interpret_floatarith (DERIV_floatarith n f) (vs[n := x])"
  2679    (is "DERIV (?i f) x :> _")
  2680 using isDERIV proof (induct f arbitrary: x)
  2681      case (Inverse a) thus ?case
  2682     by (auto intro!: DERIV_intros
  2683              simp add: algebra_simps power2_eq_square)
  2684 next case (Cos a) thus ?case
  2685   by (auto intro!: DERIV_intros
  2686            simp del: interpret_floatarith.simps(5)
  2687            simp add: interpret_floatarith_sin interpret_floatarith.simps(5)[of a])
  2688 next case (Power a n) thus ?case
  2689   by (cases n, auto intro!: DERIV_intros
  2690                     simp del: power_Suc simp add: real_eq_of_nat)
  2691 next case (Ln a) thus ?case
  2692     by (auto intro!: DERIV_intros simp add: divide_inverse)
  2693 next case (Var i) thus ?case using `n < length vs` by auto
  2694 qed (auto intro!: DERIV_intros)
  2695 
  2696 declare approx.simps[simp del]
  2697 
  2698 fun isDERIV_approx :: "nat \<Rightarrow> nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> bool" where
  2699 "isDERIV_approx prec x (Add a b) vs         = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
  2700 "isDERIV_approx prec x (Mult a b) vs        = (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs)" |
  2701 "isDERIV_approx prec x (Minus a) vs         = isDERIV_approx prec x a vs" |
  2702 "isDERIV_approx prec x (Inverse a) vs       =
  2703   (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l \<or> u < 0 | None \<Rightarrow> False))" |
  2704 "isDERIV_approx prec x (Cos a) vs           = isDERIV_approx prec x a vs" |
  2705 "isDERIV_approx prec x (Arctan a) vs        = isDERIV_approx prec x a vs" |
  2706 "isDERIV_approx prec x (Min a b) vs         = False" |
  2707 "isDERIV_approx prec x (Max a b) vs         = False" |
  2708 "isDERIV_approx prec x (Abs a) vs           = False" |
  2709 "isDERIV_approx prec x Pi vs                = True" |
  2710 "isDERIV_approx prec x (Sqrt a) vs          =
  2711   (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
  2712 "isDERIV_approx prec x (Exp a) vs           = isDERIV_approx prec x a vs" |
  2713 "isDERIV_approx prec x (Ln a) vs            =
  2714   (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
  2715 "isDERIV_approx prec x (Power a 0) vs       = True" |
  2716 "isDERIV_approx prec x (Power a (Suc n)) vs = isDERIV_approx prec x a vs" |
  2717 "isDERIV_approx prec x (Num f) vs           = True" |
  2718 "isDERIV_approx prec x (Var n) vs          = True"
  2719 
  2720 lemma isDERIV_approx:
  2721   assumes "bounded_by xs vs"
  2722   and isDERIV_approx: "isDERIV_approx prec x f vs"
  2723   shows "isDERIV x f xs"
  2724 using isDERIV_approx proof (induct f)
  2725   case (Inverse a)
  2726   then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
  2727     and *: "0 < l \<or> u < 0"
  2728     by (cases "approx prec a vs", auto)
  2729   with approx[OF `bounded_by xs vs` approx_Some]
  2730   have "interpret_floatarith a xs \<noteq> 0" unfolding less_float_def by auto
  2731   thus ?case using Inverse by auto
  2732 next
  2733   case (Ln a)
  2734   then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
  2735     and *: "0 < l"
  2736     by (cases "approx prec a vs", auto)
  2737   with approx[OF `bounded_by xs vs` approx_Some]
  2738   have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
  2739   thus ?case using Ln by auto
  2740 next
  2741   case (Sqrt a)
  2742   then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
  2743     and *: "0 < l"
  2744     by (cases "approx prec a vs", auto)
  2745   with approx[OF `bounded_by xs vs` approx_Some]
  2746   have "0 < interpret_floatarith a xs" unfolding less_float_def by auto
  2747   thus ?case using Sqrt by auto
  2748 next
  2749   case (Power a n) thus ?case by (cases n, auto)
  2750 qed auto
  2751 
  2752 lemma bounded_by_update_var:
  2753   assumes "bounded_by xs vs" and "vs ! i = Some (l, u)"
  2754   and bnd: "x \<in> { real l .. real u }"
  2755   shows "bounded_by (xs[i := x]) vs"
  2756 proof (cases "i < length xs")
  2757   case False thus ?thesis using `bounded_by xs vs` by auto
  2758 next
  2759   let ?xs = "xs[i := x]"
  2760   case True hence "i < length ?xs" by auto
  2761 { fix j
  2762   assume "j < length vs"
  2763   have "case vs ! j of None \<Rightarrow> True | Some (l, u) \<Rightarrow> ?xs ! j \<in> { real l .. real u }"
  2764   proof (cases "vs ! j")
  2765     case (Some b)
  2766     thus ?thesis
  2767     proof (cases "i = j")
  2768       case True
  2769       thus ?thesis using `vs ! i = Some (l, u)` Some and bnd `i < length ?xs`
  2770         by auto
  2771     next
  2772       case False
  2773       thus ?thesis using `bounded_by xs vs`[THEN bounded_byE, OF `j < length vs`] Some
  2774         by auto
  2775     qed
  2776   qed auto }
  2777   thus ?thesis unfolding bounded_by_def by auto
  2778 qed
  2779 
  2780 lemma isDERIV_approx':
  2781   assumes "bounded_by xs vs"
  2782   and vs_x: "vs ! x = Some (l, u)" and X_in: "X \<in> { real l .. real u }"
  2783   and approx: "isDERIV_approx prec x f vs"
  2784   shows "isDERIV x f (xs[x := X])"
  2785 proof -
  2786   note bounded_by_update_var[OF `bounded_by xs vs` vs_x X_in] approx
  2787   thus ?thesis by (rule isDERIV_approx)
  2788 qed
  2789 
  2790 lemma DERIV_approx:
  2791   assumes "n < length xs" and bnd: "bounded_by xs vs"
  2792   and isD: "isDERIV_approx prec n f vs"
  2793   and app: "Some (l, u) = approx prec (DERIV_floatarith n f) vs" (is "_ = approx _ ?D _")
  2794   shows "\<exists>x. real l \<le> x \<and> x \<le> real u \<and>
  2795              DERIV (\<lambda> x. interpret_floatarith f (xs[n := x])) (xs!n) :> x"
  2796          (is "\<exists> x. _ \<and> _ \<and> DERIV (?i f) _ :> _")
  2797 proof (rule exI[of _ "?i ?D (xs!n)"], rule conjI[OF _ conjI])
  2798   let "?i f x" = "interpret_floatarith f (xs[n := x])"
  2799   from approx[OF bnd app]
  2800   show "real l \<le> ?i ?D (xs!n)" and "?i ?D (xs!n) \<le> real u"
  2801     using `n < length xs` by auto
  2802   from DERIV_floatarith[OF `n < length xs`, of f "xs!n"] isDERIV_approx[OF bnd isD]
  2803   show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))" by simp
  2804 qed
  2805 
  2806 fun lift_bin :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float) option) \<Rightarrow> (float * float) option" where
  2807 "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2" |
  2808 "lift_bin a b f = None"
  2809 
  2810 lemma lift_bin:
  2811   assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
  2812   obtains l1 u1 l2 u2
  2813   where "a = Some (l1, u1)"
  2814   and "b = Some (l2, u2)"
  2815   and "f l1 u1 l2 u2 = Some (l, u)"
  2816 using assms by (cases a, simp, cases b, simp, auto)
  2817 
  2818 fun approx_tse where
  2819 "approx_tse prec n 0 c k f bs = approx prec f bs" |
  2820 "approx_tse prec n (Suc s) c k f bs =
  2821   (if isDERIV_approx prec n f bs then
  2822     lift_bin (approx prec f (bs[n := Some (c,c)]))
  2823              (approx_tse prec n s c (Suc k) (DERIV_floatarith n f) bs)
  2824              (\<lambda> l1 u1 l2 u2. approx prec
  2825                  (Add (Var 0)
  2826                       (Mult (Inverse (Num (Float (int k) 0)))
  2827                                  (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
  2828                                        (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), bs!n])
  2829   else approx prec f bs)"
  2830 
  2831 lemma bounded_by_Cons:
  2832   assumes bnd: "bounded_by xs vs"
  2833   and x: "x \<in> { real l .. real u }"
  2834   shows "bounded_by (x#xs) ((Some (l, u))#vs)"
  2835 proof -
  2836   { fix i assume *: "i < length ((Some (l, u))#vs)"
  2837     have "case ((Some (l,u))#vs) ! i of Some (l, u) \<Rightarrow> (x#xs)!i \<in> { real l .. real u } | None \<Rightarrow> True"
  2838     proof (cases i)
  2839       case 0 with x show ?thesis by auto
  2840     next
  2841       case (Suc i) with * have "i < length vs" by auto
  2842       from bnd[THEN bounded_byE, OF this]
  2843       show ?thesis unfolding Suc nth_Cons_Suc .
  2844     qed }
  2845   thus ?thesis by (auto simp add: bounded_by_def)
  2846 qed
  2847 
  2848 lemma approx_tse_generic:
  2849   assumes "bounded_by xs vs"
  2850   and bnd_c: "bounded_by (xs[x := real c]) vs" and "x < length vs" and "x < length xs"
  2851   and bnd_x: "vs ! x = Some (lx, ux)"
  2852   and ate: "Some (l, u) = approx_tse prec x s c k f vs"
  2853   shows "\<exists> n. (\<forall> m < n. \<forall> z \<in> {real lx .. real ux}.
  2854       DERIV (\<lambda> y. interpret_floatarith ((DERIV_floatarith x ^^ m) f) (xs[x := y])) z :>
  2855             (interpret_floatarith ((DERIV_floatarith x ^^ (Suc m)) f) (xs[x := z])))
  2856    \<and> (\<forall> t \<in> {real lx .. real ux}.  (\<Sum> i = 0..<n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) *
  2857                   interpret_floatarith ((DERIV_floatarith x ^^ i) f) (xs[x := real c]) *
  2858                   (xs!x - real c)^i) +
  2859       inverse (real (\<Prod> j \<in> {k..<k+n}. j)) *
  2860       interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := t]) *
  2861       (xs!x - real c)^n \<in> {real l .. real u})" (is "\<exists> n. ?taylor f k l u n")
  2862 using ate proof (induct s arbitrary: k f l u)
  2863   case 0
  2864   { fix t assume "t \<in> {real lx .. real ux}"
  2865     note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
  2866     from approx[OF this 0[unfolded approx_tse.simps]]
  2867     have "(interpret_floatarith f (xs[x := t])) \<in> {real l .. real u}"
  2868       by (auto simp add: algebra_simps)
  2869   } thus ?case by (auto intro!: exI[of _ 0])
  2870 next
  2871   case (Suc s)
  2872   show ?case
  2873   proof (cases "isDERIV_approx prec x f vs")
  2874     case False
  2875     note ap = Suc.prems[unfolded approx_tse.simps if_not_P[OF False]]
  2876 
  2877     { fix t assume "t \<in> {real lx .. real ux}"
  2878       note bounded_by_update_var[OF `bounded_by xs vs` bnd_x this]
  2879       from approx[OF this ap]
  2880       have "(interpret_floatarith f (xs[x := t])) \<in> {real l .. real u}"
  2881         by (auto simp add: algebra_simps)
  2882     } thus ?thesis by (auto intro!: exI[of _ 0])
  2883   next
  2884     case True
  2885     with Suc.prems
  2886     obtain l1 u1 l2 u2
  2887       where a: "Some (l1, u1) = approx prec f (vs[x := Some (c,c)])"
  2888       and ate: "Some (l2, u2) = approx_tse prec x s c (Suc k) (DERIV_floatarith x f) vs"
  2889       and final: "Some (l, u) = approx prec
  2890         (Add (Var 0)
  2891              (Mult (Inverse (Num (Float (int k) 0)))
  2892                    (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
  2893                          (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), vs!x]"
  2894       by (auto elim!: lift_bin) blast
  2895 
  2896     from bnd_c `x < length xs`
  2897     have bnd: "bounded_by (xs[x:=real c]) (vs[x:= Some (c,c)])"
  2898       by (auto intro!: bounded_by_update)
  2899 
  2900     from approx[OF this a]
  2901     have f_c: "interpret_floatarith ((DERIV_floatarith x ^^ 0) f) (xs[x := real c]) \<in> { real l1 .. real u1 }"
  2902               (is "?f 0 (real c) \<in> _")
  2903       by auto
  2904 
  2905     { fix f :: "'a \<Rightarrow> 'a" fix n :: nat fix x :: 'a
  2906       have "(f ^^ Suc n) x = (f ^^ n) (f x)"
  2907         by (induct n, auto) }
  2908     note funpow_Suc = this[symmetric]
  2909     from Suc.hyps[OF ate, unfolded this]
  2910     obtain n
  2911       where DERIV_hyp: "\<And> m z. \<lbrakk> m < n ; z \<in> { real lx .. real ux } \<rbrakk> \<Longrightarrow> DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z"
  2912       and hyp: "\<forall> t \<in> {real lx .. real ux}. (\<Sum> i = 0..<n. inverse (real (\<Prod> j \<in> {Suc k..<Suc k + i}. j)) * ?f (Suc i) (real c) * (xs!x - real c)^i) +
  2913            inverse (real (\<Prod> j \<in> {Suc k..<Suc k + n}. j)) * ?f (Suc n) t * (xs!x - real c)^n \<in> {real l2 .. real u2}"
  2914           (is "\<forall> t \<in> _. ?X (Suc k) f n t \<in> _")
  2915       by blast
  2916 
  2917     { fix m z
  2918       assume "m < Suc n" and bnd_z: "z \<in> { real lx .. real ux }"
  2919       have "DERIV (?f m) z :> ?f (Suc m) z"
  2920       proof (cases m)
  2921         case 0
  2922         with DERIV_floatarith[OF `x < length xs` isDERIV_approx'[OF `bounded_by xs vs` bnd_x bnd_z True]]
  2923         show ?thesis by simp
  2924       next
  2925         case (Suc m')
  2926         hence "m' < n" using `m < Suc n` by auto
  2927         from DERIV_hyp[OF this bnd_z]
  2928         show ?thesis using Suc by simp
  2929       qed } note DERIV = this
  2930 
  2931     have "\<And> k i. k < i \<Longrightarrow> {k ..< i} = insert k {Suc k ..< i}" by auto
  2932     hence setprod_head_Suc: "\<And> k i. \<Prod> {k ..< k + Suc i} = k * \<Prod> {Suc k ..< Suc k + i}" by auto
  2933     have setsum_move0: "\<And> k F. setsum F {0..<Suc k} = F 0 + setsum (\<lambda> k. F (Suc k)) {0..<k}"
  2934       unfolding setsum_shift_bounds_Suc_ivl[symmetric]
  2935       unfolding setsum_head_upt_Suc[OF zero_less_Suc] ..
  2936     def C \<equiv> "xs!x - real c"
  2937 
  2938     { fix t assume t: "t \<in> {real lx .. real ux}"
  2939       hence "bounded_by [xs!x] [vs!x]"
  2940         using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`]
  2941         by (cases "vs!x", auto simp add: bounded_by_def)
  2942 
  2943       with hyp[THEN bspec, OF t] f_c
  2944       have "bounded_by [?f 0 (real c), ?X (Suc k) f n t, xs!x] [Some (l1, u1), Some (l2, u2), vs!x]"
  2945         by (auto intro!: bounded_by_Cons)
  2946       from approx[OF this final, unfolded atLeastAtMost_iff[symmetric]]
  2947       have "?X (Suc k) f n t * (xs!x - real c) * inverse (real k) + ?f 0 (real c) \<in> {real l .. real u}"
  2948         by (auto simp add: algebra_simps)
  2949       also have "?X (Suc k) f n t * (xs!x - real c) * inverse (real k) + ?f 0 (real c) =
  2950                (\<Sum> i = 0..<Suc n. inverse (real (\<Prod> j \<in> {k..<k+i}. j)) * ?f i (real c) * (xs!x - real c)^i) +
  2951                inverse (real (\<Prod> j \<in> {k..<k+Suc n}. j)) * ?f (Suc n) t * (xs!x - real c)^Suc n" (is "_ = ?T")
  2952         unfolding funpow_Suc C_def[symmetric] setsum_move0 setprod_head_Suc
  2953         by (auto simp add: algebra_simps)
  2954           (simp only: mult_left_commute [of _ "inverse (real k)"] setsum_right_distrib [symmetric])
  2955       finally have "?T \<in> {real l .. real u}" . }
  2956     thus ?thesis using DERIV by blast
  2957   qed
  2958 qed
  2959 
  2960 lemma setprod_fact: "\<Prod> {1..<1 + k} = fact (k :: nat)"
  2961 proof (induct k)
  2962   case (Suc k)
  2963   have "{ 1 ..< Suc (Suc k) } = insert (Suc k) { 1 ..< Suc k }" by auto
  2964   hence "\<Prod> { 1 ..< Suc (Suc k) } = (Suc k) * \<Prod> { 1 ..< Suc k }" by auto
  2965   thus ?case using Suc by auto
  2966 qed simp
  2967 
  2968 lemma approx_tse:
  2969   assumes "bounded_by xs vs"
  2970   and bnd_x: "vs ! x = Some (lx, ux)" and bnd_c: "real c \<in> {real lx .. real ux}"
  2971   and "x < length vs" and "x < length xs"
  2972   and ate: "Some (l, u) = approx_tse prec x s c 1 f vs"
  2973   shows "interpret_floatarith f xs \<in> { real l .. real u }"
  2974 proof -
  2975   def F \<equiv> "\<lambda> n z. interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])"
  2976   hence F0: "F 0 = (\<lambda> z. interpret_floatarith f (xs[x := z]))" by auto
  2977 
  2978   hence "bounded_by (xs[x := real c]) vs" and "x < length vs" "x < length xs"
  2979     using `bounded_by xs vs` bnd_x bnd_c `x < length vs` `x < length xs`
  2980     by (auto intro!: bounded_by_update_var)
  2981 
  2982   from approx_tse_generic[OF `bounded_by xs vs` this bnd_x ate]
  2983   obtain n
  2984     where DERIV: "\<forall> m z. m < n \<and> real lx \<le> z \<and> z \<le> real ux \<longrightarrow> DERIV (F m) z :> F (Suc m) z"
  2985     and hyp: "\<And> t. t \<in> {real lx .. real ux} \<Longrightarrow>
  2986            (\<Sum> j = 0..<n. inverse (real (fact j)) * F j (real c) * (xs!x - real c)^j) +
  2987              inverse (real (fact n)) * F n t * (xs!x - real c)^n
  2988              \<in> {real l .. real u}" (is "\<And> t. _ \<Longrightarrow> ?taylor t \<in> _")
  2989     unfolding F_def atLeastAtMost_iff[symmetric] setprod_fact by blast
  2990 
  2991   have bnd_xs: "xs ! x \<in> { real lx .. real ux }"
  2992     using `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
  2993 
  2994   show ?thesis
  2995   proof (cases n)
  2996     case 0 thus ?thesis using hyp[OF bnd_xs] unfolding F_def by auto
  2997   next
  2998     case (Suc n')
  2999     show ?thesis
  3000     proof (cases "xs ! x = real c")
  3001       case True
  3002       from True[symmetric] hyp[OF bnd_xs] Suc show ?thesis
  3003         unfolding F_def Suc setsum_head_upt_Suc[OF zero_less_Suc] setsum_shift_bounds_Suc_ivl by auto
  3004     next
  3005       case False
  3006 
  3007       have "real lx \<le> real c" "real c \<le> real ux" "real lx \<le> xs!x" "xs!x \<le> real ux"
  3008         using Suc bnd_c `bounded_by xs vs`[THEN bounded_byE, OF `x < length vs`] bnd_x by auto
  3009       from Taylor.taylor[OF zero_less_Suc, of F, OF F0 DERIV[unfolded Suc] this False]
  3010       obtain t where t_bnd: "if xs ! x < real c then xs ! x < t \<and> t < real c else real c < t \<and> t < xs ! x"
  3011         and fl_eq: "interpret_floatarith f (xs[x := xs ! x]) =
  3012            (\<Sum>m = 0..<Suc n'. F m (real c) / real (fact m) * (xs ! x - real c) ^ m) +
  3013            F (Suc n') t / real (fact (Suc n')) * (xs ! x - real c) ^ Suc n'"
  3014         by blast
  3015 
  3016       from t_bnd bnd_xs bnd_c have *: "t \<in> {real lx .. real ux}"
  3017         by (cases "xs ! x < real c", auto)
  3018 
  3019       have "interpret_floatarith f (xs[x := xs ! x]) = ?taylor t"
  3020         unfolding fl_eq Suc by (auto simp add: algebra_simps divide_inverse)
  3021       also have "\<dots> \<in> {real l .. real u}" using * by (rule hyp)
  3022       finally show ?thesis by simp
  3023     qed
  3024   qed
  3025 qed
  3026 
  3027 fun approx_tse_form' where
  3028 "approx_tse_form' prec t f 0 l u cmp =
  3029   (case approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)]
  3030      of Some (l, u) \<Rightarrow> cmp l u | None \<Rightarrow> False)" |
  3031 "approx_tse_form' prec t f (Suc s) l u cmp =
  3032   (let m = (l + u) * Float 1 -1
  3033    in (if approx_tse_form' prec t f s l m cmp then
  3034       approx_tse_form' prec t f s m u cmp else False))"
  3035 
  3036 lemma approx_tse_form':
  3037   assumes "approx_tse_form' prec t f s l u cmp" and "x \<in> {real l .. real u}"
  3038   shows "\<exists> l' u' ly uy. x \<in> { real l' .. real u' } \<and> real l \<le> real l' \<and> real u' \<le> real u \<and> cmp ly uy \<and>
  3039                   approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)"
  3040 using assms proof (induct s arbitrary: l u)
  3041   case 0
  3042   then obtain ly uy
  3043     where *: "approx_tse prec 0 t ((l + u) * Float 1 -1) 1 f [Some (l, u)] = Some (ly, uy)"
  3044     and **: "cmp ly uy" by (auto elim!: option_caseE)
  3045   with 0 show ?case by (auto intro!: exI)
  3046 next
  3047   case (Suc s)
  3048   let ?m = "(l + u) * Float 1 -1"
  3049   from Suc.prems
  3050   have l: "approx_tse_form' prec t f s l ?m cmp"
  3051     and u: "approx_tse_form' prec t f s ?m u cmp"
  3052     by (auto simp add: Let_def lazy_conj)
  3053 
  3054   have m_l: "real l \<le> real ?m" and m_u: "real ?m \<le> real u"
  3055     unfolding le_float_def using Suc.prems by auto
  3056 
  3057   with `x \<in> { real l .. real u }`
  3058   have "x \<in> { real l .. real ?m} \<or> x \<in> { real ?m .. real u }" by auto
  3059   thus ?case
  3060   proof (rule disjE)
  3061     assume "x \<in> { real l .. real ?m}"
  3062     from Suc.hyps[OF l this]
  3063     obtain l' u' ly uy
  3064       where "x \<in> { real l' .. real u' } \<and> real l \<le> real l' \<and> real u' \<le> real ?m \<and> cmp ly uy \<and>
  3065                   approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
  3066     with m_u show ?thesis by (auto intro!: exI)
  3067   next
  3068     assume "x \<in> { real ?m .. real u }"
  3069     from Suc.hyps[OF u this]
  3070     obtain l' u' ly uy
  3071       where "x \<in> { real l' .. real u' } \<and> real ?m \<le> real l' \<and> real u' \<le> real u \<and> cmp ly uy \<and>
  3072                   approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 f [Some (l', u')] = Some (ly, uy)" by blast
  3073     with m_u show ?thesis by (auto intro!: exI)
  3074   qed
  3075 qed
  3076 
  3077 lemma approx_tse_form'_less:
  3078   assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 < l)"
  3079   and x: "x \<in> {real l .. real u}"
  3080   shows "interpret_floatarith b [x] < interpret_floatarith a [x]"
  3081 proof -
  3082   from approx_tse_form'[OF tse x]
  3083   obtain l' u' ly uy
  3084     where x': "x \<in> { real l' .. real u' }" and "real l \<le> real l'"
  3085     and "real u' \<le> real u" and "0 < ly"
  3086     and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
  3087     by blast
  3088 
  3089   hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
  3090 
  3091   from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
  3092   have "real ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
  3093     by (auto simp add: diff_minus)
  3094   from order_less_le_trans[OF `0 < ly`[unfolded less_float_def] this]
  3095   show ?thesis by auto
  3096 qed
  3097 
  3098 lemma approx_tse_form'_le:
  3099   assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\<lambda> l u. 0 \<le> l)"
  3100   and x: "x \<in> {real l .. real u}"
  3101   shows "interpret_floatarith b [x] \<le> interpret_floatarith a [x]"
  3102 proof -
  3103   from approx_tse_form'[OF tse x]
  3104   obtain l' u' ly uy
  3105     where x': "x \<in> { real l' .. real u' }" and "real l \<le> real l'"
  3106     and "real u' \<le> real u" and "0 \<le> ly"
  3107     and tse: "approx_tse prec 0 t ((l' + u') * Float 1 -1) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
  3108     by blast
  3109 
  3110   hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
  3111 
  3112   from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
  3113   have "real ly \<le> interpret_floatarith a [x] - interpret_floatarith b [x]"
  3114     by (auto simp add: diff_minus)
  3115   from order_trans[OF `0 \<le> ly`[unfolded le_float_def] this]
  3116   show ?thesis by auto
  3117 qed
  3118 
  3119 definition
  3120 "approx_tse_form prec t s f =
  3121   (case f
  3122    of (Bound x a b f) \<Rightarrow> x = Var 0 \<and>
  3123      (case (approx prec a [None], approx prec b [None])
  3124       of (Some (l, u), Some (l', u')) \<Rightarrow>
  3125         (case f
  3126          of Less lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)
  3127           | LessEqual lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)
  3128           | AtLeastAtMost x lf rt \<Rightarrow>
  3129             (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l) then
  3130             approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l) else False)
  3131           | _ \<Rightarrow> False)
  3132        | _ \<Rightarrow> False)
  3133    | _ \<Rightarrow> False)"
  3134 
  3135 lemma approx_tse_form:
  3136   assumes "approx_tse_form prec t s f"
  3137   shows "interpret_form f [x]"
  3138 proof (cases f)
  3139   case (Bound i a b f') note f_def = this
  3140   with assms obtain l u l' u'
  3141     where a: "approx prec a [None] = Some (l, u)"
  3142     and b: "approx prec b [None] = Some (l', u')"
  3143     unfolding approx_tse_form_def by (auto elim!: option_caseE)
  3144 
  3145   from Bound assms have "i = Var 0" unfolding approx_tse_form_def by auto
  3146   hence i: "interpret_floatarith i [x] = x" by auto
  3147 
  3148   { let "?f z" = "interpret_floatarith z [x]"
  3149     assume "?f i \<in> { ?f a .. ?f b }"
  3150     with approx[OF _ a[symmetric], of "[x]"] approx[OF _ b[symmetric], of "[x]"]
  3151     have bnd: "x \<in> { real l .. real u'}" unfolding bounded_by_def i by auto
  3152 
  3153     have "interpret_form f' [x]"
  3154     proof (cases f')
  3155       case (Less lf rt)
  3156       with Bound a b assms
  3157       have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)"
  3158         unfolding approx_tse_form_def by auto
  3159       from approx_tse_form'_less[OF this bnd]
  3160       show ?thesis using Less by auto
  3161     next
  3162       case (LessEqual lf rt)
  3163       with Bound a b assms
  3164       have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
  3165         unfolding approx_tse_form_def by auto
  3166       from approx_tse_form'_le[OF this bnd]
  3167       show ?thesis using LessEqual by auto
  3168     next
  3169       case (AtLeastAtMost x lf rt)
  3170       with Bound a b assms
  3171       have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l)"
  3172         and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
  3173         unfolding approx_tse_form_def lazy_conj by auto
  3174       from approx_tse_form'_le[OF this(1) bnd] approx_tse_form'_le[OF this(2) bnd]
  3175       show ?thesis using AtLeastAtMost by auto
  3176     next
  3177       case (Bound x a b f') with assms
  3178       show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
  3179     next
  3180       case (Assign x a f') with assms
  3181       show ?thesis by (auto elim!: option_caseE simp add: f_def approx_tse_form_def)
  3182     qed } thus ?thesis unfolding f_def by auto
  3183 next case Assign with assms show ?thesis by (auto simp add: approx_tse_form_def)
  3184 next case LessEqual with assms show ?thesis by (auto simp add: approx_tse_form_def)
  3185 next case Less with assms show ?thesis by (auto simp add: approx_tse_form_def)
  3186 next case AtLeastAtMost with assms show ?thesis by (auto simp add: approx_tse_form_def)
  3187 qed
  3188 
  3189 text {* @{term approx_form_eval} is only used for the {\tt value}-command. *}
  3190 
  3191 fun approx_form_eval :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option list" where
  3192 "approx_form_eval prec (Bound (Var n) a b f) bs =
  3193    (case (approx prec a bs, approx prec b bs)
  3194    of (Some (l, _), Some (_, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
  3195     | _ \<Rightarrow> bs)" |
  3196 "approx_form_eval prec (Assign (Var n) a f) bs =
  3197    (case (approx prec a bs)
  3198    of (Some (l, u)) \<Rightarrow> approx_form_eval prec f (bs[n := Some (l, u)])
  3199     | _ \<Rightarrow> bs)" |
  3200 "approx_form_eval prec (Less a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
  3201 "approx_form_eval prec (LessEqual a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
  3202 "approx_form_eval prec (AtLeastAtMost x a b) bs =
  3203    bs @ [approx prec x bs, approx prec a bs, approx prec b bs]" |
  3204 "approx_form_eval _ _ bs = bs"
  3205 
  3206 subsection {* Implement proof method \texttt{approximation} *}
  3207 
  3208 lemmas interpret_form_equations = interpret_form.simps interpret_floatarith.simps interpret_floatarith_num
  3209   interpret_floatarith_divide interpret_floatarith_diff interpret_floatarith_tan interpret_floatarith_powr interpret_floatarith_log
  3210   interpret_floatarith_sin
  3211 
  3212 ML {*
  3213 structure Float_Arith =
  3214 struct
  3215 
  3216 @{code_datatype float = Float}
  3217 @{code_datatype floatarith = Add | Minus | Mult | Inverse | Cos | Arctan
  3218                            | Abs | Max | Min | Pi | Sqrt | Exp | Ln | Power | Var | Num }
  3219 @{code_datatype form = Bound | Assign | Less | LessEqual | AtLeastAtMost}
  3220 
  3221 val approx_form = @{code approx_form}
  3222 val approx_tse_form = @{code approx_tse_form}
  3223 val approx' = @{code approx'}
  3224 val approx_form_eval = @{code approx_form_eval}
  3225 
  3226 end
  3227 *}
  3228 
  3229 code_reserved Eval Float_Arith
  3230 
  3231 code_type float (Eval "Float'_Arith.float")
  3232 code_const Float (Eval "Float'_Arith.Float/ (_,/ _)")
  3233 
  3234 code_type floatarith (Eval "Float'_Arith.floatarith")
  3235 code_const Add and Minus and Mult and Inverse and Cos and Arctan and Abs and Max and Min and
  3236            Pi and Sqrt  and Exp and Ln and Power and Var and Num
  3237   (Eval "Float'_Arith.Add/ (_,/ _)" and "Float'_Arith.Minus" and "Float'_Arith.Mult/ (_,/ _)" and
  3238         "Float'_Arith.Inverse" and "Float'_Arith.Cos" and
  3239         "Float'_Arith.Arctan" and "Float'_Arith.Abs" and "Float'_Arith.Max/ (_,/ _)" and
  3240         "Float'_Arith.Min/ (_,/ _)" and "Float'_Arith.Pi" and "Float'_Arith.Sqrt" and
  3241         "Float'_Arith.Exp" and "Float'_Arith.Ln" and "Float'_Arith.Power/ (_,/ _)" and
  3242         "Float'_Arith.Var" and "Float'_Arith.Num")
  3243 
  3244 code_type form (Eval "Float'_Arith.form")
  3245 code_const Bound and Assign and Less and LessEqual and AtLeastAtMost
  3246       (Eval "Float'_Arith.Bound/ (_,/ _,/ _,/ _)" and "Float'_Arith.Assign/ (_,/ _,/ _)" and
  3247             "Float'_Arith.Less/ (_,/ _)" and "Float'_Arith.LessEqual/ (_,/ _)"  and
  3248             "Float'_Arith.AtLeastAtMost/ (_,/ _,/ _)")
  3249 
  3250 code_const approx_form (Eval "Float'_Arith.approx'_form")
  3251 code_const approx_tse_form (Eval "Float'_Arith.approx'_tse'_form")
  3252 code_const approx' (Eval "Float'_Arith.approx'")
  3253 
  3254 ML {*
  3255   fun reorder_bounds_tac prems i =
  3256     let
  3257       fun variable_of_bound (Const ("Trueprop", _) $
  3258                              (Const (@{const_name "op :"}, _) $
  3259                               Free (name, _) $ _)) = name
  3260         | variable_of_bound (Const ("Trueprop", _) $
  3261                              (Const ("op =", _) $
  3262                               Free (name, _) $ _)) = name
  3263         | variable_of_bound t = raise TERM ("variable_of_bound", [t])
  3264 
  3265       val variable_bounds
  3266         = map (` (variable_of_bound o prop_of)) prems
  3267 
  3268       fun add_deps (name, bnds)
  3269         = Graph.add_deps_acyclic (name,
  3270             remove (op =) name (Term.add_free_names (prop_of bnds) []))
  3271 
  3272       val order = Graph.empty
  3273                   |> fold Graph.new_node variable_bounds
  3274                   |> fold add_deps variable_bounds
  3275                   |> Graph.strong_conn |> map the_single |> rev
  3276                   |> map_filter (AList.lookup (op =) variable_bounds)
  3277 
  3278       fun prepend_prem th tac
  3279         = tac THEN rtac (th RSN (2, @{thm mp})) i
  3280     in
  3281       fold prepend_prem order all_tac
  3282     end
  3283 
  3284   (* Should be in HOL.thy ? *)
  3285   fun gen_eval_tac conv ctxt = CONVERSION (Conv.params_conv (~1) (K (Conv.concl_conv (~1) conv)) ctxt)
  3286                                THEN' rtac TrueI
  3287 
  3288   val form_equations = PureThy.get_thms @{theory} "interpret_form_equations";
  3289 
  3290   fun rewrite_interpret_form_tac ctxt prec splitting taylor i st = let
  3291       fun lookup_splitting (Free (name, typ))
  3292         = case AList.lookup (op =) splitting name
  3293           of SOME s => HOLogic.mk_number @{typ nat} s
  3294            | NONE => @{term "0 :: nat"}
  3295       val vs = nth (prems_of st) (i - 1)
  3296                |> Logic.strip_imp_concl
  3297                |> HOLogic.dest_Trueprop
  3298                |> Term.strip_comb |> snd |> List.last
  3299                |> HOLogic.dest_list
  3300       val p = prec
  3301               |> HOLogic.mk_number @{typ nat}
  3302               |> Thm.cterm_of (ProofContext.theory_of ctxt)
  3303     in case taylor
  3304     of NONE => let
  3305          val n = vs |> length
  3306                  |> HOLogic.mk_number @{typ nat}
  3307                  |> Thm.cterm_of (ProofContext.theory_of ctxt)
  3308          val s = vs
  3309                  |> map lookup_splitting
  3310                  |> HOLogic.mk_list @{typ nat}
  3311                  |> Thm.cterm_of (ProofContext.theory_of ctxt)
  3312        in
  3313          (rtac (Thm.instantiate ([], [(@{cpat "?n::nat"}, n),
  3314                                      (@{cpat "?prec::nat"}, p),
  3315                                      (@{cpat "?ss::nat list"}, s)])
  3316               @{thm "approx_form"}) i
  3317           THEN simp_tac @{simpset} i) st
  3318        end
  3319 
  3320      | SOME t => if length vs <> 1 then raise (TERM ("More than one variable used for taylor series expansion", [prop_of st]))
  3321        else let
  3322          val t = t
  3323               |> HOLogic.mk_number @{typ nat}
  3324               |> Thm.cterm_of (ProofContext.theory_of ctxt)
  3325          val s = vs |> map lookup_splitting |> hd
  3326               |> Thm.cterm_of (ProofContext.theory_of ctxt)
  3327        in
  3328          rtac (Thm.instantiate ([], [(@{cpat "?s::nat"}, s),
  3329                                      (@{cpat "?t::nat"}, t),
  3330                                      (@{cpat "?prec::nat"}, p)])
  3331               @{thm "approx_tse_form"}) i st
  3332        end
  3333     end
  3334 
  3335   (* copied from Tools/induct.ML should probably in args.ML *)
  3336   val free = Args.context -- Args.term >> (fn (_, Free (n, t)) => n | (ctxt, t) =>
  3337     error ("Bad free variable: " ^ Syntax.string_of_term ctxt t));
  3338 
  3339 *}
  3340 
  3341 lemma intervalE: "a \<le> x \<and> x \<le> b \<Longrightarrow> \<lbrakk> x \<in> { a .. b } \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
  3342   by auto
  3343 
  3344 lemma meta_eqE: "x \<equiv> a \<Longrightarrow> \<lbrakk> x = a \<Longrightarrow> P\<rbrakk> \<Longrightarrow> P"
  3345   by auto
  3346 
  3347 method_setup approximation = {*
  3348   Scan.lift (OuterParse.nat)
  3349   --
  3350   Scan.optional (Scan.lift (Args.$$$ "splitting" |-- Args.colon)
  3351     |-- OuterParse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift OuterParse.nat)) []
  3352   --
  3353   Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon)
  3354     |-- (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift OuterParse.nat))
  3355   >>
  3356   (fn ((prec, splitting), taylor) => fn ctxt =>
  3357     SIMPLE_METHOD' (fn i =>
  3358       REPEAT (FIRST' [etac @{thm intervalE},
  3359                       etac @{thm meta_eqE},
  3360                       rtac @{thm impI}] i)
  3361       THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems i) @{context} i
  3362       THEN DETERM (TRY (filter_prems_tac (K false) i))
  3363       THEN DETERM (Reflection.genreify_tac ctxt form_equations NONE i)
  3364       THEN rewrite_interpret_form_tac ctxt prec splitting taylor i
  3365       THEN gen_eval_tac eval_oracle ctxt i))
  3366  *} "real number approximation"
  3367 
  3368 ML {*
  3369   fun calculated_subterms (@{const Trueprop} $ t) = calculated_subterms t
  3370     | calculated_subterms (@{const "op -->"} $ _ $ t) = calculated_subterms t
  3371     | calculated_subterms (@{term "op <= :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
  3372     | calculated_subterms (@{term "op < :: real \<Rightarrow> real \<Rightarrow> bool"} $ t1 $ t2) = [t1, t2]
  3373     | calculated_subterms (@{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ t1 $ 
  3374                            (@{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ t2 $ t3)) = [t1, t2, t3]
  3375     | calculated_subterms t = raise TERM ("calculated_subterms", [t])
  3376 
  3377   fun dest_interpret_form (@{const "interpret_form"} $ b $ xs) = (b, xs)
  3378     | dest_interpret_form t = raise TERM ("dest_interpret_form", [t])
  3379 
  3380   fun dest_interpret (@{const "interpret_floatarith"} $ b $ xs) = (b, xs)
  3381     | dest_interpret t = raise TERM ("dest_interpret", [t])
  3382 
  3383 
  3384   fun dest_float (@{const "Float"} $ m $ e) = (snd (HOLogic.dest_number m), snd (HOLogic.dest_number e))
  3385   fun dest_ivl (Const (@{const_name "Some"}, _) $
  3386                 (Const (@{const_name "Pair"}, _) $ u $ l)) = SOME (dest_float u, dest_float l)
  3387     | dest_ivl (Const (@{const_name "None"}, _)) = NONE
  3388     | dest_ivl t = raise TERM ("dest_result", [t])
  3389 
  3390   fun mk_approx' prec t = (@{const "approx'"}
  3391                          $ HOLogic.mk_number @{typ nat} prec
  3392                          $ t $ @{term "[] :: (float * float) option list"})
  3393 
  3394   fun mk_approx_form_eval prec t xs = (@{const "approx_form_eval"}
  3395                          $ HOLogic.mk_number @{typ nat} prec
  3396                          $ t $ xs)
  3397 
  3398   fun float2_float10 prec round_down (m, e) = (
  3399     let
  3400       val (m, e) = (if e < 0 then (m,e) else (m * Integer.pow e 2, 0))
  3401 
  3402       fun frac c p 0 digits cnt = (digits, cnt, 0)
  3403         | frac c 0 r digits cnt = (digits, cnt, r)
  3404         | frac c p r digits cnt = (let
  3405           val (d, r) = Integer.div_mod (r * 10) (Integer.pow (~e) 2)
  3406         in frac (c orelse d <> 0) (if d <> 0 orelse c then p - 1 else p) r
  3407                 (digits * 10 + d) (cnt + 1)
  3408         end)
  3409 
  3410       val sgn = Int.sign m
  3411       val m = abs m
  3412 
  3413       val round_down = (sgn = 1 andalso round_down) orelse
  3414                        (sgn = ~1 andalso not round_down)
  3415 
  3416       val (x, r) = Integer.div_mod m (Integer.pow (~e) 2)
  3417 
  3418       val p = ((if x = 0 then prec else prec - (IntInf.log2 x + 1)) * 3) div 10 + 1
  3419 
  3420       val (digits, e10, r) = if p > 0 then frac (x <> 0) p r 0 0 else (0,0,0)
  3421 
  3422       val digits = if round_down orelse r = 0 then digits else digits + 1
  3423 
  3424     in (sgn * (digits + x * (Integer.pow e10 10)), ~e10)
  3425     end)
  3426 
  3427   fun mk_result prec (SOME (l, u)) = (let
  3428       fun mk_float10 rnd x = (let val (m, e) = float2_float10 prec rnd x
  3429                          in if e = 0 then HOLogic.mk_number @{typ real} m
  3430                        else if e = 1 then @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
  3431                                           HOLogic.mk_number @{typ real} m $
  3432                                           @{term "10"}
  3433                                      else @{term "divide :: real \<Rightarrow> real \<Rightarrow> real"} $
  3434                                           HOLogic.mk_number @{typ real} m $
  3435                                           (@{term "power 10 :: nat \<Rightarrow> real"} $
  3436                                            HOLogic.mk_number @{typ nat} (~e)) end)
  3437       in @{term "atLeastAtMost :: real \<Rightarrow> real \<Rightarrow> real set"} $ mk_float10 true l $ mk_float10 false u end)
  3438     | mk_result prec NONE = @{term "UNIV :: real set"}
  3439 
  3440   fun realify t = let
  3441       val t = Logic.varify_global t
  3442       val m = map (fn (name, sort) => (name, @{typ real})) (Term.add_tvars t [])
  3443       val t = Term.subst_TVars m t
  3444     in t end
  3445 
  3446   fun converted_result t =
  3447           prop_of t
  3448        |> HOLogic.dest_Trueprop
  3449        |> HOLogic.dest_eq |> snd
  3450 
  3451   fun apply_tactic context term tactic = cterm_of context term
  3452     |> Goal.init
  3453     |> SINGLE tactic
  3454     |> the |> prems_of |> hd
  3455 
  3456   fun prepare_form context term = apply_tactic context term (
  3457       REPEAT (FIRST' [etac @{thm intervalE}, etac @{thm meta_eqE}, rtac @{thm impI}] 1)
  3458       THEN Subgoal.FOCUS (fn {prems, ...} => reorder_bounds_tac prems 1) @{context} 1
  3459       THEN DETERM (TRY (filter_prems_tac (K false) 1)))
  3460 
  3461   fun reify_form context term = apply_tactic context term
  3462      (Reflection.genreify_tac @{context} form_equations NONE 1)
  3463 
  3464   fun approx_form prec ctxt t =
  3465           realify t
  3466        |> prepare_form (ProofContext.theory_of ctxt)
  3467        |> (fn arith_term =>
  3468           reify_form (ProofContext.theory_of ctxt) arith_term
  3469        |> HOLogic.dest_Trueprop |> dest_interpret_form
  3470        |> (fn (data, xs) =>
  3471           mk_approx_form_eval prec data (HOLogic.mk_list @{typ "(float * float) option"}
  3472             (map (fn _ => @{term "None :: (float * float) option"}) (HOLogic.dest_list xs)))
  3473        |> Codegen.eval_term @{theory}
  3474        |> HOLogic.dest_list
  3475        |> curry ListPair.zip (HOLogic.dest_list xs @ calculated_subterms arith_term)
  3476        |> map (fn (elem, s) => @{term "op : :: real \<Rightarrow> real set \<Rightarrow> bool"} $ elem $ mk_result prec (dest_ivl s))
  3477        |> foldr1 HOLogic.mk_conj))
  3478 
  3479   fun approx_arith prec ctxt t = realify t
  3480        |> Reflection.genreif ctxt form_equations
  3481        |> prop_of
  3482        |> HOLogic.dest_Trueprop
  3483        |> HOLogic.dest_eq |> snd
  3484        |> dest_interpret |> fst
  3485        |> mk_approx' prec
  3486        |> Codegen.eval_term @{theory}
  3487        |> dest_ivl
  3488        |> mk_result prec
  3489 
  3490    fun approx prec ctxt t = if type_of t = @{typ prop} then approx_form prec ctxt t
  3491      else if type_of t = @{typ bool} then approx_form prec ctxt (@{const Trueprop} $ t)
  3492      else approx_arith prec ctxt t
  3493 *}
  3494 
  3495 setup {*
  3496   Value.add_evaluator ("approximate", approx 30)
  3497 *}
  3498 
  3499 end