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