author Manuel Eberl Tue Jan 19 07:59:29 2016 +0100 (2016-01-19 ago) changeset 62200 67792e4a5486 parent 62199 fc55a4e3f439 child 62201 eca7b38c8ee5 child 62202 e5bc7cbb0bcc
Made Approximation work for powr again
```     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)
```