Theory Affine_Arithmetic.Affine_Approximation

 section ‹Approximation with Affine Forms›
theory Affine_Approximation
imports
  "HOL-Decision_Procs.Approximation"
  "HOL-Library.Monad_Syntax"
  "HOL-Library.Mapping"
  Executable_Euclidean_Space
  Affine_Form
  Straight_Line_Program
begin
text ‹\label{sec:approxaffine}›

lemma convex_on_imp_above_tangent:― ‹TODO: generalizes @{thm convex_on_imp_above_tangent}
  assumes convex: "convex_on A f" and connected: "connected A"
  assumes c: "c  A" and x : "x  A"
  assumes deriv: "(f has_field_derivative f') (at c within A)"
  shows   "f x - f c  f' * (x - c)"
proof (cases x c rule: linorder_cases)
  assume xc: "x > c"
  let ?A' = "{c<..<x}"
  have subs: "?A'  A" using xc x c
    by (simp add: connected connected_contains_Ioo)
  have "at c within ?A'  bot"
    using xc
    by (simp add: at_within_eq_bot_iff)
  moreover from deriv have "((λy. (f y - f c) / (y - c))  f') (at c within ?A')"
    unfolding has_field_derivative_iff using subs
    by (blast intro: tendsto_mono at_le)
  moreover from eventually_at_right_real[OF xc]
    have "eventually (λy. (f y - f c) / (y - c)  (f x - f c) / (x - c)) (at_right c)"
  proof eventually_elim
    fix y assume y: "y  {c<..<x}"
    with convex connected x c have "f y  (f x - f c) / (x - c) * (y - c) + f c"
      using interior_subset[of A]
      by (intro convex_onD_Icc' convex_on_subset[OF convex] connected_contains_Icc) auto
    hence "f y - f c  (f x - f c) / (x - c) * (y - c)" by simp
    thus "(f y - f c) / (y - c)  (f x - f c) / (x - c)" using y xc by (simp add: divide_simps)
  qed
  hence "eventually (λy. (f y - f c) / (y - c)  (f x - f c) / (x - c)) (at c within ?A')"
    by (simp add: eventually_at_filter eventually_mono)
  ultimately have "f'  (f x - f c) / (x - c)" by (simp add: tendsto_upperbound)
  thus ?thesis using xc by (simp add: field_simps)
next
  assume xc: "x < c"
  let ?A' = "{x<..<c}"
  have subs: "?A'  A" using xc x c
    by (simp add: connected connected_contains_Ioo)
  have "at c within ?A'  bot"
    using xc
    by (simp add: at_within_eq_bot_iff)
  moreover from deriv have "((λy. (f y - f c) / (y - c))  f') (at c within ?A')"
    unfolding has_field_derivative_iff using subs
    by (blast intro: tendsto_mono at_le)
  moreover from eventually_at_left_real[OF xc]
    have "eventually (λy. (f y - f c) / (y - c)  (f x - f c) / (x - c)) (at_left c)"
  proof eventually_elim
    fix y assume y: "y  {x<..<c}"
    with convex connected x c have "f y  (f x - f c) / (c - x) * (c - y) + f c"
      using interior_subset[of A]
      by (intro convex_onD_Icc'' convex_on_subset[OF convex] connected_contains_Icc) auto
    hence "f y - f c  (f x - f c) * ((c - y) / (c - x))" by simp
    also have "(c - y) / (c - x) = (y - c) / (x - c)" using y xc by (simp add: field_simps)
    finally show "(f y - f c) / (y - c)  (f x - f c) / (x - c)" using y xc
      by (simp add: divide_simps)
  qed
  hence "eventually (λy. (f y - f c) / (y - c)  (f x - f c) / (x - c)) (at c within ?A')"
    by (simp add: eventually_at_filter eventually_mono)
  ultimately have "f'  (f x - f c) / (x - c)" by (simp add: tendsto_lowerbound)
  thus ?thesis using xc by (simp add: field_simps)
qed simp_all


text ‹Approximate operations on affine forms.›

lemma Affine_notempty[intro, simp]: "Affine X  {}"
  by (auto simp: Affine_def valuate_def)

lemma truncate_up_lt: "x < y  x < truncate_up prec y"
  by (rule less_le_trans[OF _ truncate_up])

lemma truncate_up_pos_eq[simp]: "0 < truncate_up p x  0 < x"
  by (auto simp: truncate_up_lt) (metis (poly_guards_query) not_le truncate_up_nonpos)

lemma inner_scaleR_pdevs_0: "inner_scaleR_pdevs 0 One_pdevs = zero_pdevs"
  unfolding inner_scaleR_pdevs_def
  by transfer (auto simp: unop_pdevs_raw_def)

lemma Affine_aform_of_point_eq[simp]: "Affine (aform_of_point p) = {p}"
  by (simp add: Affine_aform_of_ivl aform_of_point_def)

lemma mem_Affine_aform_of_point: "x  Affine (aform_of_point x)"
  by simp

lemma
  aform_val_aform_of_ivl_innerE:
  assumes "e  UNIV  {-1 .. 1}"
  assumes "a  b" "c  Basis"
  obtains f where "aform_val e (aform_of_ivl a b)  c = aform_val f (aform_of_ivl (a  c) (b  c))"
    "f  UNIV  {-1 .. 1}"
proof -
  have [simp]: "a  c  b  c"
    using assms by (auto simp: eucl_le[where 'a='a])
  have "(λx. x  c) ` Affine (aform_of_ivl a b) = Affine (aform_of_ivl (a  c) (b  c))"
    using assms
    by (auto simp: Affine_aform_of_ivl eucl_le[where 'a='a]
      image_eqI[where x="iBasis. (if i = c then x else a  i) *R i" for x])
  then obtain f where
      "aform_val e (aform_of_ivl a b)  c = aform_val f (aform_of_ivl (a  c) (b  c))"
      "f  UNIV  {-1 .. 1}"
    using assms
    by (force simp: Affine_def valuate_def)
  thus ?thesis ..
qed

lift_definition coord_pdevs::"nat  real pdevs" is "λn i. if i = n then 1 else 0" by auto

lemma pdevs_apply_coord_pdevs [simp]: "pdevs_apply (coord_pdevs i) x = (if x = i then 1 else 0)"
  by transfer simp

lemma degree_coord_pdevs[simp]: "degree (coord_pdevs i) = Suc i"
  by (auto intro!: degree_eqI)

lemma pdevs_val_coord_pdevs[simp]: "pdevs_val e (coord_pdevs i) = e i"
  by (auto simp: pdevs_val_sum if_distrib sum.delta cong: if_cong)

definition "aforms_of_ivls ls us = map
    (λ(i, (l, u)). ((l + u)/2, scaleR_pdevs ((u - l)/2) (coord_pdevs i)))
    (zip [0..<length ls] (zip ls us))"

lemma
  aforms_of_ivls:
  assumes "length ls = length us" "length xs = length ls"
  assumes "i. i < length xs  xs ! i  {ls ! i .. us ! i}"
  shows "xs  Joints (aforms_of_ivls ls us)"
proof -
  {
    fix i assume "i < length xs"
    then have "e. e  {-1 .. 1}  xs ! i = (ls ! i + us ! i) / 2 + e * (us ! i - ls ! i) / 2"
      using assms
      by (force intro!: exI[where x="(xs ! i - (ls ! i + us ! i) / 2) / (us ! i - ls ! i) * 2"]
          simp: divide_simps algebra_simps)
  } then obtain e where e: "e i  {-1 .. 1}"
    "xs ! i = (ls ! i + us ! i) / 2 + e i * (us ! i - ls ! i) / 2" 
    if "i < length xs" for i
    using that by metis
  define e' where "e' i = (if i < length xs then e i else 0)" for i
  show ?thesis
    using e assms
    by (auto simp: aforms_of_ivls_def Joints_def valuate_def e'_def aform_val_def
        intro!: image_eqI[where x=e'] nth_equalityI)
qed


subsection ‹Approximate Operations›

definition "max_pdev x = fold (λx y. if infnorm (snd x)  infnorm (snd y) then x else y) (list_of_pdevs x) (0, 0)"


subsubsection ‹set of generated endpoints›

fun points_of_list where
  "points_of_list x0 [] = [x0]"
| "points_of_list x0 ((i, x)#xs) = (points_of_list (x0 + x) xs @ points_of_list (x0 - x) xs)"

primrec points_of_aform where
  "points_of_aform (x, xs) = points_of_list x (list_of_pdevs xs)"


subsubsection ‹Approximate total deviation›

definition sum_list'::"nat  'a list  'a::executable_euclidean_space"
  where "sum_list' p xs = fold (λa b. eucl_truncate_up p (a + b)) xs 0"

definition "tdev' p x = sum_list' p (map (abs o snd) (list_of_pdevs x))"

lemma
  eucl_fold_mono:
  fixes f::"'a::ordered_euclidean_space'a'a"
  assumes mono: "w x y z. w  x  y  z  f w y  f x z"
  shows "x  y  fold f xs x  fold f xs y"
  by (induct xs arbitrary: x y) (auto simp: mono)

lemma sum_list_add_le_fold_eucl_truncate_up:
  fixes z::"'a::executable_euclidean_space"
  shows "sum_list xs + z  fold (λx y. eucl_truncate_up p (x + y)) xs z"
proof (induct xs arbitrary: z)
  case (Cons x xs)
  have "sum_list (x # xs) + z = sum_list xs + (z + x)"
    by simp
  also have "  fold (λx y. eucl_truncate_up p (x + y)) xs (z + x)"
    using Cons by simp
  also have "  fold (λx y. eucl_truncate_up p (x + y)) xs (eucl_truncate_up p (x + z))"
    by (auto intro!: add_mono eucl_fold_mono eucl_truncate_up eucl_truncate_up_mono simp: ac_simps)
  finally show ?case by simp
qed simp

lemma sum_list_le_sum_list':
  "sum_list xs  sum_list' p xs"
  unfolding sum_list'_def
  using sum_list_add_le_fold_eucl_truncate_up[of xs 0] by simp

lemma sum_list'_sum_list_le:
  "y  sum_list xs  y  sum_list' p xs"
  by (metis sum_list_le_sum_list' order.trans)

lemma tdev': "tdev x  tdev' p x"
  unfolding tdev'_def
proof -
  have "tdev x = (i = 0 ..< degree x. ¦pdevs_apply x i¦)"
    by (auto intro!: sum.mono_neutral_cong_left simp: tdev_def)
  also have " = (i  rev [0 ..< degree x]. ¦pdevs_apply x i¦)"
    by (metis atLeastLessThan_upt sum_list_rev rev_map sum_set_upt_conv_sum_list_nat)
  also have
    " = sum_list (map (λxa. ¦pdevs_apply x xa¦) [xarev [0..<degree x] . pdevs_apply x xa  0])"
    unfolding filter_map map_map o_def
    by (subst sum_list_map_filter) auto
  also note sum_list_le_sum_list'[of _ p]
  also have "[xarev [0..<degree x] . pdevs_apply x xa  0] =
      rev (sorted_list_of_set (pdevs_domain x))"
    by (subst rev_is_rev_conv[symmetric])
      (auto simp: filter_map rev_filter intro!: sorted_distinct_set_unique
        sorted_filter[of "λx. x", simplified] degree_gt)
  finally
  show "tdev x  sum_list' p (map (abs  snd) (list_of_pdevs x))"
    by (auto simp: list_of_pdevs_def o_def rev_map filter_map rev_filter)
qed

lemma tdev'_le: "x  tdev y  x  tdev' p y"
  by (metis order.trans tdev')

lemmas abs_pdevs_val_le_tdev' = tdev'_le[OF abs_pdevs_val_le_tdev]

lemma tdev'_uminus_pdevs[simp]: "tdev' p (uminus_pdevs x) = tdev' p x"
  by (auto simp: tdev'_def o_def rev_map filter_map rev_filter list_of_pdevs_def pdevs_domain_def)

abbreviation Radius::"'a::ordered_euclidean_space aform  'a"
  where "Radius X  tdev (snd X)"

abbreviation Radius'::"nat'a::executable_euclidean_space aform  'a"
  where "Radius' p X  tdev' p (snd X)"

lemma Radius'_uminus_aform[simp]: "Radius' p (uminus_aform X) = Radius' p X"
  by (auto simp: uminus_aform_def)


subsubsection ‹truncate partial deviations›

definition trunc_pdevs_raw::"nat  (nat  'a)  nat  'a::executable_euclidean_space"
  where "trunc_pdevs_raw p x i = eucl_truncate_down p (x i)"

lemma nonzeros_trunc_pdevs_raw:
  "{i. trunc_pdevs_raw r x i  0}  {i. x i  0}"
  by (auto simp: trunc_pdevs_raw_def[abs_def])

lift_definition trunc_pdevs::"nat  'a::executable_euclidean_space pdevs  'a pdevs"
  is trunc_pdevs_raw
  by (auto intro!: finite_subset[OF nonzeros_trunc_pdevs_raw])

definition trunc_err_pdevs_raw::"nat  (nat  'a)  nat  'a::executable_euclidean_space"
  where "trunc_err_pdevs_raw p x i = trunc_pdevs_raw p x i - x i"

lemma nonzeros_trunc_err_pdevs_raw:
  "{i. trunc_err_pdevs_raw r x i  0}  {i. x i  0}"
  by (auto simp: trunc_pdevs_raw_def trunc_err_pdevs_raw_def[abs_def])

lift_definition trunc_err_pdevs::"nat  'a::executable_euclidean_space pdevs  'a pdevs"
  is trunc_err_pdevs_raw
  by (auto intro!: finite_subset[OF nonzeros_trunc_err_pdevs_raw])

term float_plus_down

lemma pdevs_apply_trunc_pdevs[simp]:
  fixes x y::"'a::euclidean_space"
  shows "pdevs_apply (trunc_pdevs p X) n = eucl_truncate_down p (pdevs_apply X n)"
  by transfer (simp add: trunc_pdevs_raw_def)

lemma pdevs_apply_trunc_err_pdevs[simp]:
  fixes x y::"'a::euclidean_space"
  shows "pdevs_apply (trunc_err_pdevs p X) n =
    eucl_truncate_down p (pdevs_apply X n) - (pdevs_apply X n)"
  by transfer (auto simp: trunc_err_pdevs_raw_def trunc_pdevs_raw_def)

lemma pdevs_val_trunc_pdevs:
  fixes x y::"'a::euclidean_space"
  shows "pdevs_val e (trunc_pdevs p X) = pdevs_val e X + pdevs_val e (trunc_err_pdevs p X)"
proof -
  have "pdevs_val e X + pdevs_val e (trunc_err_pdevs p X) =
      pdevs_val e (add_pdevs X (trunc_err_pdevs p X))"
    by simp
  also have " = pdevs_val e (trunc_pdevs p X)"
    by (auto simp: pdevs_val_def trunc_pdevs_raw_def trunc_err_pdevs_raw_def)
  finally show ?thesis by simp
qed

lemma pdevs_val_trunc_err_pdevs:
  fixes x y::"'a::euclidean_space"
  shows "pdevs_val e (trunc_err_pdevs p X) = pdevs_val e (trunc_pdevs p X) - pdevs_val e X"
  by (simp add: pdevs_val_trunc_pdevs)

definition truncate_aform::"nat  'a aform  'a::executable_euclidean_space aform"
  where "truncate_aform p x = (eucl_truncate_down p (fst x), trunc_pdevs p (snd x))"

definition truncate_error_aform::"nat  'a aform  'a::executable_euclidean_space aform"
  where "truncate_error_aform p x =
    (eucl_truncate_down p (fst x) - fst x, trunc_err_pdevs p (snd x))"

lemma
  abs_aform_val_le:
  assumes "e  UNIV  {- 1..1}"
  shows "abs (aform_val e X)  eucl_truncate_up p (¦fst X¦ + tdev' p (snd X))"
proof -
  have "abs (aform_val e X)  ¦fst X¦ + ¦pdevs_val e (snd X)¦"
    by (auto simp: aform_val_def intro!: abs_triangle_ineq)
  also have "¦pdevs_val e (snd X)¦  tdev (snd X)"
    using assms by (rule abs_pdevs_val_le_tdev)
  also note tdev'
  also note eucl_truncate_up
  finally show ?thesis by simp
qed


subsubsection ‹truncation with error bound›

definition "trunc_bound_eucl p s =
  (let
    d = eucl_truncate_down p s;
    ed = abs (d - s) in
  (d, eucl_truncate_up p ed))"

lemma trunc_bound_euclE:
  obtains err where
  "¦err¦  snd (trunc_bound_eucl p x)"
  "fst (trunc_bound_eucl p x) = x + err"
proof atomize_elim
  have "fst (trunc_bound_eucl p x) = x + (eucl_truncate_down p x - x)"
    (is "_ = _ + ?err")
    by (simp_all add: trunc_bound_eucl_def Let_def)
  moreover have "abs ?err  snd (trunc_bound_eucl p x)"
    by (simp add: trunc_bound_eucl_def Let_def eucl_truncate_up)
  ultimately show "err. ¦err¦  snd (trunc_bound_eucl p x)  fst (trunc_bound_eucl p x) = x + err"
    by auto
qed

definition "trunc_bound_pdevs p x = (trunc_pdevs p x, tdev' p (trunc_err_pdevs p x))"

lemma pdevs_apply_fst_trunc_bound_pdevs[simp]: "pdevs_apply (fst (trunc_bound_pdevs p x)) =
  pdevs_apply (trunc_pdevs p x)"
  by (simp add: trunc_bound_pdevs_def)

lemma trunc_bound_pdevsE:
  assumes "e  UNIV  {- 1..1}"
  obtains err where
  "¦err¦  snd (trunc_bound_pdevs p x)"
  "pdevs_val e (fst ((trunc_bound_pdevs p x))) = pdevs_val e x + err"
proof atomize_elim
  have "pdevs_val e (fst (trunc_bound_pdevs p x)) = pdevs_val e x +
    pdevs_val e (add_pdevs (trunc_pdevs p x) (uminus_pdevs x))"
    (is "_ = _ + ?err")
    by (simp_all add: trunc_bound_pdevs_def Let_def)
  moreover have "abs ?err  snd (trunc_bound_pdevs p x)"
    using assms
    by (auto simp add: pdevs_val_trunc_pdevs trunc_bound_pdevs_def Let_def eucl_truncate_up
      intro!: order_trans[OF abs_pdevs_val_le_tdev tdev'])
  ultimately show "err. ¦err¦  snd (trunc_bound_pdevs p x) 
      pdevs_val e (fst ((trunc_bound_pdevs p x))) = pdevs_val e x + err"
    by auto
qed

lemma
  degree_add_pdevs_le:
  assumes "degree X  n"
  assumes "degree Y  n"
  shows "degree (add_pdevs X Y)  n"
  using assms
  by (auto intro!: degree_le)


lemma truncate_aform_error_aform_cancel:
  "aform_val e (truncate_aform p z) = aform_val e z + aform_val e (truncate_error_aform p z) "
  by (simp add: truncate_aform_def aform_val_def truncate_error_aform_def pdevs_val_trunc_pdevs)

lemma error_absE:
  assumes "abs err  k"
  obtains e::real where "err = e * k" "e  {-1 .. 1}"
  using assms
  by atomize_elim
    (safe intro!: exI[where x="err / abs k"] divide_atLeastAtMost_1_absI, auto)

lemma eucl_truncate_up_nonneg_eq_zero_iff:
  "x  0  eucl_truncate_up p x = 0  x = 0"
  by (metis (poly_guards_query) eq_iff eucl_truncate_up eucl_truncate_up_zero)

lemma
  aform_val_consume_error:
  assumes "abs err  abs (pdevs_apply (snd X) n)"
  shows "aform_val (e(n := 0)) X + err = aform_val (e(n := err/pdevs_apply (snd X) n)) X"
  using assms
  by (auto simp add: aform_val_def)

lemma
  aform_val_consume_errorE:
  fixes X::"real aform"
  assumes "abs err  abs (pdevs_apply (snd X) n)"
  obtains err' where "aform_val (e(n := 0)) X + err = aform_val (e(n := err')) X" "err'  {-1 .. 1}"
  by atomize_elim
    (rule aform_val_consume_error assms aform_val_consume_error exI conjI
      divide_atLeastAtMost_1_absI)+

lemma
  degree_trunc_pdevs_le:
  assumes "degree X  n"
  shows "degree (trunc_pdevs p X)  n"
  using assms
  by (auto intro!: degree_le)

lemma pdevs_val_sum_less_degree:
  "pdevs_val e X = (i<d. e i *R pdevs_apply X i)" if "degree X  d"
  unfolding pdevs_val_pdevs_domain
  apply (rule sum.mono_neutral_cong_left)
  using that
  by force+


subsubsection ‹general affine operation›

definition "affine_binop (X::real aform) Y a b c d k =
  (a * fst X + b * fst Y + c,
    pdev_upd (add_pdevs (scaleR_pdevs a (snd X)) (scaleR_pdevs b (snd Y))) k d)"

lemma pdevs_domain_One_pdevs[simp]: "pdevs_domain (One_pdevs::'a::executable_euclidean_space pdevs) =
  {0..<DIM('a)}"
  apply (auto simp: length_Basis_list split: if_splits)
  subgoal for i
    using nth_Basis_list_in_Basis[of i, where 'a='a]
    by (auto simp: length_Basis_list)
  done

lemma pdevs_val_One_pdevs:
  "pdevs_val e (One_pdevs::'a::executable_euclidean_space pdevs) = (i<DIM('a). e i *R Basis_list ! i)"
  by (auto simp: pdevs_val_pdevs_domain length_Basis_list intro!:sum.cong)

lemma affine_binop:
  assumes "degree_aforms [X, Y]  k"
  shows "aform_val e (affine_binop X Y a b c d k) =
    a * aform_val e X + b * aform_val e Y + c + e k * d"
  using assms
  by (auto simp: aform_val_def affine_binop_def degrees_def
      pdevs_val_msum_pdevs degree_add_pdevs_le pdevs_val_One_pdevs Basis_list_real_def
      algebra_simps)

definition "affine_binop' p (X::real aform) Y a b c d k =
  (let
    ― ‹TODO: more round-off operations here?›
    (r, e1) = trunc_bound_eucl p (a * fst X + b * fst Y + c);
    (Z, e2) = trunc_bound_pdevs p (add_pdevs (scaleR_pdevs a (snd X)) (scaleR_pdevs b (snd Y)))
  in
    (r, pdev_upd Z k (sum_list' p [e1, e2, d]))
  )"

lemma sum_list'_noneg_eq_zero_iff: "sum_list' p xs = 0  (xset xs. x = 0)" if "x. x  set xs  x  0"
proof safe
  fix x assume x: "sum_list' p xs = 0" "x  set xs"
  from that have "0  sum_list xs" by (auto intro!: sum_list_nonneg)
  with that x have "sum_list xs = 0"
    by (metis antisym sum_list_le_sum_list')
  then have "(i<length xs.  xs ! i) = 0"
    by (auto simp: sum_list_sum_nth atLeast0LessThan)
  then show "x = 0" using x(2) that
    by (subst (asm) sum_nonneg_eq_0_iff) (auto simp: in_set_conv_nth)
next
  show "xset xs. x = 0  sum_list' p xs = 0"
    by (induction xs) (auto simp: sum_list'_def)
qed

lemma affine_binop'E:
  assumes deg: "degree_aforms [X, Y]  k"
  assumes e: "e  UNIV  {- 1..1}"
  assumes d: "abs u  d"
  obtains ek where
    "a * aform_val e X + b * aform_val e Y + c + u = aform_val (e(k:=ek)) (affine_binop' p X Y a b c d k)"
    "ek  {-1 .. 1}"
proof -
  have "a * aform_val e X + b * aform_val e Y + c + u =
    (a * fst X + b * fst Y + c) + pdevs_val e (add_pdevs (scaleR_pdevs a (snd X)) (scaleR_pdevs b (snd Y))) + u"
    (is "_ = ?c + pdevs_val _ ?ps + _")
    by (auto simp: aform_val_def algebra_simps)

  from trunc_bound_euclE[of p ?c] obtain ec where ec: "abs ec  snd (trunc_bound_eucl p ?c)"
    "fst (trunc_bound_eucl p ?c) - ec = ?c"
    by (auto simp: algebra_simps)

  moreover

  from trunc_bound_pdevsE[OF e, of p ?ps]
  obtain eps where eps: "¦eps¦  snd (trunc_bound_pdevs p ?ps)"
    "pdevs_val e (fst (trunc_bound_pdevs p ?ps)) - eps = pdevs_val e ?ps"
    by (auto simp: algebra_simps)

  moreover
  define ek where "ek = (u - ec - eps)/
        sum_list' p [snd (trunc_bound_eucl p ?c), snd (trunc_bound_pdevs p ?ps), d]"
  have "degree (fst (trunc_bound_pdevs p ?ps)) 
      degree_aforms [X, Y]"
    by (auto simp: trunc_bound_pdevs_def degrees_def intro!: degree_trunc_pdevs_le degree_add_pdevs_le)
  moreover
  from this have "pdevs_apply (fst (trunc_bound_pdevs p ?ps)) k = 0"
    using deg order_trans by blast
  ultimately have "a * aform_val e X + b * aform_val e Y + c + u =
    aform_val (e(k:=ek)) (affine_binop' p X Y a b c d k)"
    apply (auto simp: affine_binop'_def algebra_simps aform_val_def split: prod.splits)
    subgoal for x y z
      apply (cases "sum_list' p [x, z, d] = 0")
      subgoal
        apply simp
        apply (subst (asm) sum_list'_noneg_eq_zero_iff)
        using d deg
        by auto
      subgoal
        apply (simp add: divide_simps algebra_simps ek_def)
        using pdevs_apply (fst (trunc_bound_pdevs p (add_pdevs (scaleR_pdevs a (snd X)) (scaleR_pdevs b (snd Y))))) k = 0 by auto
      done
    done
  moreover have "ek  {-1 .. 1}"
    unfolding ek_def
    apply (rule divide_atLeastAtMost_1_absI)
    apply (rule abs_triangle_ineq4[THEN order_trans])
    apply (rule order_trans)
     apply (rule add_right_mono)
     apply (rule abs_triangle_ineq4)
    using ec(1) eps(1)
    by (auto simp: sum_list'_def eucl_truncate_up_real_def add.assoc
        intro!: order_trans[OF _ abs_ge_self] order_trans[OF _ truncate_up_le] add_mono d )
  ultimately show ?thesis ..
qed

subsubsection ‹Inf/Sup›

definition "Inf_aform' p X = eucl_truncate_down p (fst X - tdev' p (snd X))"

definition "Sup_aform' p X = eucl_truncate_up p (fst X + tdev' p (snd X))"

lemma Inf_aform':
  shows "Inf_aform' p X  Inf_aform X"
  unfolding Inf_aform_def Inf_aform'_def
  by (auto intro!: eucl_truncate_down_le add_left_mono tdev')

lemma Sup_aform':
  shows "Sup_aform X  Sup_aform' p X"
  unfolding Sup_aform_def Sup_aform'_def
  by (rule eucl_truncate_up_le add_left_mono tdev')+

lemma Inf_aform_le_Sup_aform[intro]:
  "Inf_aform X  Sup_aform X"
  by (simp add: Inf_aform_def Sup_aform_def algebra_simps)

lemma Inf_aform'_le_Sup_aform'[intro]:
  "Inf_aform' p X  Sup_aform' p X"
  by (metis Inf_aform' Inf_aform_le_Sup_aform Sup_aform' order.trans)

definition
  "ivls_of_aforms prec = map (λa. Interval' (float_of (Inf_aform' prec a)) (float_of(Sup_aform' prec a)))"

lemma
  assumes "i. e'' i  1"
  assumes "i. -1  e'' i"
  shows Inf_aform'_le: "Inf_aform' p r  aform_val e'' r"
    and Sup_aform'_le: "aform_val e'' r  Sup_aform' p r"
  by (auto intro!: order_trans[OF Inf_aform'] order_trans[OF _ Sup_aform'] Inf_aform Sup_aform
    simp: Affine_def valuate_def intro!: image_eqI[where x=e''] assms)


lemma InfSup_aform'_in_float[intro, simp]:
  "Inf_aform' p X  float" "Sup_aform' p X  float"
  by (auto simp: Inf_aform'_def eucl_truncate_down_real_def
      Sup_aform'_def eucl_truncate_up_real_def)

theorem ivls_of_aforms: "xs  Joints XS  bounded_by xs (ivls_of_aforms prec XS)"
  by (auto simp: bounded_by_def ivls_of_aforms_def Affine_def valuate_def Pi_iff set_of_eq
      intro!: Inf_aform'_le Sup_aform'_le
      dest!: nth_in_AffineI split: Interval'_splits)

definition "isFDERIV_aform prec N xs fas AS = isFDERIV_approx prec N xs fas (ivls_of_aforms prec AS)"

theorem isFDERIV_aform:
  assumes "isFDERIV_aform prec N xs fas AS"
  assumes "vs  Joints AS"
  shows "isFDERIV N xs fas vs"
  apply (rule isFDERIV_approx)
  apply (rule ivls_of_aforms)
  apply (rule assms)
  apply (rule assms[unfolded isFDERIV_aform_def])
  done

definition "env_len env l = (xs  env. length xs = l)"

lemma env_len_takeI: "env_len xs d1  d1  d  env_len (take d ` xs) d"
  by (auto simp: env_len_def)

subsection ‹Min Range approximation›

lemma
  linear_lower:
  fixes x::real
  assumes "x. x  {a .. b}  (f has_field_derivative f' x) (at x within {a .. b})"
  assumes "x. x  {a .. b}  f' x  u"
  assumes "x  {a .. b}"
  shows "f b + u * (x - b)  f x"
proof -
  from assms(2-)
    mvt_very_simple[of x b f "λx. (*) (f' x)",
      rule_format,
      OF _ has_derivative_subset[OF assms(1)[simplified has_field_derivative_def]]]
  obtain y where "y  {x .. b}"  "f b - f x = (b - x) * f' y"
    by (auto simp: Bex_def ac_simps)
  moreover hence "f' y  u" using assms by auto
  ultimately have "f b - f x  (b - x) * u"
    by (auto intro!: mult_left_mono)
  thus ?thesis by (simp add: algebra_simps)
qed

lemma
  linear_lower2:
  fixes x::real
  assumes "x. x  {a .. b}  (f has_field_derivative f' x) (at x within {a .. b})"
  assumes "x. x  {a .. b}  l  f' x"
  assumes "x  {a .. b}"
  shows "f x  f a + l * (x - a)"
proof -
  from assms(2-)
    mvt_very_simple[of a x f "λx. (*) (f' x)",
      rule_format,
      OF _ has_derivative_subset[OF assms(1)[simplified has_field_derivative_def]]]
  obtain y where "y  {a .. x}"  "f x - f a = (x - a) * f' y"
    by (auto simp: Bex_def ac_simps)
  moreover hence "l  f' y" using assms by auto
  ultimately have "(x - a) * l  f x - f a"
    by (auto intro!: mult_left_mono)
  thus ?thesis by (simp add: algebra_simps)
qed

lemma
  linear_upper:
  fixes x::real
  assumes "x. x  {a .. b}  (f has_field_derivative f' x) (at x within {a .. b})"
  assumes "x. x  {a .. b}  f' x  u"
  assumes "x  {a .. b}"
  shows "f x  f a + u * (x - a)"
proof -
  from assms(2-)
    mvt_very_simple[of a x f "λx. (*) (f' x)",
      rule_format,
      OF _ has_derivative_subset[OF assms(1)[simplified has_field_derivative_def]]]
  obtain y where "y  {a .. x}"  "f x - f a = (x - a) * f' y"
    by (auto simp: Bex_def ac_simps)
  moreover hence "f' y  u" using assms by auto
  ultimately have "(x - a) * u  f x - f a"
    by (auto intro!: mult_left_mono)
  thus ?thesis by (simp add: algebra_simps)
qed

lemma
  linear_upper2:
  fixes x::real
  assumes "x. x  {a .. b}  (f has_field_derivative f' x) (at x within {a .. b})"
  assumes "x. x  {a .. b}  l  f' x"
  assumes "x  {a .. b}"
  shows "f x  f b + l * (x - b)"
proof -
  from assms(2-)
    mvt_very_simple[of x b f "λx. (*) (f' x)",
      rule_format,
      OF _ has_derivative_subset[OF assms(1)[simplified has_field_derivative_def]]]
  obtain y where "y  {x .. b}"  "f b - f x = (b - x) * f' y"
    by (auto simp: Bex_def ac_simps)
  moreover hence "l  f' y" using assms by auto
  ultimately have "f b - f x  (b - x) * l"
    by (auto intro!: mult_left_mono)
  thus ?thesis by (simp add: algebra_simps)
qed

lemma linear_enclosure:
  fixes x::real
  assumes "x. x  {a .. b}  (f has_field_derivative f' x) (at x within {a .. b})"
  assumes "x. x  {a .. b}  f' x  u"
  assumes "x  {a .. b}"
  shows "f x  {f b + u * (x - b) .. f a + u * (x - a)}"
  using linear_lower[OF assms] linear_upper[OF assms]
  by auto

definition "mid_err ivl = ((lower ivl + upper ivl::float)/2, (upper ivl - lower ivl)/2)"

lemma degree_aform_uminus_aform[simp]: "degree_aform (uminus_aform X) = degree_aform X"
  by (auto simp: uminus_aform_def)


subsubsection ‹Addition›

definition add_aform::"'a::real_vector aform  'a aform  'a aform"
  where "add_aform x y = (fst x + fst y, add_pdevs (snd x) (snd y))"

lemma aform_val_add_aform:
  shows "aform_val e (add_aform X Y) = aform_val e X + aform_val e Y"
  by (auto simp: add_aform_def aform_val_def)

type_synonym aform_err = "real aform × real"

definition add_aform'::"nat  aform_err  aform_err  aform_err"
  where "add_aform' p x y =
    (let
      z0 = trunc_bound_eucl p (fst (fst x) + fst (fst y));
      z = trunc_bound_pdevs p (add_pdevs (snd (fst x)) (snd (fst y)))
      in ((fst z0, fst z), (sum_list' p [snd z0, snd z, abs (snd x), abs (snd y)])))"

abbreviation degree_aform_err::"aform_err  nat"
  where "degree_aform_err X  degree_aform (fst X)"

lemma degree_aform_err_add_aform':
  assumes "degree_aform_err x  n"
  assumes "degree_aform_err y  n"
  shows "degree_aform_err (add_aform' p x y)  n"
  using assms
  by (auto simp: add_aform'_def Let_def trunc_bound_pdevs_def
      intro!: degree_pdev_upd_le degree_trunc_pdevs_le degree_add_pdevs_le)

definition "aform_err e Xe = {aform_val e (fst Xe) - snd Xe .. aform_val e (fst Xe) + snd Xe::real}"

lemma aform_errI: "x  aform_err e Xe"
  if "abs (x - aform_val e (fst Xe))  snd Xe"
  using that by (auto simp: aform_err_def abs_real_def algebra_simps split: if_splits)

lemma add_aform':
  assumes e: "e  UNIV  {- 1..1}"
  assumes x: "x  aform_err e X"
  assumes y: "y  aform_err e Y"
  shows "x + y  aform_err e (add_aform' p X Y)"
proof -
  let ?t1 = "trunc_bound_eucl p (fst (fst X) + fst (fst Y))"
  from trunc_bound_euclE
  obtain e1 where abs_e1: "¦e1¦  snd ?t1" and e1: "fst ?t1 = fst (fst X) + fst (fst Y) + e1"
    by blast
  let ?t2 = "trunc_bound_pdevs p (add_pdevs (snd (fst X)) (snd (fst Y)))"
  from trunc_bound_pdevsE[OF e, of p "add_pdevs (snd (fst X)) (snd (fst Y))"]
  obtain e2 where abs_e2: "¦e2¦  snd (?t2)"
    and e2: "pdevs_val e (fst ?t2) = pdevs_val e (add_pdevs (snd (fst X)) (snd (fst Y))) + e2"
    by blast

  have e_le: "¦e1 + e2 + snd X + snd Y¦  snd (add_aform' p (X) Y)"
    apply (auto simp: add_aform'_def Let_def )
    apply (rule sum_list'_sum_list_le)
    apply (simp add: add.assoc)
    by (intro order.trans[OF abs_triangle_ineq] add_mono abs_e1 abs_e2 order_refl)
  then show ?thesis
    apply (intro aform_errI)
    using x y abs_e1 abs_e2
    apply (simp add: aform_val_def aform_err_def add_aform_def add_aform'_def Let_def e1 e2 assms)
    by (auto intro!: order_trans[OF _ sum_list_le_sum_list'] )
qed


subsubsection ‹Scaling›

definition aform_scaleR::"real aform  'a::real_vector  'a aform"
  where "aform_scaleR x y = (fst x *R y, pdevs_scaleR (snd x) y)"

lemma aform_val_scaleR_aform[simp]:
  shows "aform_val e (aform_scaleR X y) = aform_val e X *R y"
  by (auto simp: aform_scaleR_def aform_val_def scaleR_left_distrib)


subsubsection ‹Multiplication›

lemma aform_val_mult_exact:
  "aform_val e x * aform_val e y =
    fst x * fst y +
    pdevs_val e (add_pdevs (scaleR_pdevs (fst y) (snd x)) (scaleR_pdevs (fst x) (snd y))) +
    (i<d. e i *R pdevs_apply (snd x) i)*(i<d. e i *R pdevs_apply (snd y) i)"
   if "degree (snd x)  d" "degree (snd y)  d"
   using that
  by (auto simp: pdevs_val_sum_less_degree[where d=d] aform_val_def algebra_simps)

lemma sum_times_bound:― ‹TODO: this gives better bounds for the remainder of multiplication›
  "(i<d. e i * f i::real) * (i<d. e i * g i) =
   (i<d. (e i)2 * (f i * g i)) +
   ((i, j) | i < j  j < d. (e i * e j) * (f j * g i + f i * g j))" for d::nat
proof -
  have "(i<d. e i * f i)*(i<d. e i * g i) = ((i, j){..<d} × {..<d}. e i * f i * (e j * g j))"
    unfolding sum_product sum.cartesian_product ..
  also have " = ((i, j){..<d} × {..<d}  {(i, j). i = j}. e i * f i * (e j * g j)) +
    (((i, j){..<d} × {..<d}  {(i, j). i < j}. e i * f i * (e j * g j)) +
    ((i, j){..<d} × {..<d}  {(i, j). j < i}. e i * f i * (e j * g j)))"
    (is "_ = ?a + (?b + ?c)")
    by (subst sum.union_disjoint[symmetric], force, force, force)+ (auto intro!: sum.cong)
  also have "?c = ((i, j){..<d} × {..<d}  {(i, j). i < j}. e i * f j * (e j * g i))"
    by (rule sum.reindex_cong[of "λ(x, y). (y, x)"]) (auto intro!: inj_onI)
  also have "?b +  = ((i, j){..<d} × {..<d}  {(i, j). i < j}. (e i * e j) * (f j * g i + f i * g j))"
    by (auto simp: algebra_simps sum.distrib split_beta')
  also have " = ((i, j) | i < j  j < d. (e i * e j) * (f j * g i + f i * g j))"
    by (rule sum.cong) auto
  also have "?a = (i<d. (e i)2 * (f i * g i))"
    by (rule sum.reindex_cong[of "λi. (i, i)"]) (auto simp: power2_eq_square intro!: inj_onI)
  finally show ?thesis by simp
qed

definition mult_aform::"aform_err  aform_err  aform_err"
  where "mult_aform x y = ((fst (fst x) * fst (fst y),
    (add_pdevs (scaleR_pdevs (fst (fst y)) (snd (fst x))) (scaleR_pdevs (fst (fst x)) (snd (fst y))))),
     (tdev (snd (fst x)) * tdev (snd (fst y)) +
      abs (snd x) * (abs (fst (fst y)) + Radius (fst y)) +
      abs (snd y) * (abs (fst (fst x)) + Radius (fst x)) + abs (snd x) * abs (snd y)
     ))"

lemma mult_aformE:
  fixes X Y::"aform_err"
  assumes e: "e  UNIV  {- 1..1}"
  assumes x: "x  aform_err e X"
  assumes y: "y  aform_err e Y"
  shows "x * y  aform_err e (mult_aform X Y)"
proof -
  define ex where "ex  x - aform_val e (fst X)"
  define ey where "ey  y - aform_val e (fst Y)"

  have [intro, simp]: "¦ex¦  ¦snd X¦" "¦ey¦  ¦snd Y¦"
    using x y
    by (auto simp: ex_def ey_def aform_err_def)
  have "x * y =
    fst (fst X) * fst (fst Y) +
    fst (fst Y) * pdevs_val e (snd (fst X)) +
    fst (fst X) * pdevs_val e (snd (fst Y)) +

    (pdevs_val e (snd (fst X)) * pdevs_val e (snd (fst Y)) +
    ex * (fst (fst Y) + pdevs_val e (snd (fst Y))) +
    ey * (fst (fst X) + pdevs_val e (snd (fst X))) +
    ex * ey)"
    (is "_ = ?c + ?d + ?e + ?err")
    by (auto simp: ex_def ey_def algebra_simps aform_val_def)

  have abs_err: "abs ?err  snd (mult_aform X Y)"
    by (auto simp: mult_aform_def abs_mult
        intro!: abs_triangle_ineq[THEN order_trans] add_mono mult_mono
          abs_pdevs_val_le_tdev e)
  show ?thesis
    apply (auto simp: intro!: aform_errI order_trans[OF _ abs_err])
    apply (subst mult_aform_def)
    apply (auto simp: aform_val_def ex_def ey_def algebra_simps)
    done
qed

definition mult_aform'::"nat  aform_err  aform_err  aform_err"
  where "mult_aform' p x y = (
    let
      (fx, sx) = x;
      (fy, sy) = y;
      ex = abs sx;
      ey = abs sy;
      z0 = trunc_bound_eucl p (fst fx * fst fy);
      u = trunc_bound_pdevs p (scaleR_pdevs (fst fy) (snd fx));
      v = trunc_bound_pdevs p (scaleR_pdevs (fst fx) (snd fy));
      w = trunc_bound_pdevs p (add_pdevs (fst u) (fst v));
      tx = tdev' p (snd fx);
      ty = tdev' p (snd fy);
      l = truncate_up p (tx * ty);
      ee = truncate_up p (ex * ey);
      e1 = truncate_up p (ex * truncate_up p (abs (fst fy) + ty));
      e2 = truncate_up p (ey * truncate_up p (abs (fst fx) + tx))
    in
      ((fst z0, (fst w)), (sum_list' p [ee, e1, e2, l, snd z0, snd u, snd v, snd w])))"

lemma aform_errE:
  "abs (x - aform_val e (fst X))  snd X"
  if "x  aform_err e X"
  using that by (auto simp: aform_err_def)

lemma mult_aform'E:
  fixes X Y::"aform_err"
  assumes e: "e  UNIV  {- 1..1}"
  assumes x: "x  aform_err e X"
  assumes y: "y  aform_err e Y"
  shows "x * y  aform_err e (mult_aform' p X Y)"
proof -
  let ?z0 = "trunc_bound_eucl p (fst (fst X) * fst (fst Y))"
  from trunc_bound_euclE
  obtain e1 where abs_e1: "¦e1¦  snd ?z0" and e1: "fst ?z0 = fst (fst X) * fst (fst Y) + e1"
    by blast
  let ?u = "trunc_bound_pdevs p (scaleR_pdevs (fst (fst Y)) (snd (fst X)))"
  from trunc_bound_pdevsE[OF e]
  obtain e2 where abs_e2: "¦e2¦  snd (?u)"
    and e2: "pdevs_val e (fst ?u) = pdevs_val e (scaleR_pdevs (fst (fst Y)) (snd (fst X))) + e2"
    by blast
  let ?v = "trunc_bound_pdevs p (scaleR_pdevs (fst (fst X)) (snd (fst Y)))"
  from trunc_bound_pdevsE[OF e]
  obtain e3 where abs_e3: "¦e3¦  snd (?v)"
    and e3: "pdevs_val e (fst ?v) = pdevs_val e (scaleR_pdevs (fst (fst X)) (snd (fst Y))) + e3"
    by blast
  let ?w = "trunc_bound_pdevs p (add_pdevs (fst ?u) (fst ?v))"
  from trunc_bound_pdevsE[OF e]
  obtain e4 where abs_e4: "¦e4¦  snd (?w)"
    and e4: "pdevs_val e (fst ?w) = pdevs_val e (add_pdevs (fst ?u) (fst ?v)) + e4"
    by blast
  let ?tx = "tdev' p (snd (fst X))" and ?ty = "tdev' p (snd (fst Y))"
  let ?l = "truncate_up p (?tx * ?ty)"
  let ?ee = "truncate_up p (abs (snd X) * abs (snd Y))"
  let ?e1 = "truncate_up p (abs (snd X) * truncate_up p (¦fst (fst Y)¦ + ?ty))"
  let ?e2 = "truncate_up p (abs (snd Y) * truncate_up p (¦fst (fst X)¦ + ?tx))"

  let ?e0 = "x * y - fst (fst X) * fst (fst Y) -
      fst (fst X) * pdevs_val e (snd (fst Y)) -
      fst (fst Y) * pdevs_val e (snd (fst X))"
  let ?err = "?e0 - (e1 + e2  + e3 + e4)"
  have "abs ?err  abs ?e0 + abs e1 + abs e2 + abs e3 + abs e4"
    by arith
  also have "  abs ?e0 + snd ?z0 + snd ?u + snd ?v + snd ?w"
    unfolding abs_mult
    by (auto intro!: add_mono mult_mono e abs_pdevs_val_le_tdev' abs_ge_zero abs_e1 abs_e2 abs_e3
      abs_e4 intro: tdev'_le)
  also
  have asdf: "snd (mult_aform X Y)  tdev' p (snd (fst X)) * tdev' p (snd (fst Y)) + ?e1 + ?e2 + ?ee"
    by (auto simp: mult_aform_def intro!: add_mono mult_mono order_trans[OF _ tdev'] truncate_up_le)
  have "abs ?e0  ?ee + ?e1 + ?e2 + tdev' p (snd (fst X)) * tdev' p (snd (fst Y))"
    using mult_aformE[OF e x y, THEN aform_errE, THEN order_trans, OF asdf]
    by (simp add: aform_val_def mult_aform_def) arith
  also have "tdev' p (snd (fst X)) * tdev' p (snd (fst Y))  ?l"
    by (auto intro!: truncate_up_le)
  also have "?ee + ?e1 + ?e2 + ?l + snd ?z0 + snd ?u + snd ?v + snd ?w 
      sum_list' p [?ee, ?e1, ?e2, ?l, snd ?z0, snd ?u, snd ?v, snd ?w]"
    by (rule order_trans[OF _ sum_list_le_sum_list']) simp
  also have "  (snd (mult_aform' p X Y))"
    by (auto simp: mult_aform'_def Let_def assms split: prod.splits)
  finally have err_le: "abs ?err  (snd (mult_aform' p X Y))" by arith

  show ?thesis
    apply (rule aform_errI[OF order_trans[OF _ err_le]])
    apply (subst mult_aform'_def)
    using e1 e2 e3 e4
    apply (auto simp: aform_val_def Let_def assms split: prod.splits)
    done
qed

lemma degree_aform_mult_aform':
  assumes "degree_aform_err x  n"
  assumes "degree_aform_err y  n"
  shows "degree_aform_err (mult_aform' p x y)  n"
  using assms
  by (auto simp: mult_aform'_def Let_def trunc_bound_pdevs_def split: prod.splits
      intro!: degree_pdev_upd_le degree_trunc_pdevs_le degree_add_pdevs_le)

lemma
  fixes x a b::real
  assumes "a > 0"
  assumes "x  {a ..b}"
  assumes "- inverse (b*b)  alpha"
  shows inverse_linear_lower: "inverse b + alpha * (x - b)  inverse x" (is ?lower)
    and inverse_linear_upper: "inverse x  inverse a + alpha * (x - a)" (is ?upper)
proof -
  have deriv_inv:
    "x. x  {a .. b}  (inverse has_field_derivative - inverse (x*x)) (at x within {a .. b})"
    using assms
    by (auto intro!: derivative_eq_intros)
  show ?lower
    using assms
    by (intro linear_lower[OF deriv_inv])
        (auto simp: mult_mono intro!:  order_trans[OF _ assms(3)])
  show ?upper
    using assms
    by (intro linear_upper[OF deriv_inv])
        (auto simp: mult_mono intro!:  order_trans[OF _ assms(3)])
qed


subsubsection ‹Inverse›

definition inverse_aform'::"nat  real aform  real aform × real" where
  "inverse_aform' p X = (
    let l = Inf_aform' p X in
    let u = Sup_aform' p X in
    let a = min (abs l) (abs u) in
    let b = max (abs l) (abs u) in
    let sq = truncate_up p (b * b) in
    let alpha = - real_divl p 1 sq in
    let dmax = truncate_up p (real_divr p 1 a - alpha * a) in
    let dmin = truncate_down p (real_divl p 1 b - alpha * b) in
    let zeta' = truncate_up p ((dmin + dmax) / 2) in
    let zeta = if l < 0 then - zeta' else zeta' in
    let delta = truncate_up p (zeta - dmin) in
    let res1 = trunc_bound_eucl p (alpha * fst X) in
    let res2 = trunc_bound_eucl p (fst res1 + zeta) in
    let zs = trunc_bound_pdevs p (scaleR_pdevs alpha (snd X)) in
    ((fst res2, fst zs), (sum_list' p [delta, snd res1, snd res2, snd zs])))"

lemma inverse_aform'E:
  fixes X::"real aform"
  assumes e: "e  UNIV  {-1 .. 1}"
  assumes Inf_pos: "Inf_aform' p X > 0"
  assumes "x = aform_val e X"
  shows "inverse x  aform_err e (inverse_aform' p X)"
proof -
  define l where "l = Inf_aform' p X"
  define u where "u = Sup_aform' p X"
  define a where "a = min (abs l) (abs u)"
  define b where "b = max (abs l) (abs u)"
  define sq where "sq = truncate_up p (b * b)"
  define alpha where "alpha = - (real_divl p 1 sq)"
  define d_max' where "d_max' = truncate_up p (real_divr p 1 a - alpha * a)"
  define d_min' where "d_min' = truncate_down p (real_divl p 1 b - alpha * b)"
  define zeta where "zeta = truncate_up p ((d_min' + d_max') / 2)"
  define delta where "delta = truncate_up p (zeta - d_min')"
  note vars = l_def u_def a_def b_def sq_def alpha_def d_max'_def d_min'_def zeta_def delta_def
  let ?x = "aform_val e X"

  have "0 < l" using assms by (auto simp add: l_def Inf_aform_def)
  have "l  u" by (auto simp: l_def u_def)

  hence a_def': "a = l" and b_def': "b = u" and "0 < a" "0 < b"
    using 0 < l by (simp_all add: a_def b_def)
  have "0 < ?x"
    by (rule less_le_trans[OF Inf_pos order.trans[OF Inf_aform' Inf_aform], OF e])
  have "a  ?x"
    by (metis order.trans Inf_aform e Inf_aform' a_def' l_def)
  have "?x  b"
    by (metis order.trans Sup_aform e Sup_aform' b_def' u_def)
  hence "?x  {?x .. b}"
    by simp

  have "- inverse (b * b)  alpha"
    by (auto simp add: alpha_def inverse_mult_distrib[symmetric] inverse_eq_divide sq_def
      intro!: order_trans[OF real_divl] divide_left_mono truncate_up mult_pos_pos 0 < b)

  {
    note 0 < a
    moreover
    have "?x  {a .. b}" using a  ?x ?x  b by simp
    moreover
    note - inverse (b * b)  alpha
    ultimately have "inverse ?x  inverse a + alpha * (?x - a)"
      by (rule inverse_linear_upper)
    also have " = alpha * ?x + (inverse a - alpha * a)"
      by (simp add: algebra_simps)
    also have "inverse a - (alpha * a)  (real_divr p 1 a - alpha * a)"
      by (auto simp: inverse_eq_divide real_divr)
    also have "  (truncate_down p (real_divl p 1 b - alpha * b) +
          (real_divr p 1 a - alpha * a)) / 2 +
        (truncate_up p (real_divr p 1 a - alpha * a) -
          truncate_down p (real_divl p 1 b - alpha * b)) / 2"
      (is "_  (truncate_down p ?lb + ?ra) / 2 + (truncate_up p ?ra - truncate_down p ?lb) / 2")
      by (auto simp add: field_simps
        intro!: order_trans[OF _ add_left_mono[OF mult_left_mono[OF truncate_up]]])
    also have "(truncate_down p ?lb + ?ra) / 2 
        truncate_up p ((truncate_down p ?lb + truncate_up p ?ra) / 2)"
      by (intro truncate_up_le divide_right_mono add_left_mono truncate_up) auto
    also
    have "(truncate_up p ?ra - truncate_down p ?lb) / 2 + truncate_down p ?lb 
        (truncate_up p ((truncate_down p ?lb + truncate_up p ?ra) / 2))"
      by (rule truncate_up_le) (simp add: field_simps)
    hence "(truncate_up p ?ra - truncate_down p ?lb) / 2  truncate_up p
        (truncate_up p ((truncate_down p ?lb + truncate_up p ?ra) / 2) - truncate_down p ?lb)"
      by (intro truncate_up_le) (simp add: field_simps)
    finally have "inverse ?x  alpha * ?x + zeta + delta"
      by (auto simp: zeta_def delta_def d_min'_def d_max'_def right_diff_distrib ac_simps)
  } note upper = this

  {
    have "alpha * b + truncate_down p (real_divl p 1 b - alpha * b)  inverse b"
      by (rule order_trans[OF add_left_mono[OF truncate_down]])
        (auto simp: inverse_eq_divide real_divl)
    hence "zeta + alpha * b  delta + inverse b"
      by (auto simp: zeta_def delta_def d_min'_def d_max'_def right_diff_distrib
        intro!: order_trans[OF _ add_right_mono[OF truncate_up]])
    hence "alpha * ?x + zeta - delta  inverse b + alpha * (?x - b)"
      by (simp add: algebra_simps)
    also
    {
      note 0 < aform_val e X
      moreover
      note aform_val e X  {aform_val e X .. b}
      moreover

      note - inverse (b * b)  alpha
      ultimately
      have "inverse b + alpha * (aform_val e X - b)  inverse (aform_val e X)"
        by (rule inverse_linear_lower)
    }
    finally have "alpha * (aform_val e X) + zeta - delta  inverse (aform_val e X)" .
  } note lower = this


  have "inverse (aform_val e X) = alpha * (aform_val e X) + zeta +
      (inverse (aform_val e X) - alpha * (aform_val e X) - zeta)"
    (is "_ = _ + ?linerr")
    by simp
  also
  have "?linerr  {- delta .. delta}"
    using lower upper by simp
  hence linerr_le: "abs ?linerr  delta"
    by auto

  let ?z0 = "trunc_bound_eucl p (alpha * fst X)"
  from trunc_bound_euclE
  obtain e1 where abs_e1: "¦e1¦  snd ?z0" and e1: "fst ?z0 = alpha * fst X + e1"
    by blast
  let ?z1 = "trunc_bound_eucl p (fst ?z0 + zeta)"
  from trunc_bound_euclE
  obtain e1' where abs_e1': "¦e1'¦  snd ?z1" and e1': "fst ?z1 = fst ?z0 + zeta + e1'"
    by blast

  let ?zs = "trunc_bound_pdevs p (scaleR_pdevs alpha (snd X))"
  from trunc_bound_pdevsE[OF e]
  obtain e2 where abs_e2: "¦e2¦  snd (?zs)"
    and e2: "pdevs_val e (fst ?zs) = pdevs_val e (scaleR_pdevs alpha (snd X)) + e2"
    by blast

  have "alpha * (aform_val e X) + zeta =
      aform_val e (fst (inverse_aform' p X)) + (- e1 - e1' - e2)"
    unfolding inverse_aform'_def Let_def vars[symmetric]
    using 0 < l
    by (simp add: aform_val_def assms e1') (simp add: e1 e2 algebra_simps)
  also
  let ?err = "(- e1 - e1' - e2 + inverse (aform_val e X) - alpha * aform_val e X - zeta)"
  {
    have "abs ?err  abs ?linerr + abs e1 + abs e1' + abs e2"
      by simp
    also have "  delta + snd ?z0 + snd ?z1 + snd ?zs"
      by (blast intro: add_mono linerr_le abs_e1 abs_e1' abs_e2)
    also have "  (snd (inverse_aform' p X))"
      unfolding inverse_aform'_def Let_def vars[symmetric]
      using 0 < l
      by (auto simp add: inverse_aform'_def pdevs_apply_trunc_pdevs assms vars[symmetric]
        intro!: order.trans[OF _ sum_list'_sum_list_le])
    finally have "abs ?err  snd (inverse_aform' p X)" by simp
  } note err_le = this
  have "aform_val (e) (fst (inverse_aform' p X)) + (- e1 - e1' - e2) +
    (inverse (aform_val e X) - alpha * aform_val e X - zeta) =
    aform_val e (fst (inverse_aform' p X)) + ?err"
    by simp
  finally
  show ?thesis
    apply (intro aform_errI)
    using err_le
    by (auto simp: assms)
qed

definition "inverse_aform p a =
  do {
    let l = Inf_aform' p a;
    let u = Sup_aform' p a;
    if (l  0  0  u) then None
    else if (l  0) then (Some (apfst uminus_aform (inverse_aform' p (uminus_aform a))))
    else Some (inverse_aform' p a)
  }"

lemma eucl_truncate_up_eq_eucl_truncate_down:
  "eucl_truncate_up p x = - (eucl_truncate_down p (- x))"
  by (auto simp: eucl_truncate_up_def eucl_truncate_down_def truncate_up_eq_truncate_down sum_negf)

lemma inverse_aformE:
  fixes X::"real aform"
  assumes e: "e  UNIV  {-1 .. 1}"
    and disj: "Inf_aform' p X > 0  Sup_aform' p X < 0"
  obtains Y where
    "inverse_aform p X = Some Y"
    "inverse (aform_val e X)  aform_err e Y"
proof -
  {
    assume neg: "Sup_aform' p X < 0"
    from neg have [simp]: "Inf_aform' p X  0"
      by (metis Inf_aform'_le_Sup_aform' dual_order.strict_trans1 less_asym not_less)
    from neg disj have "0 < Inf_aform' p (uminus_aform X)"
      by (auto simp: Inf_aform'_def Sup_aform'_def eucl_truncate_up_eq_eucl_truncate_down ac_simps)
    from inverse_aform'E[OF e(1) this]
    have iin: "inverse (aform_val e (uminus_aform X))  aform_err e (inverse_aform' p (uminus_aform X))"
      by simp
    let ?Y = "apfst uminus_aform (inverse_aform' p (uminus_aform X))"
    have "inverse_aform p X = Some ?Y"
      "inverse (aform_val e X)  aform_err e ?Y"
      using neg iin by (auto simp: inverse_aform_def aform_err_def)
    then have ?thesis ..
  } moreover {
    assume pos: "Inf_aform' p X > 0"
    from pos have eq: "inverse_aform p X = Some (inverse_aform' p X)"
      by (auto simp: inverse_aform_def)
    moreover
    from inverse_aform'E[OF e(1) pos refl]
    have "inverse (aform_val e X)  aform_err e (inverse_aform' p X)" .
    ultimately have ?thesis ..
  } ultimately show ?thesis
    using assms by auto
qed

definition aform_err_to_aform::"aform_err  nat  real aform"
  where "aform_err_to_aform X n = (fst (fst X),  pdev_upd (snd (fst X)) n (snd X))"

lemma aform_err_to_aformE:
  assumes "x  aform_err e X"
  assumes deg: "degree_aform_err X  n"
  obtains err where "x = aform_val (e(n:=err)) (aform_err_to_aform X n)"
    "-1  err" "err  1"
proof -
  from aform_errE[OF assms(1)] have "¦x - aform_val e (fst X)¦  snd X" by auto
  from error_absE[OF this] obtain err where err:
    "x - aform_val e (fst X) = err * snd X" "err  {- 1..1}"
    by auto
  have "x = aform_val (e(n:=err)) (aform_err_to_aform X n)"
    "-1  err" "err  1"
    using err deg
    by (auto simp: aform_val_def aform_err_to_aform_def)
  then show ?thesis ..
qed

definition aform_to_aform_err::"real aform  nat  aform_err"
  where "aform_to_aform_err X n = ((fst X,  pdev_upd (snd X) n 0), abs (pdevs_apply (snd X) n))"

lemma aform_to_aform_err: "aform_val e X  aform_err e (aform_to_aform_err X n)"
  if "e  UNIV  {-1 .. 1}"
proof -
  from that have abs_e[simp]: "i. ¦e i¦  1" by (auto simp: abs_real_def)
  have "- e n * pdevs_apply (snd X) n  ¦pdevs_apply (snd X) n¦"
  proof -
    have "- e n * pdevs_apply (snd X) n  ¦- e n * pdevs_apply (snd X) n¦"
      by auto
    also have "  abs (pdevs_apply (snd X) n)"
      using that
      by (auto simp: abs_mult intro!: mult_left_le_one_le)
    finally show ?thesis .
  qed
  moreover have "e n * pdevs_apply (snd X) n  ¦pdevs_apply (snd X) n¦"
  proof -
    have "e n * pdevs_apply (snd X) n  ¦e n * pdevs_apply (snd X) n¦"
      by auto
    also have "  abs (pdevs_apply (snd X) n)"
      using that
      by (auto simp: abs_mult intro!: mult_left_le_one_le)
    finally show ?thesis .
  qed
  ultimately
  show ?thesis
    by (auto simp: aform_to_aform_err_def aform_err_def aform_val_def)
qed

definition "acc_err p x e  (fst x, truncate_up p (snd x + e))"

definition ivl_err :: "real interval  (real × real pdevs) × real"
  where "ivl_err ivl  (((upper ivl + lower ivl)/