Theory Approximation

theory Approximation
imports Code_Target_Numeral Approximation_Bounds
 (* Author:     Johannes Hoelzl, TU Muenchen
   Coercions removed by Dmitriy Traytel *)

theory Approximation
imports
  Complex_Main
  "HOL-Library.Code_Target_Numeral"
  Approximation_Bounds
keywords "approximate" :: diag
begin

section "Implement floatarith"

subsection "Define syntax and semantics"

datatype floatarith
  = Add floatarith floatarith
  | Minus floatarith
  | Mult floatarith floatarith
  | Inverse floatarith
  | Cos floatarith
  | Arctan floatarith
  | Abs floatarith
  | Max floatarith floatarith
  | Min floatarith floatarith
  | Pi
  | Sqrt floatarith
  | Exp floatarith
  | Powr floatarith floatarith
  | Ln floatarith
  | Power floatarith nat
  | Floor floatarith
  | Var nat
  | Num float

fun interpret_floatarith :: "floatarith ⇒ real list ⇒ real" where
"interpret_floatarith (Add a b) vs   = (interpret_floatarith a vs) + (interpret_floatarith b vs)" |
"interpret_floatarith (Minus a) vs    = - (interpret_floatarith a vs)" |
"interpret_floatarith (Mult a b) vs   = (interpret_floatarith a vs) * (interpret_floatarith b vs)" |
"interpret_floatarith (Inverse a) vs  = inverse (interpret_floatarith a vs)" |
"interpret_floatarith (Cos a) vs      = cos (interpret_floatarith a vs)" |
"interpret_floatarith (Arctan a) vs   = arctan (interpret_floatarith a vs)" |
"interpret_floatarith (Min a b) vs    = min (interpret_floatarith a vs) (interpret_floatarith b vs)" |
"interpret_floatarith (Max a b) vs    = max (interpret_floatarith a vs) (interpret_floatarith b vs)" |
"interpret_floatarith (Abs a) vs      = ¦interpret_floatarith a vs¦" |
"interpret_floatarith Pi vs           = pi" |
"interpret_floatarith (Sqrt a) vs     = sqrt (interpret_floatarith a vs)" |
"interpret_floatarith (Exp a) vs      = exp (interpret_floatarith a vs)" |
"interpret_floatarith (Powr a b) vs   = interpret_floatarith a vs powr interpret_floatarith b vs" |
"interpret_floatarith (Ln a) vs       = ln (interpret_floatarith a vs)" |
"interpret_floatarith (Power a n) vs  = (interpret_floatarith a vs)^n" |
"interpret_floatarith (Floor a) vs      = floor (interpret_floatarith a vs)" |
"interpret_floatarith (Num f) vs      = f" |
"interpret_floatarith (Var n) vs     = vs ! n"

lemma interpret_floatarith_divide:
  "interpret_floatarith (Mult a (Inverse b)) vs =
    (interpret_floatarith a vs) / (interpret_floatarith b vs)"
  unfolding divide_inverse interpret_floatarith.simps ..

lemma interpret_floatarith_diff:
  "interpret_floatarith (Add a (Minus b)) vs =
    (interpret_floatarith a vs) - (interpret_floatarith b vs)"
  unfolding interpret_floatarith.simps by simp

lemma interpret_floatarith_sin:
  "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 (- 1)))) (Minus a))) vs =
    sin (interpret_floatarith a vs)"
  unfolding sin_cos_eq interpret_floatarith.simps
    interpret_floatarith_divide interpret_floatarith_diff
  by auto


subsection "Implement approximation function"

fun lift_bin :: "(float * float) option ⇒ (float * float) option ⇒ (float ⇒ float ⇒ float ⇒ float ⇒ (float * float) option) ⇒ (float * float) option" where
"lift_bin (Some (l1, u1)) (Some (l2, u2)) f = f l1 u1 l2 u2" |
"lift_bin a b f = None"

fun lift_bin' :: "(float * float) option ⇒ (float * float) option ⇒ (float ⇒ float ⇒ float ⇒ float ⇒ (float * float)) ⇒ (float * float) option" where
"lift_bin' (Some (l1, u1)) (Some (l2, u2)) f = Some (f l1 u1 l2 u2)" |
"lift_bin' a b f = None"

fun lift_un :: "(float * float) option ⇒ (float ⇒ float ⇒ ((float option) * (float option))) ⇒ (float * float) option" where
"lift_un (Some (l1, u1)) f = (case (f l1 u1) of (Some l, Some u) ⇒ Some (l, u)
                                             | t ⇒ None)" |
"lift_un b f = None"

fun lift_un' :: "(float * float) option ⇒ (float ⇒ float ⇒ (float * float)) ⇒ (float * float) option" where
"lift_un' (Some (l1, u1)) f = Some (f l1 u1)" |
"lift_un' b f = None"

definition bounded_by :: "real list ⇒ (float × float) option list ⇒ bool" where 
  "bounded_by xs vs ⟷
  (∀ i < length vs. case vs ! i of None ⇒ True
         | Some (l, u) ⇒ xs ! i ∈ { real_of_float l .. real_of_float u })"
                                                                     
lemma bounded_byE:
  assumes "bounded_by xs vs"
  shows "⋀ i. i < length vs ⟹ case vs ! i of None ⇒ True
         | Some (l, u) ⇒ xs ! i ∈ { real_of_float l .. real_of_float u }"
  using assms bounded_by_def by blast

lemma bounded_by_update:
  assumes "bounded_by xs vs"
    and bnd: "xs ! i ∈ { real_of_float l .. real_of_float u }"
  shows "bounded_by xs (vs[i := Some (l,u)])"
proof -
  {
    fix j
    let ?vs = "vs[i := Some (l,u)]"
    assume "j < length ?vs"
    hence [simp]: "j < length vs" by simp
    have "case ?vs ! j of None ⇒ True | Some (l, u) ⇒ xs ! j ∈ { real_of_float l .. real_of_float u }"
    proof (cases "?vs ! j")
      case (Some b)
      thus ?thesis
      proof (cases "i = j")
        case True
        thus ?thesis using ‹?vs ! j = Some b› and bnd by auto
      next
        case False
        thus ?thesis using ‹bounded_by xs vs› unfolding bounded_by_def by auto
      qed
    qed auto
  }
  thus ?thesis unfolding bounded_by_def by auto
qed

lemma bounded_by_None: "bounded_by xs (replicate (length xs) None)"
  unfolding bounded_by_def by auto

fun approx approx' :: "nat ⇒ floatarith ⇒ (float * float) option list ⇒ (float * float) option" where
"approx' prec a bs          = (case (approx prec a bs) of Some (l, u) ⇒ Some (float_round_down prec l, float_round_up prec u) | None ⇒ None)" |
"approx prec (Add a b) bs   =
  lift_bin' (approx' prec a bs) (approx' prec b bs)
    (λ l1 u1 l2 u2. (float_plus_down prec l1 l2, float_plus_up prec u1 u2))" |
"approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (λ l u. (-u, -l))" |
"approx prec (Mult a b) bs  =
  lift_bin' (approx' prec a bs) (approx' prec b bs) (bnds_mult prec)" |
"approx prec (Inverse a) bs = lift_un (approx' prec a bs) (λ l u. if (0 < l ∨ u < 0) then (Some (float_divl prec 1 u), Some (float_divr prec 1 l)) else (None, None))" |
"approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |
"approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |
"approx prec (Min a b) bs   = lift_bin' (approx' prec a bs) (approx' prec b bs) (λ l1 u1 l2 u2. (min l1 l2, min u1 u2))" |
"approx prec (Max a b) bs   = lift_bin' (approx' prec a bs) (approx' prec b bs) (λ l1 u1 l2 u2. (max l1 l2, max u1 u2))" |
"approx prec (Abs a) bs     = lift_un' (approx' prec a bs) (λl u. (if l < 0 ∧ 0 < u then 0 else min ¦l¦ ¦u¦, max ¦l¦ ¦u¦))" |
"approx prec (Arctan a) bs  = lift_un' (approx' prec a bs) (λ l u. (lb_arctan prec l, ub_arctan prec u))" |
"approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (λ l u. (lb_sqrt prec l, ub_sqrt prec u))" |
"approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (λ l u. (lb_exp prec l, ub_exp prec u))" |
"approx prec (Powr a b) bs  = lift_bin (approx' prec a bs) (approx' prec b bs) (bnds_powr prec)" |
"approx prec (Ln a) bs      = lift_un (approx' prec a bs) (λ l u. (lb_ln prec l, ub_ln prec u))" |
"approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds prec n)" |
"approx prec (Floor a) bs = lift_un' (approx' prec a bs) (λ l u. (floor_fl l, floor_fl u))" |
"approx prec (Num f) bs     = Some (f, f)" |
"approx prec (Var i) bs    = (if i < length bs then bs ! i else None)"

lemma approx_approx':
  assumes Pa: "⋀l u. Some (l, u) = approx prec a vs ⟹
      l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
    and approx': "Some (l, u) = approx' prec a vs"
  shows "l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
proof -
  obtain l' u' where S: "Some (l', u') = approx prec a vs"
    using approx' unfolding approx'.simps by (cases "approx prec a vs") auto
  have l': "l = float_round_down prec l'" and u': "u = float_round_up prec u'"
    using approx' unfolding approx'.simps S[symmetric] by auto
  show ?thesis unfolding l' u'
    using order_trans[OF Pa[OF S, THEN conjunct2] float_round_up[of u']]
    using order_trans[OF float_round_down[of _ l'] Pa[OF S, THEN conjunct1]] by auto
qed

lemma lift_bin_ex:
  assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
  shows "∃ l1 u1 l2 u2. Some (l1, u1) = a ∧ Some (l2, u2) = b"
proof (cases a)
  case None
  hence "None = lift_bin a b f"
    unfolding None lift_bin.simps ..
  thus ?thesis
    using lift_bin_Some by auto
next
  case (Some a')
  show ?thesis
  proof (cases b)
    case None
    hence "None = lift_bin a b f"
      unfolding None lift_bin.simps ..
    thus ?thesis using lift_bin_Some by auto
  next
    case (Some b')
    obtain la ua where a': "a' = (la, ua)"
      by (cases a') auto
    obtain lb ub where b': "b' = (lb, ub)"
      by (cases b') auto
    thus ?thesis
      unfolding ‹a = Some a'› ‹b = Some b'› a' b' by auto
  qed
qed

lemma lift_bin_f:
  assumes lift_bin_Some: "Some (l, u) = lift_bin (g a) (g b) f"
    and Pa: "⋀l u. Some (l, u) = g a ⟹ P l u a"
    and Pb: "⋀l u. Some (l, u) = g b ⟹ P l u b"
  shows "∃ l1 u1 l2 u2. P l1 u1 a ∧ P l2 u2 b ∧ Some (l, u) = f l1 u1 l2 u2"
proof -
  obtain l1 u1 l2 u2
    where Sa: "Some (l1, u1) = g a"
      and Sb: "Some (l2, u2) = g b"
    using lift_bin_ex[OF assms(1)] by auto
  have lu: "Some (l, u) = f l1 u1 l2 u2"
    using lift_bin_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin.simps] by auto
  thus ?thesis
    using Pa[OF Sa] Pb[OF Sb] by auto
qed

lemma lift_bin:
  assumes lift_bin_Some: "Some (l, u) = lift_bin (approx' prec a bs) (approx' prec b bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
      real_of_float l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ real_of_float u" (is "⋀l u. _ = ?g a ⟹ ?P l u a")
    and Pb: "⋀l u. Some (l, u) = approx prec b bs ⟹
      real_of_float l ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ real_of_float u"
  shows "∃l1 u1 l2 u2. (real_of_float l1 ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ real_of_float u1) ∧
                       (real_of_float l2 ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ real_of_float u2) ∧
                       Some (l, u) = (f l1 u1 l2 u2)"
proof -
  { fix l u assume "Some (l, u) = approx' prec a bs"
    with approx_approx'[of prec a bs, OF _ this] Pa
    have "l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u" by auto } note Pa = this
  { fix l u assume "Some (l, u) = approx' prec b bs"
    with approx_approx'[of prec b bs, OF _ this] Pb
    have "l ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ u" by auto } note Pb = this

  from lift_bin_f[where g="λa. approx' prec a bs" and P = ?P, OF lift_bin_Some, OF Pa Pb]
  show ?thesis by auto
qed

lemma lift_bin'_ex:
  assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f"
  shows "∃ l1 u1 l2 u2. Some (l1, u1) = a ∧ Some (l2, u2) = b"
proof (cases a)
  case None
  hence "None = lift_bin' a b f"
    unfolding None lift_bin'.simps ..
  thus ?thesis
    using lift_bin'_Some by auto
next
  case (Some a')
  show ?thesis
  proof (cases b)
    case None
    hence "None = lift_bin' a b f"
      unfolding None lift_bin'.simps ..
    thus ?thesis using lift_bin'_Some by auto
  next
    case (Some b')
    obtain la ua where a': "a' = (la, ua)"
      by (cases a') auto
    obtain lb ub where b': "b' = (lb, ub)"
      by (cases b') auto
    thus ?thesis
      unfolding ‹a = Some a'› ‹b = Some b'› a' b' by auto
  qed
qed

lemma lift_bin'_f:
  assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f"
    and Pa: "⋀l u. Some (l, u) = g a ⟹ P l u a"
    and Pb: "⋀l u. Some (l, u) = g b ⟹ P l u b"
  shows "∃ l1 u1 l2 u2. P l1 u1 a ∧ P l2 u2 b ∧ l = fst (f l1 u1 l2 u2) ∧ u = snd (f l1 u1 l2 u2)"
proof -
  obtain l1 u1 l2 u2
    where Sa: "Some (l1, u1) = g a"
      and Sb: "Some (l2, u2) = g b"
    using lift_bin'_ex[OF assms(1)] by auto
  have lu: "(l, u) = f l1 u1 l2 u2"
    using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto
  have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)"
    unfolding lu[symmetric] by auto
  thus ?thesis
    using Pa[OF Sa] Pb[OF Sb] by auto
qed

lemma lift_bin':
  assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
      l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u" (is "⋀l u. _ = ?g a ⟹ ?P l u a")
    and Pb: "⋀l u. Some (l, u) = approx prec b bs ⟹
      l ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ u"
  shows "∃l1 u1 l2 u2. (l1 ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u1) ∧
                       (l2 ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ u2) ∧
                       l = fst (f l1 u1 l2 u2) ∧ u = snd (f l1 u1 l2 u2)"
proof -
  { fix l u assume "Some (l, u) = approx' prec a bs"
    with approx_approx'[of prec a bs, OF _ this] Pa
    have "l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u" by auto } note Pa = this
  { fix l u assume "Some (l, u) = approx' prec b bs"
    with approx_approx'[of prec b bs, OF _ this] Pb
    have "l ≤ interpret_floatarith b xs ∧ interpret_floatarith b xs ≤ u" by auto } note Pb = this

  from lift_bin'_f[where g="λa. approx' prec a bs" and P = ?P, OF lift_bin'_Some, OF Pa Pb]
  show ?thesis by auto
qed

lemma lift_un'_ex:
  assumes lift_un'_Some: "Some (l, u) = lift_un' a f"
  shows "∃ l u. Some (l, u) = a"
proof (cases a)
  case None
  hence "None = lift_un' a f"
    unfolding None lift_un'.simps ..
  thus ?thesis
    using lift_un'_Some by auto
next
  case (Some a')
  obtain la ua where a': "a' = (la, ua)"
    by (cases a') auto
  thus ?thesis
    unfolding ‹a = Some a'› a' by auto
qed

lemma lift_un'_f:
  assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f"
    and Pa: "⋀l u. Some (l, u) = g a ⟹ P l u a"
  shows "∃ l1 u1. P l1 u1 a ∧ l = fst (f l1 u1) ∧ u = snd (f l1 u1)"
proof -
  obtain l1 u1 where Sa: "Some (l1, u1) = g a"
    using lift_un'_ex[OF assms(1)] by auto
  have lu: "(l, u) = f l1 u1"
    using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto
  have "l = fst (f l1 u1)" and "u = snd (f l1 u1)"
    unfolding lu[symmetric] by auto
  thus ?thesis
    using Pa[OF Sa] by auto
qed

lemma lift_un':
  assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
      l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
      (is "⋀l u. _ = ?g a ⟹ ?P l u a")
  shows "∃l1 u1. (l1 ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u1) ∧
    l = fst (f l1 u1) ∧ u = snd (f l1 u1)"
proof -
  have Pa: "l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
    if "Some (l, u) = approx' prec a bs" for l u
    using approx_approx'[of prec a bs, OF _ that] Pa
     by auto
  from lift_un'_f[where g="λa. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa]
  show ?thesis by auto
qed

lemma lift_un'_bnds:
  assumes bnds: "∀ (x::real) lx ux. (l, u) = f lx ux ∧ x ∈ { lx .. ux } ⟶ l ≤ f' x ∧ f' x ≤ u"
    and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
      l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
  shows "real_of_float l ≤ f' (interpret_floatarith a xs) ∧ f' (interpret_floatarith a xs) ≤ real_of_float u"
proof -
  from lift_un'[OF lift_un'_Some Pa]
  obtain l1 u1 where "l1 ≤ interpret_floatarith a xs"
    and "interpret_floatarith a xs ≤ u1"
    and "l = fst (f l1 u1)"
    and "u = snd (f l1 u1)"
    by blast
  hence "(l, u) = f l1 u1" and "interpret_floatarith a xs ∈ {l1 .. u1}"
    by auto
  thus ?thesis
    using bnds by auto
qed

lemma lift_un_ex:
  assumes lift_un_Some: "Some (l, u) = lift_un a f"
  shows "∃l u. Some (l, u) = a"
proof (cases a)
  case None
  hence "None = lift_un a f"
    unfolding None lift_un.simps ..
  thus ?thesis
    using lift_un_Some by auto
next
  case (Some a')
  obtain la ua where a': "a' = (la, ua)"
    by (cases a') auto
  thus ?thesis
    unfolding ‹a = Some a'› a' by auto
qed

lemma lift_un_f:
  assumes lift_un_Some: "Some (l, u) = lift_un (g a) f"
    and Pa: "⋀l u. Some (l, u) = g a ⟹ P l u a"
  shows "∃ l1 u1. P l1 u1 a ∧ Some l = fst (f l1 u1) ∧ Some u = snd (f l1 u1)"
proof -
  obtain l1 u1 where Sa: "Some (l1, u1) = g a"
    using lift_un_ex[OF assms(1)] by auto
  have "fst (f l1 u1) ≠ None ∧ snd (f l1 u1) ≠ None"
  proof (rule ccontr)
    assume "¬ (fst (f l1 u1) ≠ None ∧ snd (f l1 u1) ≠ None)"
    hence or: "fst (f l1 u1) = None ∨ snd (f l1 u1) = None" by auto
    hence "lift_un (g a) f = None"
    proof (cases "fst (f l1 u1) = None")
      case True
      then obtain b where b: "f l1 u1 = (None, b)"
        by (cases "f l1 u1") auto
      thus ?thesis
        unfolding Sa[symmetric] lift_un.simps b by auto
    next
      case False
      hence "snd (f l1 u1) = None"
        using or by auto
      with False obtain b where b: "f l1 u1 = (Some b, None)"
        by (cases "f l1 u1") auto
      thus ?thesis
        unfolding Sa[symmetric] lift_un.simps b by auto
    qed
    thus False
      using lift_un_Some by auto
  qed
  then obtain a' b' where f: "f l1 u1 = (Some a', Some b')"
    by (cases "f l1 u1") auto
  from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f]
  have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)"
    unfolding f by auto
  thus ?thesis
    unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto
qed

lemma lift_un:
  assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
        l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
      (is "⋀l u. _ = ?g a ⟹ ?P l u a")
  shows "∃l1 u1. (l1 ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u1) ∧
                  Some l = fst (f l1 u1) ∧ Some u = snd (f l1 u1)"
proof -
  have Pa: "l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
    if "Some (l, u) = approx' prec a bs" for l u
    using approx_approx'[of prec a bs, OF _ that] Pa by auto
  from lift_un_f[where g="λa. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa]
  show ?thesis by auto
qed

lemma lift_un_bnds:
  assumes bnds: "∀(x::real) lx ux. (Some l, Some u) = f lx ux ∧ x ∈ { lx .. ux } ⟶ l ≤ f' x ∧ f' x ≤ u"
    and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f"
    and Pa: "⋀l u. Some (l, u) = approx prec a bs ⟹
      l ≤ interpret_floatarith a xs ∧ interpret_floatarith a xs ≤ u"
  shows "real_of_float l ≤ f' (interpret_floatarith a xs) ∧ f' (interpret_floatarith a xs) ≤ real_of_float u"
proof -
  from lift_un[OF lift_un_Some Pa]
  obtain l1 u1 where "l1 ≤ interpret_floatarith a xs"
    and "interpret_floatarith a xs ≤ u1"
    and "Some l = fst (f l1 u1)"
    and "Some u = snd (f l1 u1)"
    by blast
  hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs ∈ {l1 .. u1}"
    by auto
  thus ?thesis
    using bnds by auto
qed

lemma approx:
  assumes "bounded_by xs vs"
    and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith")
  shows "l ≤ interpret_floatarith arith xs ∧ interpret_floatarith arith xs ≤ u" (is "?P l u arith")
  using ‹Some (l, u) = approx prec arith vs›
proof (induct arith arbitrary: l u)
  case (Add a b)
  from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
  obtain l1 u1 l2 u2 where "l = float_plus_down prec l1 l2"
    and "u = float_plus_up prec u1 u2" "l1 ≤ interpret_floatarith a xs"
    and "interpret_floatarith a xs ≤ u1" "l2 ≤ interpret_floatarith b xs"
    and "interpret_floatarith b xs ≤ u2"
    unfolding fst_conv snd_conv by blast
  thus ?case
    unfolding interpret_floatarith.simps by (auto intro!: float_plus_up_le float_plus_down_le)
next
  case (Minus a)
  from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
  obtain l1 u1 where "l = -u1" "u = -l1"
    and "l1 ≤ interpret_floatarith a xs" "interpret_floatarith a xs ≤ u1"
    unfolding fst_conv snd_conv by blast
  thus ?case
    unfolding interpret_floatarith.simps using minus_float.rep_eq by auto
next
  case (Mult a b)
  from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
  obtain l1 u1 l2 u2
    where l: "l = fst (bnds_mult prec l1 u1 l2 u2)"
    and u: "u = snd (bnds_mult prec l1 u1 l2 u2)"
    and a: "l1 ≤ interpret_floatarith a xs" "interpret_floatarith a xs ≤ u1"
    and b: "l2 ≤ interpret_floatarith b xs" "interpret_floatarith b xs ≤ u2" unfolding fst_conv snd_conv by blast
  from l u have lu: "(l, u) = bnds_mult prec l1 u1 l2 u2" by simp
  from bnds_mult[OF lu] a b show ?case by simp
next
  case (Inverse a)
  from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
  obtain l1 u1 where l': "Some l = (if 0 < l1 ∨ u1 < 0 then Some (float_divl prec 1 u1) else None)"
    and u': "Some u = (if 0 < l1 ∨ u1 < 0 then Some (float_divr prec 1 l1) else None)"
    and l1: "l1 ≤ interpret_floatarith a xs"
    and u1: "interpret_floatarith a xs ≤ u1"
    by blast
  have either: "0 < l1 ∨ u1 < 0"
  proof (rule ccontr)
    assume P: "¬ (0 < l1 ∨ u1 < 0)"
    show False
      using l' unfolding if_not_P[OF P] by auto
  qed
  moreover have l1_le_u1: "real_of_float l1 ≤ real_of_float u1"
    using l1 u1 by auto
  ultimately have "real_of_float l1 ≠ 0" and "real_of_float u1 ≠ 0"
    by auto

  have inv: "inverse u1 ≤ inverse (interpret_floatarith a xs)
           ∧ inverse (interpret_floatarith a xs) ≤ inverse l1"
  proof (cases "0 < l1")
    case True
    hence "0 < real_of_float u1" and "0 < real_of_float l1" "0 < interpret_floatarith a xs"
      using l1_le_u1 l1 by auto
    show ?thesis
      unfolding inverse_le_iff_le[OF ‹0 < real_of_float u1› ‹0 < interpret_floatarith a xs›]
        inverse_le_iff_le[OF ‹0 < interpret_floatarith a xs› ‹0 < real_of_float l1›]
      using l1 u1 by auto
  next
    case False
    hence "u1 < 0"
      using either by blast
    hence "real_of_float u1 < 0" and "real_of_float l1 < 0" "interpret_floatarith a xs < 0"
      using l1_le_u1 u1 by auto
    show ?thesis
      unfolding inverse_le_iff_le_neg[OF ‹real_of_float u1 < 0› ‹interpret_floatarith a xs < 0›]
        inverse_le_iff_le_neg[OF ‹interpret_floatarith a xs < 0› ‹real_of_float l1 < 0›]
      using l1 u1 by auto
  qed

  from l' have "l = float_divl prec 1 u1"
    by (cases "0 < l1 ∨ u1 < 0") auto
  hence "l ≤ inverse u1"
    unfolding nonzero_inverse_eq_divide[OF ‹real_of_float u1 ≠ 0›]
    using float_divl[of prec 1 u1] by auto
  also have "… ≤ inverse (interpret_floatarith a xs)"
    using inv by auto
  finally have "l ≤ inverse (interpret_floatarith a xs)" .
  moreover
  from u' have "u = float_divr prec 1 l1"
    by (cases "0 < l1 ∨ u1 < 0") auto
  hence "inverse l1 ≤ u"
    unfolding nonzero_inverse_eq_divide[OF ‹real_of_float l1 ≠ 0›]
    using float_divr[of 1 l1 prec] by auto
  hence "inverse (interpret_floatarith a xs) ≤ u"
    by (rule order_trans[OF inv[THEN conjunct2]])
  ultimately show ?case
    unfolding interpret_floatarith.simps using l1 u1 by auto
next
  case (Abs x)
  from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps
  obtain l1 u1 where l': "l = (if l1 < 0 ∧ 0 < u1 then 0 else min ¦l1¦ ¦u1¦)"
    and u': "u = max ¦l1¦ ¦u1¦"
    and l1: "l1 ≤ interpret_floatarith x xs"
    and u1: "interpret_floatarith x xs ≤ u1"
    by blast
  thus ?case
    unfolding l' u'
    by (cases "l1 < 0 ∧ 0 < u1") (auto simp add: real_of_float_min real_of_float_max)
next
  case (Min a b)
  from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps
  obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2"
    and l1: "l1 ≤ interpret_floatarith a xs" and u1: "interpret_floatarith a xs ≤ u1"
    and l1: "l2 ≤ interpret_floatarith b xs" and u1: "interpret_floatarith b xs ≤ u2"
    by blast
  thus ?case
    unfolding l' u' by (auto simp add: real_of_float_min)
next
  case (Max a b)
  from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps
  obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2"
    and l1: "l1 ≤ interpret_floatarith a xs" and u1: "interpret_floatarith a xs ≤ u1"
    and l1: "l2 ≤ interpret_floatarith b xs" and u1: "interpret_floatarith b xs ≤ u2"
    by blast
  thus ?case
    unfolding l' u' by (auto simp add: real_of_float_max)
next
  case (Cos a)
  with lift_un'_bnds[OF bnds_cos] show ?case by auto
next
  case (Arctan a)
  with lift_un'_bnds[OF bnds_arctan] show ?case by auto
next
  case Pi
  with pi_boundaries show ?case by auto
next
  case (Sqrt a)
  with lift_un'_bnds[OF bnds_sqrt] show ?case by auto
next
  case (Exp a)
  with lift_un'_bnds[OF bnds_exp] show ?case by auto
next
  case (Powr a b)
  from lift_bin[OF Powr.prems[unfolded approx.simps]] Powr.hyps
    obtain l1 u1 l2 u2 where lu: "Some (l, u) = bnds_powr prec l1 u1 l2 u2"
      and l1: "l1 ≤ interpret_floatarith a xs" and u1: "interpret_floatarith a xs ≤ u1"
      and l2: "l2 ≤ interpret_floatarith b xs" and u2: "interpret_floatarith b xs ≤ u2"
      by blast
  from bnds_powr[OF lu] l1 u1 l2 u2
    show ?case by simp
next
  case (Ln a)
  with lift_un_bnds[OF bnds_ln] show ?case by auto
next
  case (Power a n)
  with lift_un'_bnds[OF bnds_power] show ?case by auto
next
  case (Floor a)
  from lift_un'[OF Floor.prems[unfolded approx.simps] Floor.hyps]
  show ?case by (auto simp: floor_fl.rep_eq floor_mono)
next
  case (Num f)
  thus ?case by auto
next
  case (Var n)
  from this[symmetric] ‹bounded_by xs vs›[THEN bounded_byE, of n]
  show ?case by (cases "n < length vs") auto
qed

datatype form = Bound floatarith floatarith floatarith form
              | Assign floatarith floatarith form
              | Less floatarith floatarith
              | LessEqual floatarith floatarith
              | AtLeastAtMost floatarith floatarith floatarith
              | Conj form form
              | Disj form form

fun interpret_form :: "form ⇒ real list ⇒ bool" where
"interpret_form (Bound x a b f) vs = (interpret_floatarith x vs ∈ { interpret_floatarith a vs .. interpret_floatarith b vs } ⟶ interpret_form f vs)" |
"interpret_form (Assign x a f) vs  = (interpret_floatarith x vs = interpret_floatarith a vs ⟶ interpret_form f vs)" |
"interpret_form (Less a b) vs      = (interpret_floatarith a vs < interpret_floatarith b vs)" |
"interpret_form (LessEqual a b) vs = (interpret_floatarith a vs ≤ interpret_floatarith b vs)" |
"interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs ∈ { interpret_floatarith a vs .. interpret_floatarith b vs })" |
"interpret_form (Conj f g) vs ⟷ interpret_form f vs ∧ interpret_form g vs" |
"interpret_form (Disj f g) vs ⟷ interpret_form f vs ∨ interpret_form g vs"

fun approx_form' and approx_form :: "nat ⇒ form ⇒ (float * float) option list ⇒ nat list ⇒ bool" where
"approx_form' prec f 0 n l u bs ss = approx_form prec f (bs[n := Some (l, u)]) ss" |
"approx_form' prec f (Suc s) n l u bs ss =
  (let m = (l + u) * Float 1 (- 1)
   in (if approx_form' prec f s n l m bs ss then approx_form' prec f s n m u bs ss else False))" |
"approx_form prec (Bound (Var n) a b f) bs ss =
   (case (approx prec a bs, approx prec b bs)
   of (Some (l, _), Some (_, u)) ⇒ approx_form' prec f (ss ! n) n l u bs ss
    | _ ⇒ False)" |
"approx_form prec (Assign (Var n) a f) bs ss =
   (case (approx prec a bs)
   of (Some (l, u)) ⇒ approx_form' prec f (ss ! n) n l u bs ss
    | _ ⇒ False)" |
"approx_form prec (Less a b) bs ss =
   (case (approx prec a bs, approx prec b bs)
   of (Some (l, u), Some (l', u')) ⇒ float_plus_up prec u (-l') < 0
    | _ ⇒ False)" |
"approx_form prec (LessEqual a b) bs ss =
   (case (approx prec a bs, approx prec b bs)
   of (Some (l, u), Some (l', u')) ⇒ float_plus_up prec u (-l') ≤ 0
    | _ ⇒ False)" |
"approx_form prec (AtLeastAtMost x a b) bs ss =
   (case (approx prec x bs, approx prec a bs, approx prec b bs)
   of (Some (lx, ux), Some (l, u), Some (l', u')) ⇒ float_plus_up prec u (-lx) ≤ 0 ∧ float_plus_up prec ux (-l') ≤ 0
    | _ ⇒ False)" |
"approx_form prec (Conj a b) bs ss ⟷ approx_form prec a bs ss ∧ approx_form prec b bs ss" |
"approx_form prec (Disj a b) bs ss ⟷ approx_form prec a bs ss ∨ approx_form prec b bs ss" |
"approx_form _ _ _ _ = False"

lemma lazy_conj: "(if A then B else False) = (A ∧ B)" by simp

lemma approx_form_approx_form':
  assumes "approx_form' prec f s n l u bs ss"
    and "(x::real) ∈ { l .. u }"
  obtains l' u' where "x ∈ { l' .. u' }"
    and "approx_form prec f (bs[n := Some (l', u')]) ss"
using assms proof (induct s arbitrary: l u)
  case 0
  from this(1)[of l u] this(2,3)
  show thesis by auto
next
  case (Suc s)

  let ?m = "(l + u) * Float 1 (- 1)"
  have "real_of_float l ≤ ?m" and "?m ≤ real_of_float u"
    unfolding less_eq_float_def using Suc.prems by auto

  with ‹x ∈ { l .. u }›
  have "x ∈ { l .. ?m} ∨ x ∈ { ?m .. u }" by auto
  thus thesis
  proof (rule disjE)
    assume *: "x ∈ { l .. ?m }"
    with Suc.hyps[OF _ _ *] Suc.prems
    show thesis by (simp add: Let_def lazy_conj)
  next
    assume *: "x ∈ { ?m .. u }"
    with Suc.hyps[OF _ _ *] Suc.prems
    show thesis by (simp add: Let_def lazy_conj)
  qed
qed

lemma approx_form_aux:
  assumes "approx_form prec f vs ss"
    and "bounded_by xs vs"
  shows "interpret_form f xs"
using assms proof (induct f arbitrary: vs)
  case (Bound x a b f)
  then obtain n
    where x_eq: "x = Var n" by (cases x) auto

  with Bound.prems obtain l u' l' u
    where l_eq: "Some (l, u') = approx prec a vs"
    and u_eq: "Some (l', u) = approx prec b vs"
    and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
    by (cases "approx prec a vs", simp) (cases "approx prec b vs", auto)

  have "interpret_form f xs"
    if "xs ! n ∈ { interpret_floatarith a xs .. interpret_floatarith b xs }"
  proof -
    from approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq] that
    have "xs ! n ∈ { l .. u}" by auto

    from approx_form_approx_form'[OF approx_form' this]
    obtain lx ux where bnds: "xs ! n ∈ { lx .. ux }"
      and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .

    from ‹bounded_by xs vs› bnds have "bounded_by xs (vs[n := Some (lx, ux)])"
      by (rule bounded_by_update)
    with Bound.hyps[OF approx_form] show ?thesis
      by blast
  qed
  thus ?case
    using interpret_form.simps x_eq and interpret_floatarith.simps by simp
next
  case (Assign x a f)
  then obtain n where x_eq: "x = Var n"
    by (cases x) auto

  with Assign.prems obtain l u
    where bnd_eq: "Some (l, u) = approx prec a vs"
    and x_eq: "x = Var n"
    and approx_form': "approx_form' prec f (ss ! n) n l u vs ss"
    by (cases "approx prec a vs") auto

  have "interpret_form f xs"
    if bnds: "xs ! n = interpret_floatarith a xs"
  proof -
    from approx[OF Assign.prems(2) bnd_eq] bnds
    have "xs ! n ∈ { l .. u}" by auto
    from approx_form_approx_form'[OF approx_form' this]
    obtain lx ux where bnds: "xs ! n ∈ { lx .. ux }"
      and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" .

    from ‹bounded_by xs vs› bnds have "bounded_by xs (vs[n := Some (lx, ux)])"
      by (rule bounded_by_update)
    with Assign.hyps[OF approx_form] show ?thesis
      by blast
  qed
  thus ?case
    using interpret_form.simps x_eq and interpret_floatarith.simps by simp
next
  case (Less a b)
  then obtain l u l' u'
    where l_eq: "Some (l, u) = approx prec a vs"
      and u_eq: "Some (l', u') = approx prec b vs"
      and inequality: "real_of_float (float_plus_up prec u (-l')) < 0"
    by (cases "approx prec a vs", auto, cases "approx prec b vs", auto)
  from le_less_trans[OF float_plus_up inequality]
    approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
  show ?case by auto
next
  case (LessEqual a b)
  then obtain l u l' u'
    where l_eq: "Some (l, u) = approx prec a vs"
      and u_eq: "Some (l', u') = approx prec b vs"
      and inequality: "real_of_float (float_plus_up prec u (-l')) ≤ 0"
    by (cases "approx prec a vs", auto, cases "approx prec b vs", auto)
  from order_trans[OF float_plus_up inequality]
    approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
  show ?case by auto
next
  case (AtLeastAtMost x a b)
  then obtain lx ux l u l' u'
    where x_eq: "Some (lx, ux) = approx prec x vs"
    and l_eq: "Some (l, u) = approx prec a vs"
    and u_eq: "Some (l', u') = approx prec b vs"
    and inequality: "real_of_float (float_plus_up prec u (-lx)) ≤ 0" "real_of_float (float_plus_up prec ux (-l')) ≤ 0"
    by (cases "approx prec x vs", auto,
      cases "approx prec a vs", auto,
      cases "approx prec b vs", auto)
  from order_trans[OF float_plus_up inequality(1)] order_trans[OF float_plus_up inequality(2)]
    approx[OF AtLeastAtMost.prems(2) l_eq] approx[OF AtLeastAtMost.prems(2) u_eq] approx[OF AtLeastAtMost.prems(2) x_eq]
  show ?case by auto
qed auto

lemma approx_form:
  assumes "n = length xs"
    and "approx_form prec f (replicate n None) ss"
  shows "interpret_form f xs"
  using approx_form_aux[OF _ bounded_by_None] assms by auto


subsection ‹Implementing Taylor series expansion›

fun isDERIV :: "nat ⇒ floatarith ⇒ real list ⇒ bool" where
"isDERIV x (Add a b) vs         = (isDERIV x a vs ∧ isDERIV x b vs)" |
"isDERIV x (Mult a b) vs        = (isDERIV x a vs ∧ isDERIV x b vs)" |
"isDERIV x (Minus a) vs         = isDERIV x a vs" |
"isDERIV x (Inverse a) vs       = (isDERIV x a vs ∧ interpret_floatarith a vs ≠ 0)" |
"isDERIV x (Cos a) vs           = isDERIV x a vs" |
"isDERIV x (Arctan a) vs        = isDERIV x a vs" |
"isDERIV x (Min a b) vs         = False" |
"isDERIV x (Max a b) vs         = False" |
"isDERIV x (Abs a) vs           = False" |
"isDERIV x Pi vs                = True" |
"isDERIV x (Sqrt a) vs          = (isDERIV x a vs ∧ interpret_floatarith a vs > 0)" |
"isDERIV x (Exp a) vs           = isDERIV x a vs" |
"isDERIV x (Powr a b) vs        =
    (isDERIV x a vs ∧ isDERIV x b vs ∧ interpret_floatarith a vs > 0)" |
"isDERIV x (Ln a) vs            = (isDERIV x a vs ∧ interpret_floatarith a vs > 0)" |
"isDERIV x (Floor a) vs         = (isDERIV x a vs ∧ interpret_floatarith a vs ∉ ℤ)" |
"isDERIV x (Power a 0) vs       = True" |
"isDERIV x (Power a (Suc n)) vs = isDERIV x a vs" |
"isDERIV x (Num f) vs           = True" |
"isDERIV x (Var n) vs          = True"

fun DERIV_floatarith :: "nat ⇒ floatarith ⇒ floatarith" where
"DERIV_floatarith x (Add a b)         = Add (DERIV_floatarith x a) (DERIV_floatarith x b)" |
"DERIV_floatarith x (Mult a b)        = Add (Mult a (DERIV_floatarith x b)) (Mult (DERIV_floatarith x a) b)" |
"DERIV_floatarith x (Minus a)         = Minus (DERIV_floatarith x a)" |
"DERIV_floatarith x (Inverse a)       = Minus (Mult (DERIV_floatarith x a) (Inverse (Power a 2)))" |
"DERIV_floatarith x (Cos a)           = Minus (Mult (Cos (Add (Mult Pi (Num (Float 1 (- 1)))) (Minus a))) (DERIV_floatarith x a))" |
"DERIV_floatarith x (Arctan a)        = Mult (Inverse (Add (Num 1) (Power a 2))) (DERIV_floatarith x a)" |
"DERIV_floatarith x (Min a b)         = Num 0" |
"DERIV_floatarith x (Max a b)         = Num 0" |
"DERIV_floatarith x (Abs a)           = Num 0" |
"DERIV_floatarith x Pi                = Num 0" |
"DERIV_floatarith x (Sqrt a)          = (Mult (Inverse (Mult (Sqrt a) (Num 2))) (DERIV_floatarith x a))" |
"DERIV_floatarith x (Exp a)           = Mult (Exp a) (DERIV_floatarith x a)" |
"DERIV_floatarith x (Powr a b)        =
   Mult (Powr a b) (Add (Mult (DERIV_floatarith x b) (Ln a)) (Mult (Mult (DERIV_floatarith x a) b) (Inverse a)))" |
"DERIV_floatarith x (Ln a)            = Mult (Inverse a) (DERIV_floatarith x a)" |
"DERIV_floatarith x (Power a 0)       = Num 0" |
"DERIV_floatarith x (Power a (Suc n)) = Mult (Num (Float (int (Suc n)) 0)) (Mult (Power a n) (DERIV_floatarith x a))" |
"DERIV_floatarith x (Floor a)         = Num 0" |
"DERIV_floatarith x (Num f)           = Num 0" |
"DERIV_floatarith x (Var n)          = (if x = n then Num 1 else Num 0)"

lemma has_real_derivative_powr':
  fixes f g :: "real ⇒ real"
  assumes "(f has_real_derivative f') (at x)"
  assumes "(g has_real_derivative g') (at x)"
  assumes "f x > 0"
  defines "h ≡ λx. f x powr g x * (g' * ln (f x) + f' * g x / f x)"
  shows   "((λx. f x powr g x) has_real_derivative h x) (at x)"
proof (subst DERIV_cong_ev[OF refl _ refl])
  from assms have "isCont f x"
    by (simp add: DERIV_continuous)
  hence "f ─x→ f x" by (simp add: continuous_at)
  with ‹f x > 0› have "eventually (λx. f x > 0) (nhds x)"
    by (auto simp: tendsto_at_iff_tendsto_nhds dest: order_tendstoD)
  thus "eventually (λx. f x powr g x = exp (g x * ln (f x))) (nhds x)"
    by eventually_elim (simp add: powr_def)
next
  from assms show "((λx. exp (g x * ln (f x))) has_real_derivative h x) (at x)"
    by (auto intro!: derivative_eq_intros simp: h_def powr_def)
qed

lemma DERIV_floatarith:
  assumes "n < length vs"
  assumes isDERIV: "isDERIV n f (vs[n := x])"
  shows "DERIV (λ x'. interpret_floatarith f (vs[n := x'])) x :>
               interpret_floatarith (DERIV_floatarith n f) (vs[n := x])"
   (is "DERIV (?i f) x :> _")
using isDERIV
proof (induct f arbitrary: x)
  case (Inverse a)
  thus ?case
    by (auto intro!: derivative_eq_intros simp add: algebra_simps power2_eq_square)
next
  case (Cos a)
  thus ?case
    by (auto intro!: derivative_eq_intros
           simp del: interpret_floatarith.simps(5)
           simp add: interpret_floatarith_sin interpret_floatarith.simps(5)[of a])
next
  case (Power a n)
  thus ?case
    by (cases n) (auto intro!: derivative_eq_intros simp del: power_Suc)
next
  case (Floor a)
  thus ?case
    by (auto intro!: derivative_eq_intros DERIV_isCont floor_has_real_derivative)
next
  case (Ln a)
  thus ?case by (auto intro!: derivative_eq_intros simp add: divide_inverse)
next
  case (Var i)
  thus ?case using ‹n < length vs› by auto
next
  case (Powr a b)
  note [derivative_intros] = has_real_derivative_powr'
  from Powr show ?case
    by (auto intro!: derivative_eq_intros simp: field_simps)
qed (auto intro!: derivative_eq_intros)

declare approx.simps[simp del]

fun isDERIV_approx :: "nat ⇒ nat ⇒ floatarith ⇒ (float * float) option list ⇒ bool" where
"isDERIV_approx prec x (Add a b) vs         = (isDERIV_approx prec x a vs ∧ isDERIV_approx prec x b vs)" |
"isDERIV_approx prec x (Mult a b) vs        = (isDERIV_approx prec x a vs ∧ isDERIV_approx prec x b vs)" |
"isDERIV_approx prec x (Minus a) vs         = isDERIV_approx prec x a vs" |
"isDERIV_approx prec x (Inverse a) vs       =
  (isDERIV_approx prec x a vs ∧ (case approx prec a vs of Some (l, u) ⇒ 0 < l ∨ u < 0 | None ⇒ False))" |
"isDERIV_approx prec x (Cos a) vs           = isDERIV_approx prec x a vs" |
"isDERIV_approx prec x (Arctan a) vs        = isDERIV_approx prec x a vs" |
"isDERIV_approx prec x (Min a b) vs         = False" |
"isDERIV_approx prec x (Max a b) vs         = False" |
"isDERIV_approx prec x (Abs a) vs           = False" |
"isDERIV_approx prec x Pi vs                = True" |
"isDERIV_approx prec x (Sqrt a) vs          =
  (isDERIV_approx prec x a vs ∧ (case approx prec a vs of Some (l, u) ⇒ 0 < l | None ⇒ False))" |
"isDERIV_approx prec x (Exp a) vs           = isDERIV_approx prec x a vs" |
"isDERIV_approx prec x (Powr a b) vs        =
  (isDERIV_approx prec x a vs ∧ isDERIV_approx prec x b vs ∧ (case approx prec a vs of Some (l, u) ⇒ 0 < l | None ⇒ False))" |
"isDERIV_approx prec x (Ln a) vs            =
  (isDERIV_approx prec x a vs ∧ (case approx prec a vs of Some (l, u) ⇒ 0 < l | None ⇒ False))" |
"isDERIV_approx prec x (Power a 0) vs       = True" |
"isDERIV_approx prec x (Floor a) vs         =
  (isDERIV_approx prec x a vs ∧ (case approx prec a vs of Some (l, u) ⇒ l > floor u ∧ u < ceiling l | None ⇒ False))" |
"isDERIV_approx prec x (Power a (Suc n)) vs = isDERIV_approx prec x a vs" |
"isDERIV_approx prec x (Num f) vs           = True" |
"isDERIV_approx prec x (Var n) vs           = True"

lemma isDERIV_approx:
  assumes "bounded_by xs vs"
    and isDERIV_approx: "isDERIV_approx prec x f vs"
  shows "isDERIV x f xs"
  using isDERIV_approx
proof (induct f)
  case (Inverse a)
  then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
    and *: "0 < l ∨ u < 0"
    by (cases "approx prec a vs") auto
  with approx[OF ‹bounded_by xs vs› approx_Some]
  have "interpret_floatarith a xs ≠ 0" by auto
  thus ?case using Inverse by auto
next
  case (Ln a)
  then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
    and *: "0 < l"
    by (cases "approx prec a vs") auto
  with approx[OF ‹bounded_by xs vs› approx_Some]
  have "0 < interpret_floatarith a xs" by auto
  thus ?case using Ln by auto
next
  case (Sqrt a)
  then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
    and *: "0 < l"
    by (cases "approx prec a vs") auto
  with approx[OF ‹bounded_by xs vs› approx_Some]
  have "0 < interpret_floatarith a xs" by auto
  thus ?case using Sqrt by auto
next
  case (Power a n)
  thus ?case by (cases n) auto
next
  case (Powr a b)
  from Powr obtain l1 u1 where a: "Some (l1, u1) = approx prec a vs" and pos: "0 < l1"
    by (cases "approx prec a vs") auto
  with approx[OF ‹bounded_by xs vs› a]
    have "0 < interpret_floatarith a xs" by auto
  with Powr show ?case by auto
next
  case (Floor a)
  then obtain l u where approx_Some: "Some (l, u) = approx prec a vs"
    and "real_of_int ⌊real_of_float u⌋ < real_of_float l" "real_of_float u < real_of_int ⌈real_of_float l⌉"
    and "isDERIV x a xs"
    by (cases "approx prec a vs") auto
  with approx[OF ‹bounded_by xs vs› approx_Some] le_floor_iff
  show ?case
    by (force elim!: Ints_cases)
qed auto

lemma bounded_by_update_var:
  assumes "bounded_by xs vs"
    and "vs ! i = Some (l, u)"
    and bnd: "x ∈ { real_of_float l .. real_of_float u }"
  shows "bounded_by (xs[i := x]) vs"
proof (cases "i < length xs")
  case False
  thus ?thesis
    using ‹bounded_by xs vs› by auto
next
  case True
  let ?xs = "xs[i := x]"
  from True have "i < length ?xs" by auto
  have "case vs ! j of None ⇒ True | Some (l, u) ⇒ ?xs ! j ∈ {real_of_float l .. real_of_float u}"
    if "j < length vs" for j
  proof (cases "vs ! j")
    case None
    then show ?thesis by simp
  next
    case (Some b)
    thus ?thesis
    proof (cases "i = j")
      case True
      thus ?thesis using ‹vs ! i = Some (l, u)› Some and bnd ‹i < length ?xs›
        by auto
    next
      case False
      thus ?thesis
        using ‹bounded_by xs vs›[THEN bounded_byE, OF ‹j < length vs›] Some by auto
    qed
  qed
  thus ?thesis
    unfolding bounded_by_def by auto
qed

lemma isDERIV_approx':
  assumes "bounded_by xs vs"
    and vs_x: "vs ! x = Some (l, u)"
    and X_in: "X ∈ {real_of_float l .. real_of_float u}"
    and approx: "isDERIV_approx prec x f vs"
  shows "isDERIV x f (xs[x := X])"
proof -
  from bounded_by_update_var[OF ‹bounded_by xs vs› vs_x X_in] approx
  show ?thesis by (rule isDERIV_approx)
qed

lemma DERIV_approx:
  assumes "n < length xs"
    and bnd: "bounded_by xs vs"
    and isD: "isDERIV_approx prec n f vs"
    and app: "Some (l, u) = approx prec (DERIV_floatarith n f) vs" (is "_ = approx _ ?D _")
  shows "∃(x::real). l ≤ x ∧ x ≤ u ∧
             DERIV (λ x. interpret_floatarith f (xs[n := x])) (xs!n) :> x"
         (is "∃ x. _ ∧ _ ∧ DERIV (?i f) _ :> _")
proof (rule exI[of _ "?i ?D (xs!n)"], rule conjI[OF _ conjI])
  let "?i f" = "λx. interpret_floatarith f (xs[n := x])"
  from approx[OF bnd app]
  show "l ≤ ?i ?D (xs!n)" and "?i ?D (xs!n) ≤ u"
    using ‹n < length xs› by auto
  from DERIV_floatarith[OF ‹n < length xs›, of f "xs!n"] isDERIV_approx[OF bnd isD]
  show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))"
    by simp
qed

lemma lift_bin_aux:
  assumes lift_bin_Some: "Some (l, u) = lift_bin a b f"
  obtains l1 u1 l2 u2
  where "a = Some (l1, u1)"
    and "b = Some (l2, u2)"
    and "f l1 u1 l2 u2 = Some (l, u)"
  using assms by (cases a, simp, cases b, simp, auto)


fun approx_tse where
"approx_tse prec n 0 c k f bs = approx prec f bs" |
"approx_tse prec n (Suc s) c k f bs =
  (if isDERIV_approx prec n f bs then
    lift_bin (approx prec f (bs[n := Some (c,c)]))
             (approx_tse prec n s c (Suc k) (DERIV_floatarith n f) bs)
             (λ l1 u1 l2 u2. approx prec
                 (Add (Var 0)
                      (Mult (Inverse (Num (Float (int k) 0)))
                                 (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
                                       (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), bs!n])
  else approx prec f bs)"

lemma bounded_by_Cons:
  assumes bnd: "bounded_by xs vs"
    and x: "x ∈ { real_of_float l .. real_of_float u }"
  shows "bounded_by (x#xs) ((Some (l, u))#vs)"
proof -
  have "case ((Some (l,u))#vs) ! i of Some (l, u) ⇒ (x#xs)!i ∈ { real_of_float l .. real_of_float u } | None ⇒ True"
    if *: "i < length ((Some (l, u))#vs)" for i
  proof (cases i)
    case 0
    with x show ?thesis by auto
  next
    case (Suc i)
    with * have "i < length vs" by auto
    from bnd[THEN bounded_byE, OF this]
    show ?thesis unfolding Suc nth_Cons_Suc .
  qed
  thus ?thesis
    by (auto simp add: bounded_by_def)
qed

lemma approx_tse_generic:
  assumes "bounded_by xs vs"
    and bnd_c: "bounded_by (xs[x := c]) vs"
    and "x < length vs" and "x < length xs"
    and bnd_x: "vs ! x = Some (lx, ux)"
    and ate: "Some (l, u) = approx_tse prec x s c k f vs"
  shows "∃ n. (∀ m < n. ∀ (z::real) ∈ {lx .. ux}.
      DERIV (λ y. interpret_floatarith ((DERIV_floatarith x ^^ m) f) (xs[x := y])) z :>
            (interpret_floatarith ((DERIV_floatarith x ^^ (Suc m)) f) (xs[x := z])))
   ∧ (∀ (t::real) ∈ {lx .. ux}.  (∑ i = 0..<n. inverse (real (∏ j ∈ {k..<k+i}. j)) *
                  interpret_floatarith ((DERIV_floatarith x ^^ i) f) (xs[x := c]) *
                  (xs!x - c)^i) +
      inverse (real (∏ j ∈ {k..<k+n}. j)) *
      interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := t]) *
      (xs!x - c)^n ∈ {l .. u})" (is "∃ n. ?taylor f k l u n")
  using ate
proof (induct s arbitrary: k f l u)
  case 0
  {
    fix t::real assume "t ∈ {lx .. ux}"
    note bounded_by_update_var[OF ‹bounded_by xs vs› bnd_x this]
    from approx[OF this 0[unfolded approx_tse.simps]]
    have "(interpret_floatarith f (xs[x := t])) ∈ {l .. u}"
      by (auto simp add: algebra_simps)
  }
  thus ?case by (auto intro!: exI[of _ 0])
next
  case (Suc s)
  show ?case
  proof (cases "isDERIV_approx prec x f vs")
    case False
    note ap = Suc.prems[unfolded approx_tse.simps if_not_P[OF False]]
    {
      fix t::real assume "t ∈ {lx .. ux}"
      note bounded_by_update_var[OF ‹bounded_by xs vs› bnd_x this]
      from approx[OF this ap]
      have "(interpret_floatarith f (xs[x := t])) ∈ {l .. u}"
        by (auto simp add: algebra_simps)
    }
    thus ?thesis by (auto intro!: exI[of _ 0])
  next
    case True
    with Suc.prems
    obtain l1 u1 l2 u2
      where a: "Some (l1, u1) = approx prec f (vs[x := Some (c,c)])"
        and ate: "Some (l2, u2) = approx_tse prec x s c (Suc k) (DERIV_floatarith x f) vs"
        and final: "Some (l, u) = approx prec
          (Add (Var 0)
               (Mult (Inverse (Num (Float (int k) 0)))
                     (Mult (Add (Var (Suc (Suc 0))) (Minus (Num c)))
                           (Var (Suc 0))))) [Some (l1, u1), Some (l2, u2), vs!x]"
      by (auto elim!: lift_bin_aux)

    from bnd_c ‹x < length xs›
    have bnd: "bounded_by (xs[x:=c]) (vs[x:= Some (c,c)])"
      by (auto intro!: bounded_by_update)

    from approx[OF this a]
    have f_c: "interpret_floatarith ((DERIV_floatarith x ^^ 0) f) (xs[x := c]) ∈ { l1 .. u1 }"
              (is "?f 0 (real_of_float c) ∈ _")
      by auto

    have funpow_Suc[symmetric]: "(f ^^ Suc n) x = (f ^^ n) (f x)"
      for f :: "'a ⇒ 'a" and n :: nat and x :: 'a
      by (induct n) auto
    from Suc.hyps[OF ate, unfolded this] obtain n
      where DERIV_hyp: "⋀m z. ⟦ m < n ; (z::real) ∈ { lx .. ux } ⟧ ⟹
        DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z"
      and hyp: "∀t ∈ {real_of_float lx .. real_of_float ux}.
        (∑ i = 0..<n. inverse (real (∏ j ∈ {Suc k..<Suc k + i}. j)) * ?f (Suc i) c * (xs!x - c)^i) +
          inverse (real (∏ j ∈ {Suc k..<Suc k + n}. j)) * ?f (Suc n) t * (xs!x - c)^n ∈ {l2 .. u2}"
          (is "∀ t ∈ _. ?X (Suc k) f n t ∈ _")
      by blast

    have DERIV: "DERIV (?f m) z :> ?f (Suc m) z"
      if "m < Suc n" and bnd_z: "z ∈ { lx .. ux }" for m and z::real
    proof (cases m)
      case 0
      with DERIV_floatarith[OF ‹x < length xs›
        isDERIV_approx'[OF ‹bounded_by xs vs› bnd_x bnd_z True]]
      show ?thesis by simp
    next
      case (Suc m')
      hence "m' < n"
        using ‹m < Suc n› by auto
      from DERIV_hyp[OF this bnd_z] show ?thesis
        using Suc by simp
    qed

    have "⋀k i. k < i ⟹ {k ..< i} = insert k {Suc k ..< i}" by auto
    hence prod_head_Suc: "⋀k i. ∏{k ..< k + Suc i} = k * ∏{Suc k ..< Suc k + i}"
      by auto
    have sum_move0: "⋀k F. sum F {0..<Suc k} = F 0 + sum (λ k. F (Suc k)) {0..<k}"
      unfolding sum_shift_bounds_Suc_ivl[symmetric]
      unfolding sum_head_upt_Suc[OF zero_less_Suc] ..
    define C where "C = xs!x - c"

    {
      fix t::real assume t: "t ∈ {lx .. ux}"
      hence "bounded_by [xs!x] [vs!x]"
        using ‹bounded_by xs vs›[THEN bounded_byE, OF ‹x < length vs›]
        by (cases "vs!x", auto simp add: bounded_by_def)

      with hyp[THEN bspec, OF t] f_c
      have "bounded_by [?f 0 c, ?X (Suc k) f n t, xs!x] [Some (l1, u1), Some (l2, u2), vs!x]"
        by (auto intro!: bounded_by_Cons)
      from approx[OF this final, unfolded atLeastAtMost_iff[symmetric]]
      have "?X (Suc k) f n t * (xs!x - real_of_float c) * inverse k + ?f 0 c ∈ {l .. u}"
        by (auto simp add: algebra_simps)
      also have "?X (Suc k) f n t * (xs!x - real_of_float c) * inverse (real k) + ?f 0 c =
               (∑ i = 0..<Suc n. inverse (real (∏ j ∈ {k..<k+i}. j)) * ?f i c * (xs!x - c)^i) +
               inverse (real (∏ j ∈ {k..<k+Suc n}. j)) * ?f (Suc n) t * (xs!x - c)^Suc n" (is "_ = ?T")
        unfolding funpow_Suc C_def[symmetric] sum_move0 prod_head_Suc
        by (auto simp add: algebra_simps)
          (simp only: mult.left_commute [of _ "inverse (real k)"] sum_distrib_left [symmetric])
      finally have "?T ∈ {l .. u}" .
    }
    thus ?thesis using DERIV by blast
  qed
qed

lemma prod_fact: "real (∏ {1..<1 + k}) = fact (k :: nat)"
  by (simp add: fact_prod atLeastLessThanSuc_atLeastAtMost)

lemma approx_tse:
  assumes "bounded_by xs vs"
    and bnd_x: "vs ! x = Some (lx, ux)"
    and bnd_c: "real_of_float c ∈ {lx .. ux}"
    and "x < length vs" and "x < length xs"
    and ate: "Some (l, u) = approx_tse prec x s c 1 f vs"
  shows "interpret_floatarith f xs ∈ {l .. u}"
proof -
  define F where [abs_def]: "F n z =
    interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])" for n z
  hence F0: "F 0 = (λ z. interpret_floatarith f (xs[x := z]))" by auto

  hence "bounded_by (xs[x := c]) vs" and "x < length vs" "x < length xs"
    using ‹bounded_by xs vs› bnd_x bnd_c ‹x < length vs› ‹x < length xs›
    by (auto intro!: bounded_by_update_var)

  from approx_tse_generic[OF ‹bounded_by xs vs› this bnd_x ate]
  obtain n
    where DERIV: "∀ m z. m < n ∧ real_of_float lx ≤ z ∧ z ≤ real_of_float ux ⟶ DERIV (F m) z :> F (Suc m) z"
    and hyp: "⋀ (t::real). t ∈ {lx .. ux} ⟹
           (∑ j = 0..<n. inverse(fact j) * F j c * (xs!x - c)^j) +
             inverse ((fact n)) * F n t * (xs!x - c)^n
             ∈ {l .. u}" (is "⋀ t. _ ⟹ ?taylor t ∈ _")
    unfolding F_def atLeastAtMost_iff[symmetric] prod_fact
    by blast

  have bnd_xs: "xs ! x ∈ { lx .. ux }"
    using ‹bounded_by xs vs›[THEN bounded_byE, OF ‹x < length vs›] bnd_x by auto

  show ?thesis
  proof (cases n)
    case 0
    thus ?thesis
      using hyp[OF bnd_xs] unfolding F_def by auto
  next
    case (Suc n')
    show ?thesis
    proof (cases "xs ! x = c")
      case True
      from True[symmetric] hyp[OF bnd_xs] Suc show ?thesis
        unfolding F_def Suc sum_head_upt_Suc[OF zero_less_Suc] sum_shift_bounds_Suc_ivl
        by auto
    next
      case False
      have "lx ≤ real_of_float c" "real_of_float c ≤ ux" "lx ≤ xs!x" "xs!x ≤ ux"
        using Suc bnd_c ‹bounded_by xs vs›[THEN bounded_byE, OF ‹x < length vs›] bnd_x by auto
      from taylor[OF zero_less_Suc, of F, OF F0 DERIV[unfolded Suc] this False]
      obtain t::real where t_bnd: "if xs ! x < c then xs ! x < t ∧ t < c else c < t ∧ t < xs ! x"
        and fl_eq: "interpret_floatarith f (xs[x := xs ! x]) =
           (∑m = 0..<Suc n'. F m c / (fact m) * (xs ! x - c) ^ m) +
           F (Suc n') t / (fact (Suc n')) * (xs ! x - c) ^ Suc n'"
        unfolding atLeast0LessThan by blast

      from t_bnd bnd_xs bnd_c have *: "t ∈ {lx .. ux}"
        by (cases "xs ! x < c") auto

      have "interpret_floatarith f (xs[x := xs ! x]) = ?taylor t"
        unfolding fl_eq Suc by (auto simp add: algebra_simps divide_inverse)
      also have "… ∈ {l .. u}"
        using * by (rule hyp)
      finally show ?thesis
        by simp
    qed
  qed
qed

fun approx_tse_form' where
"approx_tse_form' prec t f 0 l u cmp =
  (case approx_tse prec 0 t ((l + u) * Float 1 (- 1)) 1 f [Some (l, u)]
     of Some (l, u) ⇒ cmp l u | None ⇒ False)" |
"approx_tse_form' prec t f (Suc s) l u cmp =
  (let m = (l + u) * Float 1 (- 1)
   in (if approx_tse_form' prec t f s l m cmp then
      approx_tse_form' prec t f s m u cmp else False))"

lemma approx_tse_form':
  fixes x :: real
  assumes "approx_tse_form' prec t f s l u cmp"
    and "x ∈ {l .. u}"
  shows "∃l' u' ly uy. x ∈ {l' .. u'} ∧ real_of_float l ≤ l' ∧ u' ≤ real_of_float u ∧ cmp ly uy ∧
    approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)"
  using assms
proof (induct s arbitrary: l u)
  case 0
  then obtain ly uy
    where *: "approx_tse prec 0 t ((l + u) * Float 1 (- 1)) 1 f [Some (l, u)] = Some (ly, uy)"
    and **: "cmp ly uy" by (auto elim!: case_optionE)
  with 0 show ?case by auto
next
  case (Suc s)
  let ?m = "(l + u) * Float 1 (- 1)"
  from Suc.prems
  have l: "approx_tse_form' prec t f s l ?m cmp"
    and u: "approx_tse_form' prec t f s ?m u cmp"
    by (auto simp add: Let_def lazy_conj)

  have m_l: "real_of_float l ≤ ?m" and m_u: "?m ≤ real_of_float u"
    unfolding less_eq_float_def using Suc.prems by auto
  with ‹x ∈ { l .. u }› consider "x ∈ { l .. ?m}" | "x ∈ {?m .. u}"
    by atomize_elim auto
  thus ?case
  proof cases
    case 1
    from Suc.hyps[OF l this]
    obtain l' u' ly uy where
      "x ∈ {l' .. u'} ∧ real_of_float l ≤ l' ∧ real_of_float u' ≤ ?m ∧ cmp ly uy ∧
        approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)"
      by blast
    with m_u show ?thesis
      by (auto intro!: exI)
  next
    case 2
    from Suc.hyps[OF u this]
    obtain l' u' ly uy where
      "x ∈ { l' .. u' } ∧ ?m ≤ real_of_float l' ∧ u' ≤ real_of_float u ∧ cmp ly uy ∧
        approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)"
      by blast
    with m_u show ?thesis
      by (auto intro!: exI)
  qed
qed

lemma approx_tse_form'_less:
  fixes x :: real
  assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (λ l u. 0 < l)"
    and x: "x ∈ {l .. u}"
  shows "interpret_floatarith b [x] < interpret_floatarith a [x]"
proof -
  from approx_tse_form'[OF tse x]
  obtain l' u' ly uy
    where x': "x ∈ {l' .. u'}"
    and "real_of_float l ≤ real_of_float l'"
    and "real_of_float u' ≤ real_of_float u" and "0 < ly"
    and tse: "approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
    by blast

  hence "bounded_by [x] [Some (l', u')]"
    by (auto simp add: bounded_by_def)
  from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
  have "ly ≤ interpret_floatarith a [x] - interpret_floatarith b [x]"
    by auto
  from order_less_le_trans[OF _ this, of 0] ‹0 < ly› show ?thesis
    by auto
qed

lemma approx_tse_form'_le:
  fixes x :: real
  assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (λ l u. 0 ≤ l)"
    and x: "x ∈ {l .. u}"
  shows "interpret_floatarith b [x] ≤ interpret_floatarith a [x]"
proof -
  from approx_tse_form'[OF tse x]
  obtain l' u' ly uy
    where x': "x ∈ {l' .. u'}"
    and "l ≤ real_of_float l'"
    and "real_of_float u' ≤ u" and "0 ≤ ly"
    and tse: "approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)"
    by blast

  hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def)
  from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x'
  have "ly ≤ interpret_floatarith a [x] - interpret_floatarith b [x]"
    by auto
  from order_trans[OF _ this, of 0] ‹0 ≤ ly› show ?thesis
    by auto
qed

fun approx_tse_concl where
"approx_tse_concl prec t (Less lf rt) s l u l' u' ⟷
    approx_tse_form' prec t (Add rt (Minus lf)) s l u' (λ l u. 0 < l)" |
"approx_tse_concl prec t (LessEqual lf rt) s l u l' u' ⟷
    approx_tse_form' prec t (Add rt (Minus lf)) s l u' (λ l u. 0 ≤ l)" |
"approx_tse_concl prec t (AtLeastAtMost x lf rt) s l u l' u' ⟷
    (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (λ l u. 0 ≤ l) then
      approx_tse_form' prec t (Add rt (Minus x)) s l u' (λ l u. 0 ≤ l) else False)" |
"approx_tse_concl prec t (Conj f g) s l u l' u' ⟷
    approx_tse_concl prec t f s l u l' u' ∧ approx_tse_concl prec t g s l u l' u'" |
"approx_tse_concl prec t (Disj f g) s l u l' u' ⟷
    approx_tse_concl prec t f s l u l' u' ∨ approx_tse_concl prec t g s l u l' u'" |
"approx_tse_concl _ _ _ _ _ _ _ _ ⟷ False"

definition
  "approx_tse_form prec t s f =
    (case f of
      Bound x a b f ⇒
        x = Var 0 ∧
        (case (approx prec a [None], approx prec b [None]) of
          (Some (l, u), Some (l', u')) ⇒ approx_tse_concl prec t f s l u l' u'
        | _ ⇒ False)
    | _ ⇒ False)"

lemma approx_tse_form:
  assumes "approx_tse_form prec t s f"
  shows "interpret_form f [x]"
proof (cases f)
  case f_def: (Bound i a b f')
  with assms obtain l u l' u'
    where a: "approx prec a [None] = Some (l, u)"
    and b: "approx prec b [None] = Some (l', u')"
    unfolding approx_tse_form_def by (auto elim!: case_optionE)

  from f_def assms have "i = Var 0"
    unfolding approx_tse_form_def by auto
  hence i: "interpret_floatarith i [x] = x" by auto

  {
    let ?f = "λz. interpret_floatarith z [x]"
    assume "?f i ∈ { ?f a .. ?f b }"
    with approx[OF _ a[symmetric], of "[x]"] approx[OF _ b[symmetric], of "[x]"]
    have bnd: "x ∈ { l .. u'}" unfolding bounded_by_def i by auto

    have "interpret_form f' [x]"
      using assms[unfolded f_def]
    proof (induct f')
      case (Less lf rt)
      with a b
      have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (λ l u. 0 < l)"
        unfolding approx_tse_form_def by auto
      from approx_tse_form'_less[OF this bnd]
      show ?case using Less by auto
    next
      case (LessEqual lf rt)
      with f_def a b assms
      have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (λ l u. 0 ≤ l)"
        unfolding approx_tse_form_def by auto
      from approx_tse_form'_le[OF this bnd]
      show ?case using LessEqual by auto
    next
      case (AtLeastAtMost x lf rt)
      with f_def a b assms
      have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (λ l u. 0 ≤ l)"
        and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (λ l u. 0 ≤ l)"
        unfolding approx_tse_form_def lazy_conj by (auto split: if_split_asm)
      from approx_tse_form'_le[OF this(1) bnd] approx_tse_form'_le[OF this(2) bnd]
      show ?case using AtLeastAtMost by auto
    qed (auto simp: f_def approx_tse_form_def elim!: case_optionE)
  }
  thus ?thesis unfolding f_def by auto
qed (insert assms, auto simp add: approx_tse_form_def)

text ‹@{term approx_form_eval} is only used for the {\tt value}-command.›

fun approx_form_eval :: "nat ⇒ form ⇒ (float * float) option list ⇒ (float * float) option list" where
"approx_form_eval prec (Bound (Var n) a b f) bs =
   (case (approx prec a bs, approx prec b bs)
   of (Some (l, _), Some (_, u)) ⇒ approx_form_eval prec f (bs[n := Some (l, u)])
    | _ ⇒ bs)" |
"approx_form_eval prec (Assign (Var n) a f) bs =
   (case (approx prec a bs)
   of (Some (l, u)) ⇒ approx_form_eval prec f (bs[n := Some (l, u)])
    | _ ⇒ bs)" |
"approx_form_eval prec (Less a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
"approx_form_eval prec (LessEqual a b) bs = bs @ [approx prec a bs, approx prec b bs]" |
"approx_form_eval prec (AtLeastAtMost x a b) bs =
   bs @ [approx prec x bs, approx prec a bs, approx prec b bs]" |
"approx_form_eval _ _ bs = bs"


subsection ‹Implement proof method \texttt{approximation}›

oracle approximation_oracle = ‹fn (thy, t) =>
let
  fun bad t = error ("Bad term: " ^ Syntax.string_of_term_global thy t);

  fun term_of_bool true = @{term True}
    | term_of_bool false = @{term False};

  val mk_int = HOLogic.mk_number @{typ int} o @{code integer_of_int};
  fun dest_int (@{term int_of_integer} $ j) = @{code int_of_integer} (snd (HOLogic.dest_number j))
    | dest_int i = @{code int_of_integer} (snd (HOLogic.dest_number i));

  fun term_of_float (@{code Float} (k, l)) =
    @{term Float} $ mk_int k $ mk_int l;

  fun term_of_float_float_option NONE = @{term "None :: (float × float) option"}
    | term_of_float_float_option (SOME ff) = @{term "Some :: float × float ⇒ _"}
        $ HOLogic.mk_prod (apply2 term_of_float ff);

  val term_of_float_float_option_list =
    HOLogic.mk_list @{typ "(float × float) option"} o map term_of_float_float_option;

  fun nat_of_term t = @{code nat_of_integer}
    (HOLogic.dest_nat t handle TERM _ => snd (HOLogic.dest_number t));

  fun float_of_term (@{term Float} $ k $ l) =
        @{code Float} (dest_int k, dest_int l)
    | float_of_term t = bad t;

  fun floatarith_of_term (@{term Add} $ a $ b) = @{code Add} (floatarith_of_term a, floatarith_of_term b)
    | floatarith_of_term (@{term Minus} $ a) = @{code Minus} (floatarith_of_term a)
    | floatarith_of_term (@{term Mult} $ a $ b) = @{code Mult} (floatarith_of_term a, floatarith_of_term b)
    | floatarith_of_term (@{term Inverse} $ a) = @{code Inverse} (floatarith_of_term a)
    | floatarith_of_term (@{term Cos} $ a) = @{code Cos} (floatarith_of_term a)
    | floatarith_of_term (@{term Arctan} $ a) = @{code Arctan} (floatarith_of_term a)
    | floatarith_of_term (@{term Abs} $ a) = @{code Abs} (floatarith_of_term a)
    | floatarith_of_term (@{term Max} $ a $ b) = @{code Max} (floatarith_of_term a, floatarith_of_term b)
    | floatarith_of_term (@{term Min} $ a $ b) = @{code Min} (floatarith_of_term a, floatarith_of_term b)
    | floatarith_of_term @{term Pi} = @{code Pi}
    | floatarith_of_term (@{term Sqrt} $ a) = @{code Sqrt} (floatarith_of_term a)
    | floatarith_of_term (@{term Exp} $ a) = @{code Exp} (floatarith_of_term a)
    | floatarith_of_term (@{term Powr} $ a $ b) = @{code Powr} (floatarith_of_term a, floatarith_of_term b)
    | floatarith_of_term (@{term Ln} $ a) = @{code Ln} (floatarith_of_term a)
    | floatarith_of_term (@{term Power} $ a $ n) =
        @{code Power} (floatarith_of_term a, nat_of_term n)
    | floatarith_of_term (@{term Floor} $ a) = @{code Floor} (floatarith_of_term a)
    | floatarith_of_term (@{term Var} $ n) = @{code Var} (nat_of_term n)
    | floatarith_of_term (@{term Num} $ m) = @{code Num} (float_of_term m)
    | floatarith_of_term t = bad t;

  fun form_of_term (@{term Bound} $ a $ b $ c $ p) = @{code Bound}
        (floatarith_of_term a, floatarith_of_term b, floatarith_of_term c, form_of_term p)
    | form_of_term (@{term Assign} $ a $ b $ p) = @{code Assign}
        (floatarith_of_term a, floatarith_of_term b, form_of_term p)
    | form_of_term (@{term Less} $ a $ b) = @{code Less}
        (floatarith_of_term a, floatarith_of_term b)
    | form_of_term (@{term LessEqual} $ a $ b) = @{code LessEqual}
        (floatarith_of_term a, floatarith_of_term b)
    | form_of_term (@{term Conj} $ a $ b) = @{code Conj}
        (form_of_term a, form_of_term b)
    | form_of_term (@{term Disj} $ a $ b) = @{code Disj}
        (form_of_term a, form_of_term b)
    | form_of_term (@{term AtLeastAtMost} $ a $ b $ c) = @{code AtLeastAtMost}
        (floatarith_of_term a, floatarith_of_term b, floatarith_of_term c)
    | form_of_term t = bad t;

  fun float_float_option_of_term @{term "None :: (float × float) option"} = NONE
    | float_float_option_of_term (@{term "Some :: float × float ⇒ _"} $ ff) =
        SOME (apply2 float_of_term (HOLogic.dest_prod ff))
    | float_float_option_of_term (@{term approx'} $ n $ a $ ffs) = @{code approx'}
        (nat_of_term n) (floatarith_of_term a) (float_float_option_list_of_term ffs)
    | float_float_option_of_term t = bad t
  and float_float_option_list_of_term
        (@{term "replicate :: _ ⇒ (float × float) option ⇒ _"} $ n $ @{term "None :: (float × float) option"}) =
          @{code replicate} (nat_of_term n) NONE
    | float_float_option_list_of_term (@{term approx_form_eval} $ n $ p $ ffs) =
        @{code approx_form_eval} (nat_of_term n) (form_of_term p) (float_float_option_list_of_term ffs)
    | float_float_option_list_of_term t = map float_float_option_of_term
        (HOLogic.dest_list t);

  val nat_list_of_term = map nat_of_term o HOLogic.dest_list ;

  fun bool_of_term (@{term approx_form} $ n $ p $ ffs $ ms) = @{code approx_form}
        (nat_of_term n) (form_of_term p) (float_float_option_list_of_term ffs) (nat_list_of_term ms)
    | bool_of_term (@{term approx_tse_form} $ m $ n $ q $ p) =
        @{code approx_tse_form} (nat_of_term m) (nat_of_term n) (nat_of_term q) (form_of_term p)
    | bool_of_term t = bad t;

  fun eval t = case fastype_of t
   of @{typ bool} =>
        (term_of_bool o bool_of_term) t
    | @{typ "(float × float) option"} =>
        (term_of_float_float_option o float_float_option_of_term) t
    | @{typ "(float × float) option list"} =>
        (term_of_float_float_option_list o float_float_option_list_of_term) t
    | _ => bad t;

  val normalize = eval o Envir.beta_norm o Envir.eta_long [];

in Thm.global_cterm_of thy (Logic.mk_equals (t, normalize t)) end
›

lemma intervalE: "a ≤ x ∧ x ≤ b ⟹ ⟦ x ∈ { a .. b } ⟹ P⟧ ⟹ P"
  by auto

lemma meta_eqE: "x ≡ a ⟹ ⟦ x = a ⟹ P⟧ ⟹ P"
  by auto

named_theorems approximation_preproc

lemma approximation_preproc_floatarith[approximation_preproc]:
  "0 = real_of_float 0"
  "1 = real_of_float 1"
  "0 = Float 0 0"
  "1 = Float 1 0"
  "numeral a = Float (numeral a) 0"
  "numeral a = real_of_float (numeral a)"
  "x - y = x + - y"
  "x / y = x * inverse y"
  "ceiling x = - floor (- x)"
  "log x y = ln y * inverse (ln x)"
  "sin x = cos (pi / 2 - x)"
  "tan x = sin x / cos x"
  by (simp_all add: inverse_eq_divide ceiling_def log_def sin_cos_eq tan_def real_of_float_eq)

lemma approximation_preproc_int[approximation_preproc]:
  "real_of_int 0 = real_of_float 0"
  "real_of_int 1 = real_of_float 1"
  "real_of_int (i + j) = real_of_int i + real_of_int j"
  "real_of_int (- i) = - real_of_int i"
  "real_of_int (i - j) = real_of_int i - real_of_int j"
  "real_of_int (i * j) = real_of_int i * real_of_int j"
  "real_of_int (i div j) = real_of_int (floor (real_of_int i / real_of_int j))"
  "real_of_int (min i j) = min (real_of_int i) (real_of_int j)"
  "real_of_int (max i j) = max (real_of_int i) (real_of_int j)"
  "real_of_int (abs i) = abs (real_of_int i)"
  "real_of_int (i ^ n) = (real_of_int i) ^ n"
  "real_of_int (numeral a) = real_of_float (numeral a)"
  "i mod j = i - i div j * j"
  "i = j ⟷ real_of_int i = real_of_int j"
  "i ≤ j ⟷ real_of_int i ≤ real_of_int j"
  "i < j ⟷ real_of_int i < real_of_int j"
  "i ∈ {j .. k} ⟷ real_of_int i ∈ {real_of_int j .. real_of_int k}"
  by (simp_all add: floor_divide_of_int_eq minus_div_mult_eq_mod [symmetric])

lemma approximation_preproc_nat[approximation_preproc]:
  "real 0 = real_of_float 0"
  "real 1 = real_of_float 1"
  "real (i + j) = real i + real j"
  "real (i - j) = max (real i - real j) 0"
  "real (i * j) = real i * real j"
  "real (i div j) = real_of_int (floor (real i / real j))"
  "real (min i j) = min (real i) (real j)"
  "real (max i j) = max (real i) (real j)"
  "real (i ^ n) = (real i) ^ n"
  "real (numeral a) = real_of_float (numeral a)"
  "i mod j = i - i div j * j"
  "n = m ⟷ real n = real m"
  "n ≤ m ⟷ real n ≤ real m"
  "n < m ⟷ real n < real m"
  "n ∈ {m .. l} ⟷ real n ∈ {real m .. real l}"
  by (simp_all add: real_div_nat_eq_floor_of_divide minus_div_mult_eq_mod [symmetric])

ML_file "approximation.ML"

method_setup approximation = ‹
  let
    val free =
      Args.context -- Args.term >> (fn (_, Free (n, _)) => n | (ctxt, t) =>
        error ("Bad free variable: " ^ Syntax.string_of_term ctxt t));
  in
    Scan.lift Parse.nat --
    Scan.optional (Scan.lift (Args.$$$ "splitting" |-- Args.colon)
      |-- Parse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift Parse.nat)) [] --
    Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon) |--
    (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift Parse.nat)) >>
    (fn ((prec, splitting), taylor) => fn ctxt =>
      SIMPLE_METHOD' (Approximation.approximation_tac prec splitting taylor ctxt))
  end
› "real number approximation"


section "Quickcheck Generator"

lemma approximation_preproc_push_neg[approximation_preproc]:
  fixes a b::real
  shows
    "¬ (a < b) ⟷ b ≤ a"
    "¬ (a ≤ b) ⟷ b < a"
    "¬ (a = b) ⟷ b < a ∨ a < b"
    "¬ (p ∧ q) ⟷ ¬ p ∨ ¬ q"
    "¬ (p ∨ q) ⟷ ¬ p ∧ ¬ q"
    "¬ ¬ q ⟷ q"
  by auto

ML_file "approximation_generator.ML"
setup "Approximation_Generator.setup"

end