# 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="∑i∈Basis. (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¦) [xa←rev [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 "[xa←rev [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 ⟷ (∀x∈set 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 "∀x∈set 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)/