Made Approximation work for powr again
authorManuel Eberl <eberlm@in.tum.de>
Tue Jan 19 07:59:29 2016 +0100 (2016-01-19 ago)
changeset 6220067792e4a5486
parent 62199 fc55a4e3f439
child 62201 eca7b38c8ee5
child 62202 e5bc7cbb0bcc
Made Approximation work for powr again
src/HOL/Decision_Procs/Approximation.thy
     1.1 --- a/src/HOL/Decision_Procs/Approximation.thy	Mon Jan 18 16:03:58 2016 +0100
     1.2 +++ b/src/HOL/Decision_Procs/Approximation.thy	Tue Jan 19 07:59:29 2016 +0100
     1.3 @@ -191,6 +191,62 @@
     1.4      l1 \<le> x ^ n \<and> x ^ n \<le> u1"
     1.5    using float_power_bnds by auto
     1.6  
     1.7 +section \<open>Approximation utility functions\<close>
     1.8 +
     1.9 +definition bnds_mult :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<times> float" where
    1.10 +  "bnds_mult prec a1 a2 b1 b2 =
    1.11 +      (float_plus_down prec (nprt a1 * pprt b2)
    1.12 +          (float_plus_down prec (nprt a2 * nprt b2)
    1.13 +            (float_plus_down prec (pprt a1 * pprt b1) (pprt a2 * nprt b1))),
    1.14 +        float_plus_up prec (pprt a2 * pprt b2)
    1.15 +            (float_plus_up prec (pprt a1 * nprt b2)
    1.16 +              (float_plus_up prec (nprt a2 * pprt b1) (nprt a1 * nprt b1))))"
    1.17 +
    1.18 +lemma bnds_mult:
    1.19 +  fixes prec :: nat and a1 aa2 b1 b2 :: float
    1.20 +  assumes "(l, u) = bnds_mult prec a1 a2 b1 b2"
    1.21 +  assumes "a \<in> {real_of_float a1..real_of_float a2}"
    1.22 +  assumes "b \<in> {real_of_float b1..real_of_float b2}"
    1.23 +  shows   "a * b \<in> {real_of_float l..real_of_float u}"
    1.24 +proof -
    1.25 +  from assms have "real_of_float l \<le> a * b" 
    1.26 +    by (intro order.trans[OF _ mult_ge_prts[of a1 a a2 b1 b b2]])
    1.27 +       (auto simp: bnds_mult_def intro!: float_plus_down_le)
    1.28 +  moreover from assms have "real_of_float u \<ge> a * b" 
    1.29 +    by (intro order.trans[OF mult_le_prts[of a1 a a2 b1 b b2]])
    1.30 +       (auto simp: bnds_mult_def intro!: float_plus_up_le)
    1.31 +  ultimately show ?thesis by simp
    1.32 +qed
    1.33 +
    1.34 +definition map_bnds :: "(nat \<Rightarrow> float \<Rightarrow> float) \<Rightarrow> (nat \<Rightarrow> float \<Rightarrow> float) \<Rightarrow>
    1.35 +                           nat \<Rightarrow> (float \<times> float) \<Rightarrow> (float \<times> float)" where
    1.36 +  "map_bnds lb ub prec = (\<lambda>(l,u). (lb prec l, ub prec u))"
    1.37 +
    1.38 +lemma map_bnds:
    1.39 +  assumes "(lf, uf) = map_bnds lb ub prec (l, u)"
    1.40 +  assumes "mono f"
    1.41 +  assumes "x \<in> {real_of_float l..real_of_float u}"
    1.42 +  assumes "real_of_float (lb prec l) \<le> f (real_of_float l)"
    1.43 +  assumes "real_of_float (ub prec u) \<ge> f (real_of_float u)"
    1.44 +  shows   "f x \<in> {real_of_float lf..real_of_float uf}"
    1.45 +proof -
    1.46 +  from assms have "real_of_float lf = real_of_float (lb prec l)"
    1.47 +    by (simp add: map_bnds_def)
    1.48 +  also have "real_of_float (lb prec l) \<le> f (real_of_float l)"  by fact
    1.49 +  also from assms have "\<dots> \<le> f x"
    1.50 +    by (intro monoD[OF \<open>mono f\<close>]) auto
    1.51 +  finally have lf: "real_of_float lf \<le> f x" .
    1.52 +
    1.53 +  from assms have "f x \<le> f (real_of_float u)"
    1.54 +    by (intro monoD[OF \<open>mono f\<close>]) auto
    1.55 +  also have "\<dots> \<le> real_of_float (ub prec u)" by fact
    1.56 +  also from assms have "\<dots> = real_of_float uf"
    1.57 +    by (simp add: map_bnds_def)
    1.58 +  finally have uf: "f x \<le> real_of_float uf" .
    1.59 +
    1.60 +  from lf uf show ?thesis by simp
    1.61 +qed
    1.62 +
    1.63  
    1.64  section "Square root"
    1.65  
    1.66 @@ -2648,6 +2704,101 @@
    1.67  qed
    1.68  
    1.69  
    1.70 +section \<open>Real power function\<close>
    1.71 +
    1.72 +definition bnds_powr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float \<times> float) option" where
    1.73 +  "bnds_powr prec l1 u1 l2 u2 = (
    1.74 +     if l1 = 0 \<and> u1 = 0 then
    1.75 +       Some (0, 0)
    1.76 +     else if l1 = 0 \<and> l2 \<ge> 1 then
    1.77 +       let uln = the (ub_ln prec u1)
    1.78 +       in  Some (0, ub_exp prec (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2))))
    1.79 +     else if l1 \<le> 0 then
    1.80 +       None
    1.81 +     else
    1.82 +       Some (map_bnds lb_exp ub_exp prec 
    1.83 +               (bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2)))"
    1.84 +
    1.85 +lemmas [simp del] = lb_exp.simps ub_exp.simps
    1.86 +
    1.87 +lemma mono_exp_real: "mono (exp :: real \<Rightarrow> real)"
    1.88 +  by (auto simp: mono_def)
    1.89 +
    1.90 +lemma ub_exp_nonneg: "real_of_float (ub_exp prec x) \<ge> 0"
    1.91 +proof -
    1.92 +  have "0 \<le> exp (real_of_float x)" by simp
    1.93 +  also from exp_boundaries[of x prec] 
    1.94 +    have "\<dots> \<le> real_of_float (ub_exp prec x)" by simp
    1.95 +  finally show ?thesis .
    1.96 +qed
    1.97 +
    1.98 +lemma bnds_powr:
    1.99 +  assumes lu: "Some (l, u) = bnds_powr prec l1 u1 l2 u2"
   1.100 +  assumes x: "x \<in> {real_of_float l1..real_of_float u1}"
   1.101 +  assumes y: "y \<in> {real_of_float l2..real_of_float u2}"
   1.102 +  shows   "x powr y \<in> {real_of_float l..real_of_float u}"
   1.103 +proof -
   1.104 +  consider "l1 = 0" "u1 = 0" | "l1 = 0" "u1 \<noteq> 0" "l2 \<ge> 1" | 
   1.105 +           "l1 \<le> 0" "\<not>(l1 = 0 \<and> (u1 = 0 \<or> l2 \<ge> 1))" | "l1 > 0" by force
   1.106 +  thus ?thesis
   1.107 +  proof cases
   1.108 +    assume "l1 = 0" "u1 = 0"
   1.109 +    with x lu show ?thesis by (auto simp: bnds_powr_def)
   1.110 +  next
   1.111 +    assume A: "l1 = 0" "u1 \<noteq> 0" "l2 \<ge> 1"
   1.112 +    def uln \<equiv> "the (ub_ln prec u1)"
   1.113 +    show ?thesis
   1.114 +    proof (cases "x = 0")
   1.115 +      case False
   1.116 +      with A x y have "x powr y = exp (ln x * y)" by (simp add: powr_def)
   1.117 +      also {
   1.118 +        from A x False have "ln x \<le> ln (real_of_float u1)" by simp
   1.119 +        also from ub_ln_lb_ln_bounds[of u1 prec] A y x False
   1.120 +          have "ln (real_of_float u1) \<le> real_of_float uln" by (simp add: uln_def del: lb_ln.simps)
   1.121 +        also from A x y have "\<dots> * y \<le> real_of_float uln * (if uln \<ge> 0 then u2 else l2)"
   1.122 +          by (auto intro: mult_left_mono mult_left_mono_neg)
   1.123 +        also have "\<dots> \<le> real_of_float (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2)))"
   1.124 +          by (simp add: float_round_up_le)
   1.125 +        finally have "ln x * y \<le> \<dots>" using A y by - simp
   1.126 +      }
   1.127 +      also have "exp (real_of_float (float_round_up prec (uln * (if uln \<ge> 0 then u2 else l2)))) \<le>
   1.128 +                   real_of_float (ub_exp prec (float_round_up prec
   1.129 +                       (uln * (if uln \<ge> 0 then u2 else l2))))"
   1.130 +        using exp_boundaries by simp
   1.131 +      finally show ?thesis using A x y lu 
   1.132 +        by (simp add: bnds_powr_def uln_def Let_def del: lb_ln.simps ub_ln.simps)
   1.133 +    qed (insert x y lu A, simp_all add: bnds_powr_def Let_def ub_exp_nonneg
   1.134 +                                   del: lb_ln.simps ub_ln.simps)
   1.135 +  next
   1.136 +    assume "l1 \<le> 0" "\<not>(l1 = 0 \<and> (u1 = 0 \<or> l2 \<ge> 1))"
   1.137 +    with lu show ?thesis by (simp add: bnds_powr_def split: split_if_asm)
   1.138 +  next
   1.139 +    assume l1: "l1 > 0"
   1.140 +    obtain lm um where lmum:
   1.141 +      "(lm, um) = bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2"
   1.142 +      by (cases "bnds_mult prec (the (lb_ln prec l1)) (the (ub_ln prec u1)) l2 u2") simp
   1.143 +    with l1 have "(l, u) = map_bnds lb_exp ub_exp prec (lm, um)"
   1.144 +      using lu by (simp add: bnds_powr_def del: lb_ln.simps ub_ln.simps split: split_if_asm)
   1.145 +    hence "exp (ln x * y) \<in> {real_of_float l..real_of_float u}"
   1.146 +    proof (rule map_bnds[OF _ mono_exp_real], goal_cases)
   1.147 +      case 1
   1.148 +      let ?lln = "the (lb_ln prec l1)" and ?uln = "the (ub_ln prec u1)"
   1.149 +      from ub_ln_lb_ln_bounds[of l1 prec] ub_ln_lb_ln_bounds[of u1 prec] x l1
   1.150 +        have "real_of_float ?lln \<le> ln (real_of_float l1) \<and> 
   1.151 +              ln (real_of_float u1) \<le> real_of_float ?uln"
   1.152 +        by (auto simp del: lb_ln.simps ub_ln.simps)
   1.153 +      moreover from l1 x have "ln (real_of_float l1) \<le> ln x \<and> ln x \<le> ln (real_of_float u1)"
   1.154 +        by auto
   1.155 +      ultimately have ln: "real_of_float ?lln \<le> ln x \<and> ln x \<le> real_of_float ?uln" by simp
   1.156 +      from lmum show ?case
   1.157 +        by (rule bnds_mult) (insert y ln, simp_all)
   1.158 +    qed (insert exp_boundaries[of lm prec] exp_boundaries[of um prec], simp_all)
   1.159 +    with x l1 show ?thesis
   1.160 +      by (simp add: powr_def mult_ac)
   1.161 +  qed
   1.162 +qed
   1.163 +
   1.164 +
   1.165  section "Implement floatarith"
   1.166  
   1.167  subsection "Define syntax and semantics"
   1.168 @@ -2665,6 +2816,7 @@
   1.169    | Pi
   1.170    | Sqrt floatarith
   1.171    | Exp floatarith
   1.172 +  | Powr floatarith floatarith
   1.173    | Ln floatarith
   1.174    | Power floatarith nat
   1.175    | Var nat
   1.176 @@ -2683,6 +2835,7 @@
   1.177  "interpret_floatarith Pi vs           = pi" |
   1.178  "interpret_floatarith (Sqrt a) vs     = sqrt (interpret_floatarith a vs)" |
   1.179  "interpret_floatarith (Exp a) vs      = exp (interpret_floatarith a vs)" |
   1.180 +"interpret_floatarith (Powr a b) vs   = interpret_floatarith a vs powr interpret_floatarith b vs" |
   1.181  "interpret_floatarith (Ln a) vs       = ln (interpret_floatarith a vs)" |
   1.182  "interpret_floatarith (Power a n) vs  = (interpret_floatarith a vs)^n" |
   1.183  "interpret_floatarith (Num f) vs      = f" |
   1.184 @@ -2727,6 +2880,10 @@
   1.185  
   1.186  subsection "Implement approximation function"
   1.187  
   1.188 +fun lift_bin :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float) option) \<Rightarrow> (float * float) option" where
   1.189 +"lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2" |
   1.190 +"lift_bin a b f = None"
   1.191 +
   1.192  fun lift_bin' :: "(float * float) option \<Rightarrow> (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float)) \<Rightarrow> (float * float) option" where
   1.193  "lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
   1.194  "lift_bin' a b f = None"
   1.195 @@ -2787,14 +2944,7 @@
   1.196      (\<lambda> l1 u1 l2 u2. (float_plus_down prec l1 l2, float_plus_up prec u1 u2))" |
   1.197  "approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
   1.198  "approx prec (Mult a b) bs  =
   1.199 -  lift_bin' (approx' prec a bs) (approx' prec b bs)
   1.200 -    (\<lambda> a1 a2 b1 b2.
   1.201 -      (float_plus_down prec (nprt a1 * pprt b2)
   1.202 -          (float_plus_down prec (nprt a2 * nprt b2)
   1.203 -            (float_plus_down prec (pprt a1 * pprt b1) (pprt a2 * nprt b1))),
   1.204 -        float_plus_up prec (pprt a2 * pprt b2)
   1.205 -            (float_plus_up prec (pprt a1 * nprt b2)
   1.206 -              (float_plus_up prec (nprt a2 * pprt b1) (nprt a1 * nprt b1)))))" |
   1.207 +  lift_bin' (approx' prec a bs) (approx' prec b bs) (bnds_mult prec)" |
   1.208  "approx prec (Inverse a) bs = lift_un (approx' prec a bs) (\<lambda> l u. if (0 < l \<or> u < 0) then (Some (float_divl prec 1 u), Some (float_divr prec 1 l)) else (None, None))" |
   1.209  "approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |
   1.210  "approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |
   1.211 @@ -2804,11 +2954,92 @@
   1.212  "approx prec (Arctan a) bs  = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_arctan prec l, ub_arctan prec u))" |
   1.213  "approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
   1.214  "approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
   1.215 +"approx prec (Powr a b) bs  = lift_bin (approx' prec a bs) (approx' prec b bs) (bnds_powr prec)" |
   1.216  "approx prec (Ln a) bs      = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
   1.217  "approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds prec n)" |
   1.218  "approx prec (Num f) bs     = Some (f, f)" |
   1.219  "approx prec (Var i) bs    = (if i < length bs then bs ! i else None)"
   1.220  
   1.221 +lemma approx_approx':
   1.222 +  assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow>
   1.223 +      l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
   1.224 +    and approx': "Some (l, u) = approx' prec a vs"
   1.225 +  shows "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
   1.226 +proof -
   1.227 +  obtain l' u' where S: "Some (l', u') = approx prec a vs"
   1.228 +    using approx' unfolding approx'.simps by (cases "approx prec a vs") auto
   1.229 +  have l': "l = float_round_down prec l'" and u': "u = float_round_up prec u'"
   1.230 +    using approx' unfolding approx'.simps S[symmetric] by auto
   1.231 +  show ?thesis unfolding l' u'
   1.232 +    using order_trans[OF Pa[OF S, THEN conjunct2] float_round_up[of u']]
   1.233 +    using order_trans[OF float_round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
   1.234 +qed
   1.235 +
   1.236 +lemma lift_bin_ex:
   1.237 +  assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
   1.238 +  shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
   1.239 +proof (cases a)
   1.240 +  case None
   1.241 +  hence "None = lift_bin a b f"
   1.242 +    unfolding None lift_bin.simps ..
   1.243 +  thus ?thesis
   1.244 +    using lift_bin_Some by auto
   1.245 +next
   1.246 +  case (Some a')
   1.247 +  show ?thesis
   1.248 +  proof (cases b)
   1.249 +    case None
   1.250 +    hence "None = lift_bin a b f"
   1.251 +      unfolding None lift_bin.simps ..
   1.252 +    thus ?thesis using lift_bin_Some by auto
   1.253 +  next
   1.254 +    case (Some b')
   1.255 +    obtain la ua where a': "a' = (la, ua)"
   1.256 +      by (cases a') auto
   1.257 +    obtain lb ub where b': "b' = (lb, ub)"
   1.258 +      by (cases b') auto
   1.259 +    thus ?thesis
   1.260 +      unfolding \<open>a = Some a'\<close> \<open>b = Some b'\<close> a' b' by auto
   1.261 +  qed
   1.262 +qed
   1.263 +
   1.264 +lemma lift_bin_f:
   1.265 +  assumes lift_bin_Some: "Some (l, u) = lift_bin (g a) (g b) f"
   1.266 +    and Pa: "\<And>l u. Some (l, u) = g a \<Longrightarrow> P l u a"
   1.267 +    and Pb: "\<And>l u. Some (l, u) = g b \<Longrightarrow> P l u b"
   1.268 +  shows "\<exists> l1 u1 l2 u2. P l1 u1 a \<and> P l2 u2 b \<and> Some (l, u) = f l1 u1 l2 u2"
   1.269 +proof -
   1.270 +  obtain l1 u1 l2 u2
   1.271 +    where Sa: "Some (l1, u1) = g a"
   1.272 +      and Sb: "Some (l2, u2) = g b"
   1.273 +    using lift_bin_ex[OF assms(1)] by auto
   1.274 +  have lu: "Some (l, u) = f l1 u1 l2 u2"
   1.275 +    using lift_bin_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin.simps] by auto
   1.276 +  thus ?thesis
   1.277 +    using Pa[OF Sa] Pb[OF Sb] by auto
   1.278 +qed
   1.279 +
   1.280 +lemma lift_bin:
   1.281 +  assumes lift_bin_Some: "Some (l, u) = lift_bin (approx' prec a bs) (approx' prec b bs) f"
   1.282 +    and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow>
   1.283 +      real_of_float l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real_of_float u" (is "\<And>l u. _ = ?g a \<Longrightarrow> ?P l u a")
   1.284 +    and Pb: "\<And>l u. Some (l, u) = approx prec b bs \<Longrightarrow>
   1.285 +      real_of_float l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real_of_float u"
   1.286 +  shows "\<exists>l1 u1 l2 u2. (real_of_float l1 \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> real_of_float u1) \<and>
   1.287 +                       (real_of_float l2 \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> real_of_float u2) \<and>
   1.288 +                       Some (l, u) = (f l1 u1 l2 u2)"
   1.289 +proof -
   1.290 +  { fix l u assume "Some (l, u) = approx' prec a bs"
   1.291 +    with approx_approx'[of prec a bs, OF _ this] Pa
   1.292 +    have "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u" by auto } note Pa = this
   1.293 +  { fix l u assume "Some (l, u) = approx' prec b bs"
   1.294 +    with approx_approx'[of prec b bs, OF _ this] Pb
   1.295 +    have "l \<le> interpret_floatarith b xs \<and> interpret_floatarith b xs \<le> u" by auto } note Pb = this
   1.296 +
   1.297 +  from lift_bin_f[where g="\<lambda>a. approx' prec a bs" and P = ?P, OF lift_bin_Some, OF Pa Pb]
   1.298 +  show ?thesis by auto
   1.299 +qed
   1.300 +
   1.301  lemma lift_bin'_ex:
   1.302    assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
   1.303    shows "\<exists> l1 u1 l2 u2. Some (l1, u1) = a \<and> Some (l2, u2) = b"
   1.304 @@ -2855,21 +3086,6 @@
   1.305      using Pa[OF Sa] Pb[OF Sb] by auto
   1.306  qed
   1.307  
   1.308 -lemma approx_approx':
   1.309 -  assumes Pa: "\<And>l u. Some (l, u) = approx prec a vs \<Longrightarrow>
   1.310 -      l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
   1.311 -    and approx': "Some (l, u) = approx' prec a vs"
   1.312 -  shows "l \<le> interpret_floatarith a xs \<and> interpret_floatarith a xs \<le> u"
   1.313 -proof -
   1.314 -  obtain l' u' where S: "Some (l', u') = approx prec a vs"
   1.315 -    using approx' unfolding approx'.simps by (cases "approx prec a vs") auto
   1.316 -  have l': "l = float_round_down prec l'" and u': "u = float_round_up prec u'"
   1.317 -    using approx' unfolding approx'.simps S[symmetric] by auto
   1.318 -  show ?thesis unfolding l' u'
   1.319 -    using order_trans[OF Pa[OF S, THEN conjunct2] float_round_up[of u']]
   1.320 -    using order_trans[OF float_round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
   1.321 -qed
   1.322 -
   1.323  lemma lift_bin':
   1.324    assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
   1.325      and Pa: "\<And>l u. Some (l, u) = approx prec a bs \<Longrightarrow>
   1.326 @@ -3075,21 +3291,12 @@
   1.327    case (Mult a b)
   1.328    from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
   1.329    obtain l1 u1 l2 u2
   1.330 -    where l: "l = float_plus_down prec (nprt l1 * pprt u2) (float_plus_down prec (nprt u1 * nprt u2) (float_plus_down prec (pprt l1 * pprt l2) (pprt u1 * nprt l2)))"
   1.331 -    and u: "u = float_plus_up prec (pprt u1 * pprt u2) (float_plus_up prec (pprt l1 * nprt u2) (float_plus_up prec (nprt u1 * pprt l2) (nprt l1 * nprt l2)))"
   1.332 -    and "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1"
   1.333 -    and "l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
   1.334 -  hence bnds:
   1.335 -    "nprt l1 * pprt u2 + nprt u1 * nprt u2 + pprt l1 * pprt l2 + pprt u1 * nprt l2 \<le>
   1.336 -      interpret_floatarith (Mult a b) xs" (is "?l \<le> _")
   1.337 -    "interpret_floatarith (Mult a b) xs \<le>
   1.338 -      pprt u1 * pprt u2 + pprt l1 * nprt u2 + nprt u1 * pprt l2 + nprt l1 * nprt l2" (is "_ \<le> ?u")
   1.339 -    unfolding interpret_floatarith.simps l u
   1.340 -    using mult_le_prts mult_ge_prts by auto
   1.341 -  from l u have "l \<le> ?l" "?u \<le> u"
   1.342 -    by (auto intro!: float_plus_up_le float_plus_down_le)
   1.343 -  thus ?case
   1.344 -    using bnds by simp
   1.345 +    where l: "l = fst (bnds_mult prec l1 u1 l2 u2)"
   1.346 +    and u: "u = snd (bnds_mult prec l1 u1 l2 u2)"
   1.347 +    and a: "l1 \<le> interpret_floatarith a xs" "interpret_floatarith a xs \<le> u1"
   1.348 +    and b: "l2 \<le> interpret_floatarith b xs" "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
   1.349 +  from l u have lu: "(l, u) = bnds_mult prec l1 u1 l2 u2" by simp
   1.350 +  from bnds_mult[OF lu] a b show ?case by simp
   1.351  next
   1.352    case (Inverse a)
   1.353    from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
   1.354 @@ -3194,6 +3401,15 @@
   1.355    case (Exp a)
   1.356    with lift_un'_bnds[OF bnds_exp] show ?case by auto
   1.357  next
   1.358 +  case (Powr a b)
   1.359 +  from lift_bin[OF Powr.prems[unfolded approx.simps]] Powr.hyps
   1.360 +    obtain l1 u1 l2 u2 where lu: "Some (l, u) = bnds_powr prec l1 u1 l2 u2"
   1.361 +      and l1: "l1 \<le> interpret_floatarith a xs" and u1: "interpret_floatarith a xs \<le> u1"
   1.362 +      and l2: "l2 \<le> interpret_floatarith b xs" and u2: "interpret_floatarith b xs \<le> u2"
   1.363 +      by blast
   1.364 +  from bnds_powr[OF lu] l1 u1 l2 u2
   1.365 +    show ?case by simp
   1.366 +next
   1.367    case (Ln a)
   1.368    with lift_un_bnds[OF bnds_ln] show ?case by auto
   1.369  next
   1.370 @@ -3402,6 +3618,8 @@
   1.371  "isDERIV x Pi vs                = True" |
   1.372  "isDERIV x (Sqrt a) vs          = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
   1.373  "isDERIV x (Exp a) vs           = isDERIV x a vs" |
   1.374 +"isDERIV x (Powr a b) vs        = 
   1.375 +    (isDERIV x a vs \<and> isDERIV x b vs \<and> interpret_floatarith a vs > 0)" |
   1.376  "isDERIV x (Ln a) vs            = (isDERIV x a vs \<and> interpret_floatarith a vs > 0)" |
   1.377  "isDERIV x (Power a 0) vs       = True" |
   1.378  "isDERIV x (Power a (Suc n)) vs = isDERIV x a vs" |
   1.379 @@ -3421,12 +3639,34 @@
   1.380  "DERIV_floatarith x Pi                = Num 0" |
   1.381  "DERIV_floatarith x (Sqrt a)          = (Mult (Inverse (Mult (Sqrt a) (Num 2))) (DERIV_floatarith x a))" |
   1.382  "DERIV_floatarith x (Exp a)           = Mult (Exp a) (DERIV_floatarith x a)" |
   1.383 +"DERIV_floatarith x (Powr a b)        =
   1.384 +   Mult (Powr a b) (Add (Mult (DERIV_floatarith x b) (Ln a)) (Mult (Mult (DERIV_floatarith x a) b) (Inverse a)))" |
   1.385  "DERIV_floatarith x (Ln a)            = Mult (Inverse a) (DERIV_floatarith x a)" |
   1.386  "DERIV_floatarith x (Power a 0)       = Num 0" |
   1.387  "DERIV_floatarith x (Power a (Suc n)) = Mult (Num (Float (int (Suc n)) 0)) (Mult (Power a n) (DERIV_floatarith x a))" |
   1.388  "DERIV_floatarith x (Num f)           = Num 0" |
   1.389  "DERIV_floatarith x (Var n)          = (if x = n then Num 1 else Num 0)"
   1.390  
   1.391 +lemma has_real_derivative_powr':
   1.392 +  fixes f g :: "real \<Rightarrow> real"
   1.393 +  assumes "(f has_real_derivative f') (at x)"
   1.394 +  assumes "(g has_real_derivative g') (at x)"
   1.395 +  assumes "f x > 0"
   1.396 +  defines "h \<equiv> \<lambda>x. f x powr g x * (g' * ln (f x) + f' * g x / f x)"
   1.397 +  shows   "((\<lambda>x. f x powr g x) has_real_derivative h x) (at x)"
   1.398 +proof (subst DERIV_cong_ev[OF refl _ refl])
   1.399 +  from assms have "isCont f x"
   1.400 +    by (simp add: DERIV_continuous)
   1.401 +  hence "f \<midarrow>x\<rightarrow> f x" by (simp add: continuous_at)
   1.402 +  with \<open>f x > 0\<close> have "eventually (\<lambda>x. f x > 0) (nhds x)"
   1.403 +    by (auto simp: tendsto_at_iff_tendsto_nhds dest: order_tendstoD)
   1.404 +  thus "eventually (\<lambda>x. f x powr g x = exp (g x * ln (f x))) (nhds x)"
   1.405 +    by eventually_elim (simp add: powr_def)
   1.406 +next
   1.407 +  from assms show "((\<lambda>x. exp (g x * ln (f x))) has_real_derivative h x) (at x)"
   1.408 +    by (auto intro!: derivative_eq_intros simp: h_def powr_def)
   1.409 +qed
   1.410 +
   1.411  lemma DERIV_floatarith:
   1.412    assumes "n < length vs"
   1.413    assumes isDERIV: "isDERIV n f (vs[n := x])"
   1.414 @@ -3454,6 +3694,11 @@
   1.415  next
   1.416    case (Var i)
   1.417    thus ?case using \<open>n < length vs\<close> by auto
   1.418 +next
   1.419 +  case (Powr a b)
   1.420 +  note [derivative_intros] = has_real_derivative_powr'
   1.421 +  from Powr show ?case
   1.422 +    by (auto intro!: derivative_eq_intros simp: field_simps)
   1.423  qed (auto intro!: derivative_eq_intros)
   1.424  
   1.425  declare approx.simps[simp del]
   1.426 @@ -3473,12 +3718,14 @@
   1.427  "isDERIV_approx prec x (Sqrt a) vs          =
   1.428    (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
   1.429  "isDERIV_approx prec x (Exp a) vs           = isDERIV_approx prec x a vs" |
   1.430 +"isDERIV_approx prec x (Powr a b) vs        =
   1.431 +  (isDERIV_approx prec x a vs \<and> isDERIV_approx prec x b vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
   1.432  "isDERIV_approx prec x (Ln a) vs            =
   1.433    (isDERIV_approx prec x a vs \<and> (case approx prec a vs of Some (l, u) \<Rightarrow> 0 < l | None \<Rightarrow> False))" |
   1.434  "isDERIV_approx prec x (Power a 0) vs       = True" |
   1.435  "isDERIV_approx prec x (Power a (Suc n)) vs = isDERIV_approx prec x a vs" |
   1.436  "isDERIV_approx prec x (Num f) vs           = True" |
   1.437 -"isDERIV_approx prec x (Var n) vs          = True"
   1.438 +"isDERIV_approx prec x (Var n) vs           = True"
   1.439  
   1.440  lemma isDERIV_approx:
   1.441    assumes "bounded_by xs vs"
   1.442 @@ -3512,6 +3759,13 @@
   1.443  next
   1.444    case (Power a n)
   1.445    thus ?case by (cases n) auto
   1.446 +next
   1.447 +  case (Powr a b)
   1.448 +  from Powr obtain l1 u1 where a: "Some (l1, u1) = approx prec a vs" and pos: "0 < l1"
   1.449 +    by (cases "approx prec a vs") auto
   1.450 +  with approx[OF \<open>bounded_by xs vs\<close> a]
   1.451 +    have "0 < interpret_floatarith a xs" by auto
   1.452 +  with Powr show ?case by auto
   1.453  qed auto
   1.454  
   1.455  lemma bounded_by_update_var:
   1.456 @@ -3578,13 +3832,7 @@
   1.457      by simp
   1.458  qed
   1.459  
   1.460 -fun lift_bin :: "(float * float) option \<Rightarrow>
   1.461 -    (float * float) option \<Rightarrow> (float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> (float * float) option) \<Rightarrow>
   1.462 -    (float * float) option" where
   1.463 -  "lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2"
   1.464 -| "lift_bin a b f = None"
   1.465 -
   1.466 -lemma lift_bin:
   1.467 +lemma lift_bin_aux:
   1.468    assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
   1.469    obtains l1 u1 l2 u2
   1.470    where "a = Some (l1, u1)"
   1.471 @@ -3592,6 +3840,7 @@
   1.472      and "f l1 u1 l2 u2 = Some (l, u)"
   1.473    using assms by (cases a, simp, cases b, simp, auto)
   1.474  
   1.475 +
   1.476  fun approx_tse where
   1.477  "approx_tse prec n 0 c k f bs = approx prec f bs" |
   1.478  "approx_tse prec n (Suc s) c k f bs =
   1.479 @@ -3676,7 +3925,7 @@
   1.480                 (Mult (Inverse (Num (Float (int k) 0)))
   1.481                       (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
   1.482                             (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), vs!x]"
   1.483 -      by (auto elim!: lift_bin)
   1.484 +      by (auto elim!: lift_bin_aux)
   1.485  
   1.486      from bnd_c \<open>x < length xs\<close>
   1.487      have bnd: "bounded_by (xs[x:=c]) (vs[x:= Some (c,c)])"
   1.488 @@ -4052,6 +4301,7 @@
   1.489      | floatarith_of_term @{term Pi} = @{code Pi}
   1.490      | floatarith_of_term (@{term Sqrt} $ a) = @{code Sqrt} (floatarith_of_term a)
   1.491      | floatarith_of_term (@{term Exp} $ a) = @{code Exp} (floatarith_of_term a)
   1.492 +    | floatarith_of_term (@{term Powr} $ a $ b) = @{code Powr} (floatarith_of_term a, floatarith_of_term b)
   1.493      | floatarith_of_term (@{term Ln} $ a) = @{code Ln} (floatarith_of_term a)
   1.494      | floatarith_of_term (@{term Power} $ a $ n) =
   1.495          @{code Power} (floatarith_of_term a, nat_of_term n)