Theory Refine_Rigorous_Numerics_Aform

theory Refine_Rigorous_Numerics_Aform
  imports
    Refine_Rigorous_Numerics
    "HOL-Types_To_Sets.Types_To_Sets"
begin

lemma Joints_ne_empty[simp]: "Joints xs ≠ {}" "{} ≠ Joints xs"
  by (auto simp: Joints_def valuate_def)

lemma Inf_aform_le_Affine: "x ∈ Affine X ⟹ Inf_aform X ≤ x"
  by (auto simp: Affine_def valuate_def intro!: Inf_aform)
lemma Sup_aform_ge_Affine: "x ∈ Affine X ⟹ x ≤ Sup_aform X"
  by (auto simp: Affine_def valuate_def intro!: Sup_aform)

lemmas Inf_aform'_Affine_le = order_trans[OF Inf_aform' Inf_aform_le_Affine]
lemmas Sup_aform'_Affine_ge = order_trans[OF Sup_aform_ge_Affine Sup_aform']

fun approx_form_aform :: "nat ⇒ form ⇒ real aform list ⇒ bool" where
"approx_form_aform prec (Less a b) bs =
   (case (approx_floatariths prec [a - b] bs)
   of Some [r] ⇒ Sup_aform' prec r < 0
    | _ ⇒ False)"
| "approx_form_aform prec (LessEqual a b) bs =
   (case (approx_floatariths prec [a - b] bs)
   of Some [r] ⇒ Sup_aform' prec r ≤ 0
    | _ ⇒ False)"
| "approx_form_aform prec (AtLeastAtMost a b c) bs =
   (case (approx_floatariths prec [a - b, a - c] bs)
   of Some [r, s] ⇒ 0 ≤ Inf_aform' prec r ∧  Sup_aform' prec s ≤ 0
    | _ ⇒ False)"
| "approx_form_aform prec (Bound a b c d) bs = False"
| "approx_form_aform prec (Assign a b c) bs = False"
| "approx_form_aform prec (Conj a b) bs ⟷
    approx_form_aform prec a bs ∧ approx_form_aform prec b bs"
| "approx_form_aform prec (Disj a b) bs ⟷
    approx_form_aform prec a bs ∨ approx_form_aform prec b bs"

lemma in_Joints_intervalD:
  "x ∈ {Inf_aform' p X .. Sup_aform' p X} ∧ xs ∈ Joints XS"
  if "x#xs ∈ Joints (X#XS)"
  using that
  by (auto simp: Joints_def valuate_def Affine_def
      intro!: Inf_aform'_Affine_le Sup_aform'_Affine_ge)

lemma approx_form_aform:
  "interpret_form f vs"
  if "approx_form_aform p f VS" "vs ∈ Joints VS"
  using that
  by (induction f)
    (auto split: option.splits list.splits simp: algebra_simps
      dest!: approx_floatariths_outer
      dest!: in_Joints_intervalD[where p=p])

fun msum_aforms::"nat ⇒ real aform list ⇒ real aform list ⇒ real aform list"  where
  "msum_aforms d (x#xs) (y#ys) = msum_aform d x y#msum_aforms d xs ys"
| "msum_aforms d _ _ = []"

definition [simp]: "degree_aforms_real = (degree_aforms::real aform list ⇒ nat)"

abbreviation "msum_aforms' ≡ λX Y. msum_aforms (degree_aforms_real (X @ Y)) X Y"

lemma aform_val_msum_aforms:
  assumes "degree_aforms xs ≤ d"
  shows "aform_vals e (msum_aforms d xs ys) = List.map2 (+) (aform_vals e xs) (aform_vals (λi. e (i + d)) ys)"
  using assms
proof (induction xs ys rule: msum_aforms.induct)
  case (1 d x xs y ys)
  from 1 have "degree_aforms xs ≤ d"
    by (auto simp: degrees_def)
  from 1(1)[OF this] 1
  have "aform_vals e (msum_aforms d xs ys) = List.map2 (+) (aform_vals e xs) (aform_vals (λi. e (i + d)) ys)"
    by simp
  then show ?case
    using 1
    by (simp add: aform_vals_def aform_val_msum_aform degrees_def)
qed (auto simp: aform_vals_def)

lemma Joints_msum_aforms:
  assumes "degree_aforms xs ≤ d"
  assumes "degree_aforms ys ≤ d"
  shows "Joints (msum_aforms d xs ys) = {List.map2 (+) a b |a b. a ∈ Joints xs ∧ b ∈ Joints ys}"
  apply (auto simp: Joints_def valuate_def aform_vals_def[symmetric]
      aform_val_msum_aforms assms)
   apply force
  subgoal for e e'
    apply (rule image_eqI[where x="λi. if i < d then e i else e' (i - d)"])
     apply (auto simp: Pi_iff)
  proof goal_cases
    case 1
    have "(aform_vals e xs) = aform_vals (λi. if i < d then e i else e' (i - d)) xs"
      using assms
      by (auto intro!: simp: aform_vals_def aform_val_def degrees_def intro!: pdevs_val_degree_cong)
    then show ?case
      by simp
  qed
  done

definition "split_aforms_largest_uncond_take n X =
  (let (i, x) = max_pdev (abssum_of_pdevs_list (map snd (take n X))) in split_aforms X i)"

text ‹intersection with plane›

definition
  "project_coord x b n = (∑i←Basis_list. (if i = b then n else if i = -b then -n else x ∙ i) *R i)"

lemma inner_eq_abs_times_sgn:
  "u ∙ b = abs u ∙ b * sgn (u ∙ b)" if "b ∈ Basis" for b::"'a::ordered_euclidean_space"
  by (subst sgn_mult_abs[symmetric]) (auto simp: that abs_inner )

lemma inner_Basis_eq_zero_absI: "x ∈ Basis ⟹ abs u ∈ Basis ⟹ x ≠ ¦u¦ ⟹ x ≠ 0 ⟹ u ∙ x = 0"
  for x::"'a::ordered_euclidean_space"
  apply (subst euclidean_inner)
  apply (auto simp: inner_Basis if_distribR if_distrib cong: if_cong)
  apply (subst inner_eq_abs_times_sgn)
  by (auto simp: inner_Basis)

lemma abs_in_BasisE:
  fixes u::"'a::ordered_euclidean_space"
  assumes "abs u ∈ Basis"
  obtains "abs u = u" | "abs u = - u"
proof (cases "u ∙ abs u = 1")
  case True
  have "abs u = (∑i∈Basis. (if i = abs u then (abs u ∙ i) *R i else 0))"
    using assms
    by auto
  also have "… = (∑i∈Basis. (if i = abs u then (u ∙ i) *R i else 0))"
    by (rule sum.cong) (auto simp: True)
  also have "… = (∑i∈Basis. (u ∙ i) *R i)"
    by (rule sum.cong) (auto simp: inner_Basis_eq_zero_absI assms)
  also have "… = u" by (simp add: euclidean_representation)
  finally show ?thesis ..
next
  case False
  then have F: "u ∙ ¦u¦ = -1"
    apply (subst inner_eq_abs_times_sgn[OF assms]) 
    apply (subst (asm) inner_eq_abs_times_sgn[OF assms]) 
    using assms
    apply (auto simp: inner_Basis sgn_if)
    by (metis (full_types) abs_0_eq euclidean_all_zero_iff inner_Basis_eq_zero_absI nonzero_Basis)
  have "abs u = (∑i∈Basis. (if i = abs u then (abs u ∙ i) *R i else 0))"
    using assms
    by auto
  also have "… = (∑i∈Basis. (if i = abs u then (- u ∙ i) *R i else 0))"
    by (rule sum.cong) (auto simp: F)
  also have "… = (∑i∈Basis. (- u ∙ i) *R i)"
    by (rule sum.cong) (auto simp: inner_Basis_eq_zero_absI assms)
  also have "… = - u" by (subst euclidean_representation) simp
  finally show ?thesis ..
qed
  
lemma abs_in_Basis_iff:
  fixes u::"'a::ordered_euclidean_space"
  shows "abs u ∈ Basis ⟷ u ∈ Basis ∨ - u ∈ Basis"
proof -
  have u: "u = (∑i∈Basis. (u ∙ i) *R i)"
    by (simp add: euclidean_representation)
  have u': "- u = (∑i∈Basis. (- (u ∙ i)) *R i)"
    by (subst u) (simp add: sum_negf)
  have au: "abs u = (∑i∈Basis. ¦u ∙ i¦ *R i)"
    by (simp add: eucl_abs[where 'a='a])
  have "(u ∈ Basis ∨ - u ∈ Basis)" if "(¦u¦ ∈ Basis)"
    apply (rule abs_in_BasisE[OF that])
    using that
    by auto
  show ?thesis
    apply safe
    subgoal premises prems
      using prems(1)
      apply (rule abs_in_BasisE)
      using prems
      by simp_all
    subgoal
      apply (subst au)
      apply (subst (asm) u)
      apply (subst sum.cong[OF refl])
       apply (subst abs_inner[symmetric]) apply auto
      using u apply auto[1]
      done
    subgoal
      apply (subst au)
      apply (subst (asm) u')
      apply (subst sum.cong[OF refl])
      apply (subst abs_inner[symmetric]) apply auto
      using u' apply auto
      by (metis Basis_nonneg abs_of_nonpos inner_minus_left neg_0_le_iff_le scaleR_left.minus)
    done
qed

lemma abs_inner_Basis:
  fixes u v::"'a::ordered_euclidean_space"
  assumes "abs u ∈ Basis" "abs v ∈ Basis"
  shows "inner u v = (if u = v then 1 else if u = -v then -1 else 0)"
proof -
  define i where "i ≡ if u ∈ Basis then u else -u"
  define j where "j ≡ if v ∈ Basis then v else -v"
  have u: "u = (if u ∈ Basis then i else - i)"
    and v: "v = (if v ∈ Basis then j else - j)"
    by (auto simp: i_def j_def)
  have "i ∈ Basis" "j ∈ Basis"
    using assms
     by (auto simp: i_def j_def abs_in_Basis_iff)
   then show ?thesis
     apply (subst u)
     apply (subst v)
     apply (auto simp: inner_Basis)
     apply (auto simp: j_def i_def split: if_splits)
     using Basis_nonneg abs_of_nonpos by fastforce
 qed

lemma
  project_coord_inner_Basis:
  assumes "i ∈ Basis"
  shows "(project_coord x b n) ∙ i = (if i = b then n else if i = -b then -n else x ∙ i)"
proof -
  have "project_coord x b n ∙ i =
    (∑j∈Basis. (if j = b then n else if j = -b then -n else x ∙ j) * (if j = i then 1 else 0))"
    using assms
    by (auto simp: project_coord_def inner_sum_left inner_Basis)
  also have "… = (∑j∈Basis. (if j = i then (if j = b then n else if j = -b then -n else x ∙ j) else 0))"
    by (rule sum.cong) auto
  also have "… = (if i = b then n else if i = -b then -n else x ∙ i)"
    using assms
    by (subst sum.delta) auto
  finally show ?thesis by simp
qed

lemma
  project_coord_inner:
  assumes "abs i ∈ Basis"
  shows "(project_coord x b n) ∙ i = (if i = b then n else if i = -b then -n else x ∙ i)"
proof -
  consider "i ∈ Basis" | "- i ∈ Basis"
    using assms
    by (auto simp: abs_in_Basis_iff)
  then show ?thesis
  proof cases
    case 1
    then show ?thesis by (simp add: project_coord_inner_Basis)
  next
    case 2
    have "project_coord x b n ∙ i = - (project_coord x b n ∙ - i)"
      by simp
    also have "… = - (if - i = b then n else if - i = -b then -n else x ∙ - i)"
      using 2
      by (subst project_coord_inner_Basis) simp_all
    also have "… = (if i = b then n else if i = -b then -n else x ∙ i)"
      using 2 by auto (metis Basis_nonneg abs_le_zero_iff abs_of_nonneg neg_le_0_iff_le nonzero_Basis)
    finally show ?thesis .
  qed
qed

lift_definition project_coord_pdevs::"'a::executable_euclidean_space sctn ⇒ 'a pdevs ⇒ 'a pdevs" is
  "λsctn xs i. project_coord (xs i) (normal sctn) 0"
  by (rule finite_subset) (auto simp: project_coord_def cong: if_cong)

definition "project_coord_aform sctn cxs =
  (project_coord (fst cxs) (normal sctn) (pstn sctn), project_coord_pdevs sctn (snd cxs))"

definition project_coord_aform_lv :: "real aform list ⇒ nat ⇒ real ⇒ real aform list"
  where "project_coord_aform_lv xs b n = xs[b:=(n, zero_pdevs)]"

definition "project_ncoord_aform x b n = project_coord_aform (Sctn (Basis_list ! b) n) x"

lemma
  finite_sum_subset_aux_lemma:
  assumes "finite s"
  shows " {i. (∑x∈s. f x i) ≠ 0} ⊆ {i. ∃x∈s. f x i ≠ 0}"
proof -
  have "{i. (∑x∈s. f x i) ≠ 0} = - {i. (∑x∈s. f x i) = 0}"
    by auto
  also have "… ⊆ - {i. ∀x ∈ s. f x i = 0}"
    by auto
  also have "… = {i. ∃x ∈ s. f x i ≠ 0}"
    by auto
  finally show ?thesis by simp
qed

lift_definition sum_pdevs::"('i ⇒ 'a::comm_monoid_add pdevs) ⇒ 'i set ⇒  'a pdevs"
  is "λf X. if finite X then (λi. ∑x∈X. f x i) else (λ_. 0)"
  apply auto
  apply (rule finite_subset)
   apply (rule finite_sum_subset_aux_lemma)
   apply auto
  done

lemma pdevs_apply_sum_pdevs[simp]:
  "pdevs_apply (sum_pdevs f X) i = (∑x∈X. pdevs_apply (f x) i)"
  by transfer auto

lemma sum_pdevs_empty[simp]: "sum_pdevs f {} = zero_pdevs"
  by transfer auto

lemma sum_pdevs_insert[simp]:
  "finite xs ⟹ sum_pdevs f (insert a xs) =
    (if a ∈ xs then sum_pdevs f xs else add_pdevs (f a) (sum_pdevs f xs))"
  by (auto simp: insert_absorb intro!: pdevs_eqI)

lemma sum_pdevs_infinite[simp]: "infinite X ⟹ sum_pdevs f X = zero_pdevs"
  by transfer auto

lemma compute_sum_pdevs[code]:
  "sum_pdevs f (set XS) = foldr (λa b. add_pdevs (f a) b) (remdups XS) zero_pdevs"
  by (induction XS) auto

lemma degree_sum_pdevs_le:
  "degree (sum_pdevs f X) ≤ Max (degree ` f ` X)"
  apply (rule degree_le)
  apply auto
  apply (cases "X = {}")
  subgoal by simp
  subgoal by (cases "finite X") simp_all
  done

lemma pdevs_val_sum_pdevs:
  "pdevs_val e (sum_pdevs f X) = (∑x∈X. pdevs_val e (f x))"
  apply (cases "finite X")
  subgoal
    apply (subst pdevs_val_sum_le)
     apply (rule degree_sum_pdevs_le)
    apply (auto simp: scaleR_sum_right)
    apply (subst sum.swap)
    apply (rule sum.cong[OF refl])
    apply (subst pdevs_val_sum_le[where d="Max (degree ` f ` X)"])
     apply (auto simp: Max_ge_iff)
    done
  subgoal by simp
  done

definition eucl_of_list_aform :: "(real × real pdevs) list ⇒ 'a::executable_euclidean_space aform"
where "eucl_of_list_aform xs =
  (eucl_of_list (map fst xs), sum_pdevs (λi. pdevs_scaleR (snd (xs ! index Basis_list i)) i) Basis)"

definition lv_aforms_rel::"(((real × real pdevs) list) × ('a::executable_euclidean_space aform)) set"
  where "lv_aforms_rel = br eucl_of_list_aform (λxs. length xs = DIM('a))"


definition "inner_aforms' p X Y =
  (let fas = [inner_floatariths (map floatarith.Var [0..<length X])
    (map floatarith.Var[length X..<length X + length Y])]
  in
    approx_slp_outer p (length fas) fas (X@Y)
    )"

lemma
  affine_extension_AffineD:
  assumes "affine_extension2 (λd x y. Some (F d x y)) f"
  assumes "[x, y] ∈ Joints [X, Y]"
  assumes "d ≥ degree_aform X"
  assumes "d ≥ degree_aform Y"
  shows "f x y ∈ Affine (F d X Y)"
proof -
  from assms(2) obtain e where e:
    "e ∈ funcset UNIV {-1 .. 1}"
    "x = aform_val e X"
    "y = aform_val e Y"
    by (auto simp: Joints_def valuate_def)
  from affine_extension2E[OF assms(1) refl e(1) assms(3) assms(4)]
  obtain e' where "e' ∈ funcset UNIV {- 1..1}" "f (aform_val e X) (aform_val e Y) = aform_val e' (F d X Y)"
    "⋀n. n < d ⟹ e' n = e n" "aform_val e X = aform_val e' X" "aform_val e Y = aform_val e' Y"
    by auto
  then show ?thesis by (auto simp: Affine_def valuate_def e)
qed

lemma aform_val_zero_pdevs[simp]: "aform_val e (x, zero_pdevs) = x"
  by (auto simp: aform_val_def)

lemma neg_equal_zero_eucl[simp]: "-a = a ⟷ a = 0" for a::"'a::euclidean_space"
  by (auto simp: algebra_simps euclidean_eq_iff[where 'a='a])

context includes autoref_syntax begin

lemma Option_bind_param[param, autoref_rules]:
  "((⤜), (⤜)) ∈ ⟨S⟩option_rel → (S → ⟨R⟩option_rel) → ⟨R⟩option_rel"
  unfolding Option.bind_def
  by parametricity

lemma zip_Basis_list_pat[autoref_op_pat_def]: "∑(b, m)←zip Basis_list ms. m *R b ≡ OP eucl_of_list $ ms"
proof (rule eq_reflection)
  have z: "zip ms (Basis_list::'a list) = map (λ(x, y). (y, x)) (zip Basis_list ms)"
    by (subst zip_commute) simp
  show "(∑(b, m)←zip (Basis_list::'a list) ms. m *R b) = OP eucl_of_list $ ms"
    unfolding eucl_of_list_def autoref_tag_defs
    apply (subst z)
    apply (subst map_map)
    apply (auto cong: map_cong simp: split_beta')
    done
qed

lemma length_aforms_of_ivls: "length (aforms_of_ivls a b) = min (length a) (length b)"
  by (auto simp: aforms_of_ivls_def)

lemma length_lv_rel:
  "(ys, y::'a) ∈ lv_rel ⟹ length ys = DIM('a::executable_euclidean_space)"
  by (auto simp: lv_rel_def br_def)

lemma lv_rel_nth[simp]: "(xs, x::'a::executable_euclidean_space) ∈ lv_rel ⟹ b ∈ Basis ⟹
  xs ! index (Basis_list) b = x ∙ b"
  unfolding lv_rel_def
  by (auto simp: br_def eucl_of_list_inner)

lemma aform_of_ivl_autoref[autoref_rules]:
  "(aforms_of_ivls, aform_of_ivl) ∈ lv_rel → lv_rel → lv_aforms_rel"
  apply (auto simp: lv_aforms_rel_def br_def eucl_of_list_aform_def length_aforms_of_ivls length_lv_rel)
  subgoal for xs x ys y
    apply (auto simp: aform_of_ivl_def aforms_of_ivls_def o_def eucl_of_list_inner inner_simps
        pdevs_apply_pdevs_of_ivl length_lv_rel
        intro!: euclidean_eqI[where 'a='a] pdevs_eqI)
    by (metis index_Basis_list_nth inner_not_same_Basis length_Basis_list nth_Basis_list_in_Basis)
  done

lemma bound_intersect_2d_ud[autoref_rules]:
  shows "(bound_intersect_2d_ud, bound_intersect_2d_ud) ∈ nat_rel → ((rnv_rel ×r rnv_rel) ×r Id) → rnv_rel → ⟨rnv_rel ×r rnv_rel⟩option_rel"
  by auto

lemma eucl_of_list_autoref[autoref_rules]:
  includes autoref_syntax
  assumes "SIDE_PRECOND (length xs = DIM('a::executable_euclidean_space))"
  assumes "(xsi, xs) ∈ ⟨rnv_rel⟩list_rel"
  shows "(xsi, eucl_of_list $ xs::'a) ∈ lv_rel"
  using assms
  unfolding lv_rel_def
  by (auto simp: br_def)

definition "inner2s x b c = (inner_lv_rel x b, inner_lv_rel x c)"

lemma inner_lv_rel_autoref[autoref_rules]:
  "(inner_lv_rel, (∙)) ∈ lv_rel → lv_rel → rnv_rel"
  using lv_rel_inner[unfolded inner_lv_rel_def[symmetric]]
  by auto

lemma inner_lv_rel_eq:
  "⟦length xs = DIM('a::executable_euclidean_space); (xa, x'a) ∈ lv_rel⟧ ⟹
  inner_lv_rel xs xa = eucl_of_list xs ∙ (x'a::'a)"
  using inner_lv_rel_autoref[param_fo, of "xs" "eucl_of_list xs" xa x'a]
  unfolding lv_rel_def
  by (auto simp: br_def)

definition "inner_pdevs xs ys = sum_pdevs (λi. scaleR_pdevs (ys ! i) (xs ! i)) {..<length xs}"

definition "pdevs_inner2s xs a b = prod_of_pdevs (inner_pdevs xs a) (inner_pdevs xs b)"

lemma inner2_aform_autoref[autoref_rules]:
  shows "((λxs a b. (inner2s (map fst xs) a b, pdevs_inner2s (map snd xs) a b)), inner2_aform) ∈ lv_aforms_rel → lv_rel → lv_rel → ((rnv_rel ×r rnv_rel)×rId)"
  unfolding inner2_aform_def
  apply (auto simp: lv_aforms_rel_def br_def eucl_of_list_aform_def inner2_aform_def )
   apply (auto simp: inner2s_def inner2_def inner_lv_rel_eq pdevs_inner2s_def inner_pdevs_def
      sum_Basis_sum_nth_Basis_list inner_sum_left inner_Basis mult.commute
      intro!: pdevs_eqI)
  subgoal for a b c d e f
    apply (rule sum.cong)
     apply force
    subgoal for n by (auto dest!: lv_rel_nth[where b="Basis_list ! n"] simp: inner_commute)
    done
  subgoal for a b c d e f
    apply (rule sum.cong)
     apply force
    subgoal for n by (auto dest!: lv_rel_nth[where b="Basis_list ! n"] simp: inner_commute)
    done
  done

lemma RETURN_inter_aform_plane_ortho_annot:
  "RETURN (inter_aform_plane_ortho p (Z::'a aform) n g) =
    (case those (map (λb. bound_intersect_2d_ud p (inner2_aform Z n b) g) Basis_list) of
      Some mMs ⇒
      do {
        ASSERT (length mMs = DIM('a::executable_euclidean_space));
        let l = (∑(b, m)←zip Basis_list (map fst mMs). m *R b);
        let u = (∑(b, M)←zip Basis_list (map snd mMs). M *R b);
        RETURN (Some (aform_of_ivl l u))
      }
  | None ⇒ RETURN None)"
  apply (auto simp: inter_aform_plane_ortho_def split: option.splits)
  subgoal premises prems for x2
  proof -
    have "length x2 = DIM('a)" using prems
      using map_eq_imp_length_eq by (force simp: those_eq_Some_map_Some_iff)
    then show ?thesis using prems by auto
  qed
  done

definition "inter_aform_plane_ortho_nres p Z n g = RETURN (inter_aform_plane_ortho p Z n g)"

schematic_goal inter_aform_plane_ortho_lv:
  fixes Z::"'a::executable_euclidean_space aform"
  assumes [autoref_rules_raw]: "DIM_precond TYPE('a) D"
  assumes [autoref_rules]: "(pp, p) ∈ nat_rel" "(Zi, Z) ∈ lv_aforms_rel"
    "(ni, n) ∈ lv_rel" "(gi, g) ∈ rnv_rel"
  notes [autoref_rules] = those_param param_map
  shows "(nres_of (?f::?'a dres), inter_aform_plane_ortho_nres $ p $ Z $ n $ g) ∈ ?R"
  unfolding autoref_tag_defs
  unfolding RETURN_inter_aform_plane_ortho_annot inter_aform_plane_ortho_nres_def
  including art
  by autoref_monadic
concrete_definition inter_aform_plane_ortho_lv for pp Zi ni gi uses inter_aform_plane_ortho_lv
lemmas [autoref_rules] = inter_aform_plane_ortho_lv.refine

lemma project_coord_lv[autoref_rules]:
  assumes "(xs, x::'a::executable_euclidean_space aform) ∈ lv_aforms_rel" "(ni, n) ∈ nat_rel"
  assumes "SIDE_PRECOND (n < DIM('a))"
  shows "(project_coord_aform_lv xs ni, project_ncoord_aform $ x $ n) ∈ rnv_rel → lv_aforms_rel"
  using assms
  apply (auto simp: project_coord_aform_lv_def project_ncoord_aform_def lv_rel_def
      project_coord_aform_def eucl_of_list_aform_def
      lv_aforms_rel_def br_def)
  subgoal
    apply (auto intro!: euclidean_eqI[where 'a='a])
    apply (subst project_coord_inner_Basis)
     apply (auto simp: eucl_of_list_inner )
     apply (subst nth_list_update)
      apply auto
    using in_Basis_index_Basis_list apply force
    using inner_not_same_Basis nth_Basis_list_in_Basis apply fastforce
    apply (auto simp: nth_list_update)
    done
  subgoal for i
    apply (auto intro!: pdevs_eqI simp: project_coord_pdevs.rep_eq)
    apply (auto intro!: euclidean_eqI[where 'a='a])
    apply (subst project_coord_inner_Basis)
     apply (auto simp: eucl_of_list_inner )
     apply (subst nth_list_update)
      apply auto
    using inner_not_same_Basis nth_Basis_list_in_Basis apply fastforce
    apply (auto simp: nth_list_update)
    done
  done

definition inter_aform_plane
  where "inter_aform_plane prec Xs (sctn::'a sctn) =
    do {
      cxs ← inter_aform_plane_ortho_nres (prec) (Xs) (normal sctn) (pstn sctn);
      case cxs of
        Some cxs ⇒
          (if normal sctn ∈ set Basis_list
            then do {
              ASSERT (index Basis_list (normal sctn) < DIM('a::executable_euclidean_space));
              RETURN ((project_ncoord_aform cxs (index Basis_list (normal sctn)) (pstn sctn)))
            }
          else if - normal sctn ∈ set Basis_list
            then do {
              ASSERT (index Basis_list (-normal sctn) < DIM('a));
              RETURN ((project_ncoord_aform cxs (index Basis_list (-normal sctn)) (- pstn sctn)))
            }
          else SUCCEED)
      | None ⇒ SUCCEED
    }"

lemma [autoref_rules]:
  assumes [THEN GEN_OP_D, param]: "GEN_OP (=) (=) (A → A → bool_rel)"
  shows "(index, index) ∈ ⟨A⟩list_rel → A → nat_rel"
  unfolding index_def find_index_def
  by parametricity

schematic_goal inter_aform_plane_lv:
  fixes Z::"'a::executable_euclidean_space aform"
  assumes [autoref_rules_raw]: "DIM_precond TYPE('a) D"
  assumes [autoref_rules]: "(preci, prec) ∈ nat_rel" "(Zi, Z) ∈ lv_aforms_rel"
      "(sctni, sctn) ∈ ⟨lv_rel⟩sctn_rel"
  notes [autoref_rules] = those_param param_map
  shows "(nres_of (?f::?'a dres), inter_aform_plane prec Z sctn) ∈ ?R"
  unfolding autoref_tag_defs
  unfolding inter_aform_plane_def
  including art
  by (autoref_monadic )

concrete_definition inter_aform_plane_lv for preci Zi sctni uses inter_aform_plane_lv
lemmas [autoref_rules] = inter_aform_plane_lv.refine[autoref_higher_order_rule(1)]

end


lemma Basis_not_uminus:
  fixes i::"'a::euclidean_space"
  shows "i ∈ Basis ⟹ - i ∉ Basis"
  by (metis inner_Basis inner_minus_left one_neq_neg_one zero_neq_neg_one)

lemma
  pdevs_apply_project_coord_pdevs:
  assumes "normal sctn ∈ Basis ∨ - normal sctn ∈ Basis"
  assumes "i ∈ Basis"
  shows "pdevs_apply (project_coord_pdevs sctn cxs) x ∙ i =
    (if i = abs (normal sctn) then 0 else pdevs_apply cxs x ∙ i)"
  using assms[unfolded abs_in_Basis_iff[symmetric]]
  apply (transfer fixing: sctn)
  apply (auto simp: project_coord_inner Basis_not_uminus)
  using Basis_nonneg apply force
  using Basis_nonneg assms(1) by force

lemma pdevs_domain_project_coord_pdevs_subset:
  "pdevs_domain (project_coord_pdevs sctn X) ⊆ pdevs_domain X"
  apply auto
  apply (auto simp: euclidean_eq_iff[where 'a='a] )
  by (metis add.inverse_neutral project_coord_inner_Basis project_coord_pdevs.rep_eq)

lemma pdevs_val_project_coord_pdevs_inner_Basis:
  assumes "b ∈ Basis"
  shows "pdevs_val e (project_coord_pdevs sctn X) ∙ b =
    (if b = abs (normal sctn) then 0 else pdevs_val e X ∙ b)"
  using assms
  apply auto
   apply (simp add: pdevs_val_pdevs_domain inner_sum_left )
   apply (subst sum.cong[OF refl])
    apply (subst pdevs_apply_project_coord_pdevs)
      apply (simp add: abs_in_Basis_iff)
     apply simp
    apply (rule refl)
   apply auto
  unfolding pdevs_val_pdevs_domain inner_sum_left
  apply (rule sum.mono_neutral_cong_left)
     apply simp
    apply (rule pdevs_domain_project_coord_pdevs_subset)
     apply auto
   apply (metis Basis_nonneg abs_minus_cancel abs_of_nonneg euclidean_all_zero_iff
      project_coord_inner_Basis project_coord_pdevs.rep_eq)
   apply (metis Basis_nonneg abs_minus_cancel abs_of_nonneg
      project_coord_inner_Basis project_coord_pdevs.rep_eq)
  done

lemma inner_in_Sum: "b ∈ Basis ⟹ x ∙ b = (∑i∈Basis. x ∙ i * (i ∙ b))"
  apply (subst euclidean_representation[of x, symmetric])
  unfolding inner_sum_left
  by (auto simp: intro!: )

lemma aform_val_project_coord_aform:
  "aform_val e (project_coord_aform sctn X) = project_coord (aform_val e X) (normal sctn) (pstn sctn)"
  apply (auto simp: aform_val_def project_coord_def project_coord_aform_def)
  apply (rule euclidean_eqI)
  unfolding inner_add_left inner_sum_left
  apply (subst pdevs_val_project_coord_pdevs_inner_Basis)
   apply auto
   apply (rule sum.cong)
  apply auto
   apply (metis abs_in_Basis_iff abs_inner abs_inner_Basis abs_zero inner_commute)
  apply (subst inner_in_Sum[where x="pdevs_val e (snd X)"], force)
  unfolding sum.distrib[symmetric]
  apply (rule sum.cong)
  apply auto
  apply (metis Basis_nonneg abs_inner_Basis abs_of_nonneg abs_of_nonpos inner_commute neg_0_le_iff_le)
   apply (metis Basis_nonneg abs_inner_Basis abs_of_nonneg abs_of_nonpos neg_0_le_iff_le)
  apply (auto simp: inner_Basis)
  done

lemma project_coord_on_plane:
  assumes "x ∈ plane_of (Sctn n c)"
  shows "project_coord x n c = x"
  using assms
  by (auto simp: project_coord_def plane_of_def intro!: euclidean_eqI[where 'a='a])

lemma
  mem_Affine_project_coord_aformI:
  assumes "x ∈ Affine X"
  assumes "x ∈ plane_of sctn"
  shows "x ∈ Affine (project_coord_aform sctn X)"
proof -
  have "x = project_coord x (normal sctn) (pstn sctn)"
    using assms by (intro project_coord_on_plane[symmetric]) auto
  also
  from assms obtain e where "e ∈ funcset UNIV {-1 .. 1}" "x = aform_val e X"
    by (auto simp: Affine_def valuate_def)
  note this(2)
  also have "project_coord (aform_val e X) (normal sctn) (pstn sctn) =
    aform_val e (project_coord_aform sctn X)"
    by (simp add: aform_val_project_coord_aform)
  finally have "x = aform_val e (project_coord_aform sctn X)" unfolding ‹x = aform_val e X› .
  then show ?thesis using ‹e ∈ _›
    by (auto simp: Affine_def valuate_def)
qed

lemma
  mem_Affine_project_coord_aformD:
  assumes "x ∈ Affine (project_coord_aform sctn X)"
  assumes "abs (normal sctn) ∈ Basis"
  shows "x ∈ plane_of sctn"
  using assms
  by (auto simp: Affine_def valuate_def aform_val_project_coord_aform plane_of_def
      project_coord_inner)

definition "reduce_aform prec t X =
    summarize_aforms (prec) (collect_threshold (prec) 0 t) (degree_aforms X) X"

definition "correct_girard p m N (X::real aform list) =
  (let
    Xs = map snd X;
    M = pdevs_mapping Xs;
    D = domain_pdevs Xs;
    diff = m - card D
  in if 0 < diff then (λ_ _. True) else
    let
      Ds = sorted_list_of_set D;
      ortho_indices = map fst (take (diff + N) (sort_key (λ(i, r). r) (map (λi. let xs = M i in (i, sum_list' p (map abs xs) - fold max (map abs xs) 0)) Ds)));
      _ = ()
    in (λi (xs::real list). i ∉ set ortho_indices))"

definition "reduce_aforms prec C X = summarize_aforms (prec) C (degree_aforms X) X"

definition "pdevs_of_real x = (x, zero_pdevs)"

definition aform_inf_inner where "aform_inf_inner prec X n =
  (case inner_aforms' (prec) X (map pdevs_of_real n) of
    Some Xn ⇒ Inf_aform' (prec) (hd Xn))"
definition aform_sup_inner where "aform_sup_inner prec X n =
  (case inner_aforms' (prec) X (map pdevs_of_real n) of
    Some Xn ⇒ Sup_aform' (prec) (hd Xn))"

text ‹cannot fail›

lemma approx_un_ne_None: "approx_un p (λivl. Some (f ivl)) (Some r) ≠ None"
  by (auto simp: approx_un_def split_beta')

lemma approx_un_eq_Some:
  "approx_un p (λivl. Some (f ivl)) (Some r) = Some s ⟷
  s = ivl_err (real_interval (f (ivl_of_aform_err p r)))"
  using approx_un_ne_None
  by (auto simp: approx_un_def split_beta')

lemma is_float_uminus:
  "is_float aa ⟹ is_float (- aa)"
  by (auto simp: is_float_def)

lemma is_float_uminus_iff[simp]:
  "is_float (- aa) = is_float aa"
  using is_float_uminus[of aa] is_float_uminus[of "-aa"]
  by (auto simp: is_float_def)

lemma is_float_simps[simp]:
  "is_float 0"
  "is_float 1"
  "is_float (real_of_float x)"
  "is_float (truncate_down p X)"
  "is_float (truncate_up p X)"
  "is_float (eucl_truncate_down p X)"
  "is_float (eucl_truncate_up p X)"
  by (auto simp: is_float_def eucl_truncate_down_def eucl_truncate_up_def truncate_down_def truncate_up_def)

lemma is_float_sum_list'[simp]:
  "is_float (sum_list' p xs)"
  by (induction xs) (auto simp: is_float_def)

lemma is_float_inverse_aform': "is_float (fst (fst (inverse_aform' p X)))"
  unfolding inverse_aform'_def
  by (simp add: Let_def trunc_bound_eucl_def)

lemma is_float_min_range:
  "min_range_antimono p F D E X Y = Some ((a, b), c) ⟹ is_float a"
  "min_range_antimono p F D E X Y = Some ((a, b), c) ⟹ is_float c"
  "min_range_mono p F D E X Y = Some ((a, b), c) ⟹ is_float a"
  "min_range_mono p F D E X Y = Some ((a, b), c) ⟹ is_float c"
  by (auto simp: min_range_antimono_def min_range_mono_def Let_def bind_eq_Some_conv
      prod_eq_iff trunc_bound_eucl_def
      mid_err_def affine_unop_def split: prod.splits)

lemma is_float_ivl_err:
  assumes "ivl_err x = ((a, b), c)" "is_float (lower x)" "is_float (upper x)"
  shows "is_float a" "is_float c"
proof -
  define x1 where "x1 = lower x"
  define x2 where "x2 = upper x"
  have "a = (x1 + x2) / 2" "c = (x2 - x1) / 2"
    using assms by (auto simp: ivl_err_def x1_def x2_def)
  moreover have "(x1 + x2) / 2 ∈ float" "(x2 - x1) / 2 ∈ float"
    using assms
    by (auto simp: is_float_def x1_def x2_def)
  ultimately show "is_float a" "is_float c"
    unfolding is_float_def by blast+
qed

lemma is_float_trunc_bound_eucl[simp]:
  "is_float (fst (trunc_bound_eucl p X))"
  by (auto simp: trunc_bound_eucl_def Let_def)

lemma add_eq_times_2_in_float: "a + b = c * 2 ⟹ is_float a ⟹ is_float b ⟹ is_float c"
proof goal_cases
  case 1
  then have "c = (a + b) / 2"
    by simp
  also have "… ∈ float" using 1(3,2) by (simp add: is_float_def)
  finally show ?case by (simp add: is_float_def)
qed


lemma approx_floatarith_Some_is_float:
  "approx_floatarith p fa XS = Some ((a, b), ba) ⟹
  list_all (λ((a, b), c). is_float a ∧ is_float c) XS ⟹
  is_float a ∧ is_float ba"
  apply (induction fa arbitrary: a b ba)
  subgoal by (auto simp: bind_eq_Some_conv add_aform'_def Let_def )
  subgoal by (auto simp: bind_eq_Some_conv uminus_aform_def Let_def)
  subgoal by (auto simp: bind_eq_Some_conv mult_aform'_def Let_def )
  subgoal by (auto simp: bind_eq_Some_conv inverse_aform_err_def map_aform_err_def
        prod_eq_iff is_float_inverse_aform'
        acc_err_def aform_to_aform_err_def aform_err_to_aform_def inverse_aform_def
        Let_def split: if_splits)
  subgoal by (auto simp: bind_eq_Some_conv cos_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv arctan_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal apply (auto simp: bind_eq_Some_conv uminus_aform_def Let_def is_float_min_range
        is_float_ivl_err set_of_eq real_interval_abs
        split: if_splits prod.splits)
    apply (metis is_float_ivl_err(1) is_float_simps(3) lower_real_interval real_interval_abs upper_real_interval)
    by (metis is_float_ivl_err(2) is_float_simps(3) lower_real_interval real_interval_abs upper_real_interval)
  subgoal by (auto simp: bind_eq_Some_conv max_aform_err_def Let_def is_float_min_range
        is_float_ivl_err max_def
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv min_aform_err_def Let_def is_float_min_range
        is_float_ivl_err min_def
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv arctan_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv sqrt_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv exp_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv powr_aform_err_def approx_bin_def
        exp_aform_err_def Let_def is_float_min_range
        is_float_ivl_err mult_aform'_def
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv ln_aform_err_def Let_def is_float_min_range
        is_float_ivl_err
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv power_aform_err_def Let_def is_float_min_range
        is_float_ivl_err mid_err_def dest!: add_eq_times_2_in_float
        split: if_splits prod.splits)
  subgoal by (auto simp: bind_eq_Some_conv cos_aform_err_def Let_def is_float_min_range
        approx_un_def is_float_ivl_err
        split: if_splits)
  subgoal by (auto simp: bind_eq_Some_conv cos_aform_err_def Let_def is_float_min_range
        list_all_length
        split: if_splits)
  subgoal by (auto simp: bind_eq_Some_conv num_aform_def Let_def)
  done

lemma plain_floatarith_not_None:
  assumes "plain_floatarith N fa" "N ≤ length XS"
    "list_all (λ((a, b), c). is_float a ∧ is_float c) XS"
  shows "approx_floatarith p fa XS ≠ None"
  using assms
  by (induction fa)
    (auto simp: Let_def split_beta' approx_un_eq_Some prod_eq_iff
      approx_floatarith_Some_is_float)

lemma plain_floatarith_slp_not_None:
  assumes "⋀fa. fa ∈ set fas ⟹ plain_floatarith N fa" "N ≤ length XS"
    "list_all (λ((a, b), c). is_float a ∧ is_float c) XS"
  shows "approx_slp p fas XS ≠ None"
  using assms 
  apply (induction fas arbitrary: XS)
  subgoal by simp
  subgoal for fa fas XS
    using plain_floatarith_not_None[of N fa XS p]
    by (auto simp: Let_def split_beta' approx_un_eq_Some prod_eq_iff
      bind_eq_Some_conv approx_floatarith_Some_is_float)
  done

lemma plain_floatarithE:
  assumes "plain_floatarith N fa" "N ≤ length XS"
    "list_all (λ((a, b), c). is_float a ∧ is_float c) XS"
  obtains R where "approx_floatarith p fa XS = Some R"
  using plain_floatarith_not_None[OF assms]
  by force

lemma approx_slp_outer_eq_None_iff:
  "approx_slp_outer p a b XS = None ⟷
    approx_slp p b ((map (λx. (x, 0)) XS)) = None"
  by (auto simp: approx_slp_outer_def Let_def bind_eq_None_conv)

lemma approx_slp_sing_eq_None_iff:
  "approx_slp p [b] xs = None ⟷ approx_floatarith p b xs = None"
  by (auto simp: approx_slp_outer_def Let_def bind_eq_None_conv)

lemma plain_inner_floatariths:
  "plain_floatarith N (inner_floatariths xs ys)"
  if "list_all (plain_floatarith N) xs"  "list_all (plain_floatarith N) ys"
  using that
  by (induction xs ys rule: inner_floatariths.induct) auto

lemma aiN: "approx_floatarith p (inner_floatariths xs ys) zs ≠ None"
  if "⋀x. x ∈ set xs ⟹ approx_floatarith p x zs ≠ None"
  and "⋀x. x ∈ set ys ⟹ approx_floatarith p x zs ≠ None"
  using that
  apply (induction xs ys rule: inner_floatariths.induct)
  apply (auto simp: Let_def bind_eq_Some_conv)
  by (metis old.prod.exhaust)

lemma aiVN: "approx_floatarith p
        (inner_floatariths (map floatarith.Var [0..<length a])
          (map floatarith.Var [length a..<length a + length b]))
        (map (λx. (x, 0)) a @ map (λx. (x, 0)) b) ≠ None"
  by (rule aiN) (auto simp: nth_append)

lemma iaN: "inner_aforms' p a b ≠ None"
  unfolding inner_aforms'_def Let_def approx_slp_outer_def
  using aiVN[of p a b]
  by (auto simp: Let_def bind_eq_Some_conv)

definition "support_aform prec b X = (let ia = inner_aform X b in fst X ∙ b + tdev' (prec) (snd ia))"

definition "width_aforms prec X =
  (let t = tdev' (prec) ((abssum_of_pdevs_list (map snd X))) in t)"

definition "inf_aforms prec xs = map (Inf_aform' (prec)) xs"
definition "sup_aforms prec xs = map (Sup_aform' (prec)) xs"

definition "fresh_index_aforms xs = Max (insert 0 (degree_aform ` set xs))"

definition "independent_aforms x env = (msum_aform (fresh_index_aforms env) (0, zero_pdevs) x#env)"

definition "aform_form_ivl prec form xs =
  dRETURN (approx_form prec form (ivls_of_aforms prec xs) [])"
definition "aform_form prec form xs =
  dRETURN (approx_form_aform prec form xs)"
definition "aform_slp prec n slp xs =
  dRETURN (approx_slp_outer prec n slp xs)"
definition
  "aform_isFDERIV prec n ns fas xs =
    dRETURN (isFDERIV_approx prec n ns fas (ivls_of_aforms prec xs))"

definition aform_rel_internal: "aform_rel R = br Affine top O ⟨R⟩set_rel"
lemma aform_rel_def: "⟨rnv_rel⟩aform_rel = br Affine top"
  unfolding relAPP_def
  by (auto simp: aform_rel_internal)
definition "aforms_rel = br Joints top"

definition aform_rell :: "((real × real pdevs) list × real list set) set"
where "aform_rell = aforms_rel"

definition aforms_relp_internal: "aforms_relp R = aforms_rel O ⟨R⟩set_rel"
lemma aforms_relp_def: "⟨R⟩aforms_relp = aforms_rel O ⟨R⟩set_rel"
  by (auto simp: aforms_relp_internal relAPP_def)

definition "product_aforms x y = x @ msum_aforms (degree_aforms (x)) (replicate (length y) (0, zero_pdevs)) y"

lemma
  eucl_of_list_mem_eucl_of_list_aform:
  assumes "x ∈ Joints a"
  assumes "length a = DIM('a::executable_euclidean_space)"
  shows "eucl_of_list x ∈ Affine (eucl_of_list_aform a::'a aform)"
  using assms
  by (auto simp: Joints_def Affine_def valuate_def eucl_of_list_aform_def
      aform_val_def pdevs_val_sum_pdevs inner_simps eucl_of_list_inner
      intro!: euclidean_eqI[where 'a='a])

lemma
  in_image_eucl_of_list_eucl_of_list_aform:
  assumes "length x = DIM('a::executable_euclidean_space)" "xa ∈ Affine (eucl_of_list_aform x::'a aform)"
  shows "xa ∈ eucl_of_list ` Joints x"
  using assms
  by (auto intro!: nth_equalityI image_eqI[where x="list_of_eucl xa"]
      simp: aform_val_def eucl_of_list_aform_def Affine_def valuate_def Joints_def
      inner_simps pdevs_val_sum_pdevs eucl_of_list_inner)

lemma
  eucl_of_list_image_Joints:
  assumes "length x = DIM('a::executable_euclidean_space)"
  shows "eucl_of_list ` Joints x = Affine (eucl_of_list_aform x::'a aform)"
  using assms
  by (auto intro!: eucl_of_list_mem_eucl_of_list_aform
      in_image_eucl_of_list_eucl_of_list_aform)

definition "aform_ops =
⦇
  appr_of_ivl = aforms_of_ivls,
  product_appr = product_aforms,
  msum_appr = msum_aforms',
  inf_of_appr = λoptns. inf_aforms (precision optns),
  sup_of_appr = λoptns. sup_aforms (precision optns),
  split_appr = split_aforms_largest_uncond_take,
  appr_inf_inner = λoptns. aform_inf_inner (precision optns),
  appr_sup_inner = λoptns. aform_sup_inner (precision optns),
  inter_appr_plane = λoptns xs sctn. inter_aform_plane_lv (length xs) (precision optns) xs sctn,
  reduce_appr = λoptns. reduce_aforms (precision optns),
  width_appr = λoptns. width_aforms (precision optns),
  approx_slp_dres = λoptns. aform_slp (precision optns),
  approx_euclarithform = λoptns. aform_form (precision optns),
  approx_isFDERIV = λoptns. aform_isFDERIV  (precision optns)
⦈"

lemma Affine_eq_permI:
  assumes "fst X = fst Y"
  assumes "map snd (list_of_pdevs (snd X)) <~~> map snd (list_of_pdevs (snd Y))" (is "?perm X Y")
  shows "Affine X = Affine Y"
proof -
  {
    fix X Y and e::"nat ⇒ real"
    assume perm: "?perm X Y" and e: "e ∈ funcset UNIV {- 1..1}"
    from pdevs_val_of_list_of_pdevs2[OF e, of "snd X"]
    obtain e' where e':
      "pdevs_val e (snd X) = pdevs_val e' (pdevs_of_list (map snd (list_of_pdevs (snd X))))"
      "e' ∈ funcset UNIV {- 1..1}"
      by auto
    note e'(1)
    also from pdevs_val_perm[OF perm e'(2)]
    obtain e'' where e'':
      "e'' ∈ funcset UNIV {- 1..1}"
      "pdevs_val e' (pdevs_of_list (map snd (list_of_pdevs (snd X)))) =
      pdevs_val e'' (pdevs_of_list (map snd (list_of_pdevs (snd Y))))"
      by auto
    note e''(2)
    also from pdevs_val_of_list_of_pdevs[OF e''(1), of "snd Y", simplified]
    obtain e''' where e''':
      "pdevs_val e'' (pdevs_of_list (map snd (list_of_pdevs (snd Y)))) = pdevs_val e''' (snd Y)"
      "e''' ∈ funcset UNIV {- 1..1}"
      by auto
    note e'''(1)
    finally have "∃e' ∈ funcset UNIV {-1 .. 1}. pdevs_val e (snd X) = pdevs_val e' (snd Y)"
      using e'''(2) by auto
  } note E = this
  note e1 = E[OF assms(2)]
    and e2 = E[OF perm_sym[OF assms(2)]]
  show ?thesis
    by (auto simp: Affine_def valuate_def aform_val_def assms dest: e1 e2)
qed

context includes autoref_syntax begin

lemma aform_of_ivl_refine: "x ≤ y ⟹ (aform_of_ivl x y, atLeastAtMost x y) ∈ ⟨rnv_rel⟩aform_rel"
  by (auto simp: aform_rel_def br_def Affine_aform_of_ivl)

lemma aforms_of_ivl_leI1:
  fixes en::real
  assumes "-1 ≤ en" "xsn ≤ ysn"
  shows "xsn ≤ (xsn + ysn) / 2 + (ysn - xsn) * en / 2"
proof -
  have "xsn * (1 + en) ≤ ysn * (1 + en)"
    using assms mult_right_mono by force
  then show ?thesis
    by (auto simp: algebra_simps divide_simps)
qed

lemma aforms_of_ivl_leI2:
  fixes en::real
  assumes "1 ≥ en" "xsn ≤ ysn"
  shows "(xsn + ysn) / 2 + (ysn - xsn) * en / 2 ≤ ysn"
proof -
  have "xsn * (1 - en) ≤ ysn * (1 - en)"
    using assms mult_right_mono by force
  then show ?thesis
    by (auto simp: algebra_simps divide_simps)
qed

lemma Joints_aforms_of_ivlsD1:
  "zs ∈ Joints (aforms_of_ivls xs ys) ⟹ list_all2 (≤) xs ys ⟹ list_all2 (≤) xs zs"
  by (auto simp: Joints_def valuate_def aforms_of_ivls_def aform_val_def Pi_iff
      list_all2_conv_all_nth intro!: list_all2_all_nthI aforms_of_ivl_leI1)

lemma Joints_aforms_of_ivlsD2:
  "zs ∈ Joints (aforms_of_ivls xs ys) ⟹ list_all2 (≤) xs ys ⟹ list_all2 (≤) zs ys"
  by (auto simp: Joints_def valuate_def aforms_of_ivls_def aform_val_def Pi_iff
      list_all2_conv_all_nth intro!: list_all2_all_nthI aforms_of_ivl_leI2)

lemma aforms_of_ivls_refine:
  "list_all2 (≤) xrs yrs ⟹
       (xri, xrs) ∈ ⟨rnv_rel⟩list_rel ⟹
       (yri, yrs) ∈ ⟨rnv_rel⟩list_rel ⟹ (aforms_of_ivls xri yri, lv_ivl xrs yrs) ∈ aforms_rel"
  apply (auto simp: aforms_rel_def br_def list_all2_lengthD lv_ivl_def
      Joints_aforms_of_ivlsD1 Joints_aforms_of_ivlsD2 intro!: aforms_of_ivls)
  subgoal by (simp add: list_all2_conv_all_nth)
  subgoal by (simp add: list_all2_conv_all_nth)
  done

lemma degrees_zero_pdevs[simp]: "degrees (replicate n zero_pdevs) = 0"
  by (auto simp: degrees_def intro!: Max_eqI)

lemma Joints_product_aforms:
  "Joints (product_aforms a b) = product_listset (Joints a) (Joints b)"
  apply (auto simp: Joints_def valuate_def product_listset_def product_aforms_def
      aform_vals_def[symmetric] aform_val_msum_aforms)
  subgoal for e
    apply (rule image_eqI[where
          x="(aform_vals e a,
              List.map2 (+) (aform_vals e (replicate (length b) (0, zero_pdevs))) (aform_vals (λi. e (i + degree_aforms a)) b))"])
     apply (auto simp: split_beta')
    apply (auto simp: aform_vals_def intro!: nth_equalityI image_eqI[where x="λi. e (i + degree_aforms a)"])
    done
  subgoal for e1 e2
    apply (rule image_eqI[where x="λi. if i < degree_aforms a then e1 i else e2 (i - degree_aforms a)"])
     apply (auto simp: aform_vals_def aform_val_def Pi_iff
        intro!: nth_equalityI pdevs_val_degree_cong)
    subgoal premises prems for x xs k
    proof -
      from prems
      have "degree xs ≤ degree_aforms a"
        by (auto simp: degrees_def Max_gr_iff)
      then show ?thesis using prems by auto
    qed
    done
  done

lemma product_aforms_refine:
  "(product_aforms, product_listset) ∈ aforms_rel → aforms_rel → aforms_rel"
  by (auto simp: aforms_rel_def br_def Joints_product_aforms)

lemma mem_lv_rel_set_rel_iff:
  fixes z::"'a::executable_euclidean_space set"
  shows "(y, z) ∈ ⟨lv_rel⟩set_rel ⟷ (z = eucl_of_list ` y ∧ (∀x ∈ y. length x = DIM('a)))"
  unfolding lv_rel_def
  by (auto simp: set_rel_def br_def)

lemma eucl_of_list_mem_lv_rel: "length x = DIM('a::executable_euclidean_space) ⟹
    (x, eucl_of_list x::'a) ∈ lv_rel"
  unfolding lv_rel_def
  by (auto simp:  br_def)

lemma
  mem_Joints_msum_aforms'I:
  "a ∈ Joints x ⟹ b ∈ Joints y ⟹ List.map2 (+) a b ∈ Joints (msum_aforms' x y)"
  by (auto simp: Joints_msum_aforms degrees_def)

lemma
  mem_Joints_msum_aforms'E:
  assumes "xa ∈ Joints (msum_aforms' x y)" 
  obtains a b where "xa = List.map2 (+) a b" "a ∈ Joints x" "b ∈ Joints y"
  using assms
  by (auto simp: Joints_msum_aforms degrees_def)

lemma msum_aforms'_refine_raw:
  shows "(msum_aforms' x y, {List.map2 (+) a b|a b. a ∈ Joints x ∧ b ∈ Joints y}) ∈ aforms_rel"
  unfolding aforms_rel_def br_def
  by (safe elim!: mem_Joints_msum_aforms'E intro!: mem_Joints_msum_aforms'I) (auto simp: Joints_imp_length_eq)

lemma aforms_relD: "(a, b) ∈ aforms_rel ⟹ b = Joints a"
  by (auto simp: aforms_rel_def br_def)

lemma msum_aforms'_refine:
  "(msum_aforms', λxs ys. {List.map2 (+) x y |x y. x ∈ xs ∧ y ∈ ys}) ∈ aforms_rel → aforms_rel → aforms_rel"
  by (safe dest!: aforms_relD intro!: msum_aforms'_refine_raw)

lemma length_inf_aforms[simp]: "length (inf_aforms optns x) = length x"
  and length_sup_aforms[simp]: "length (sup_aforms optns x) = length x"
  by (auto simp: inf_aforms_def sup_aforms_def)

lemma inf_aforms_refine:
  "(xi, x) ∈ aforms_rel ⟹ length xi = d ⟹ (RETURN (inf_aforms optns xi), Inf_specs d x) ∈ ⟨rl_rel⟩nres_rel"
  unfolding o_def
  apply (auto simp: br_def aforms_rel_def mem_lv_rel_set_rel_iff nres_rel_def Inf_specs_def
      intro!: RETURN_SPEC_refine)
  unfolding lv_rel_def
  by (auto simp: aforms_rel_def br_def Joints_imp_length_eq
      eucl_of_list_inner inf_aforms_def nth_in_AffineI
      intro!: Inf_aform'_Affine_le list_all2_all_nthI)

lemma sup_aforms_refine:
  "(xi, x) ∈ aforms_rel ⟹ length xi = d ⟹ (RETURN (sup_aforms optns xi), Sup_specs d x) ∈ ⟨rl_rel⟩nres_rel"
  unfolding o_def
  apply (auto simp: aforms_relp_def mem_lv_rel_set_rel_iff nres_rel_def Sup_specs_def
      intro!: RETURN_SPEC_refine)
  unfolding lv_rel_def
  by (auto simp: aforms_rel_def br_def Joints_imp_length_eq
      eucl_of_list_inner sup_aforms_def nth_in_AffineI
      intro!: Sup_aform'_Affine_ge list_all2_all_nthI)

lemma nres_of_THE_DRES_le:
  assumes "⋀v. x = Some v ⟹ RETURN v ≤ y"
  shows "nres_of (THE_DRES x) ≤ y"
  using assms by (auto simp: THE_DRES_def split: option.split)

lemma degree_le_fresh_index: "a ∈ set A ⟹ degree_aform a ≤ fresh_index_aforms A"
  by (auto simp: fresh_index_aforms_def intro!: Max_ge)

lemma zero_in_JointsI: "xs ∈ Joints XS ⟹ z = (0, zero_pdevs) ⟹ 0 # xs ∈ Joints (z # XS)"
  by (auto simp: Joints_def valuate_def)

lemma cancel_nonneg_pos_add_multI: "0 ≤ c + c * x"
  if "c ≥ 0" "1 + x ≥ 0"
  for c x::real
proof -
  have "0 ≤ c + c * x ⟷ 0 ≤ c * (1 + x)" by (auto simp: algebra_simps)
  also have "… ⟷ 0 ≤ 1 + x"
    using that
    by (auto simp: zero_le_mult_iff)
  finally show ?thesis 
    using that by simp
qed
   
lemma Joints_map_split_aform:
  fixes x::"real aform list"
  shows "Joints x ⊆ Joints (map (λa. fst (split_aform a i)) x) ∪ Joints (map (λb. snd (split_aform b i)) x)"
  unfolding subset_iff
  apply (intro allI impI)
proof goal_cases
  case (1 xs)
  then obtain e where e: "e ∈ funcset UNIV {-1 .. 1}" "xs = map (aform_val e) x"
    by (auto simp: Joints_def valuate_def)
  consider "e i ≥ 0" | "e i ≤ 0" by arith
  then show ?case
  proof cases
    case 1
    let ?e = "e(i:= 2 * e i - 1)"
    note e(2)
    also
    have "map (aform_val e) x = map (aform_val ?e) (map (λb. snd (split_aform b i)) x)"
      by (auto simp: aform_val_def split_aform_def Let_def divide_simps algebra_simps)
    also have "… ∈ Joints (map (λb. snd (split_aform b i)) x)"
      using e ‹0 ≤ e i›
      by (auto simp: Joints_def valuate_def Pi_iff intro!: image_eqI[where x = ?e])
    finally show ?thesis by simp
  next
    case 2
    let ?e = "e(i:= 2 * e i + 1)"
    note e(2)
    also
    have "map (aform_val e) x = map (aform_val ?e) (map (λb. fst (split_aform b i)) x)"
      by (auto simp: aform_val_def split_aform_def Let_def divide_simps algebra_simps)
    also have "… ∈ Joints (map (λb. fst (split_aform b i)) x)"
      using e ‹0 ≥ e i›
      by (auto simp: Joints_def valuate_def Pi_iff real_0_le_add_iff
          intro!: image_eqI[where x = ?e] cancel_nonneg_pos_add_multI)
    finally show ?thesis by simp
  qed
qed


lemma split_aforms_lemma:
  fixes x::"real aform list"
  assumes "(x, y) ∈ aforms_rel"
  assumes "z ⊆ y"
  assumes "l = (length x)"
  shows "∃a b. (split_aforms x i, a, b)
     ∈ aforms_rel ×r aforms_rel ∧ env_len a l ∧ env_len b l ∧ z ⊆ a ∪ b"
  using assms
  apply (auto simp: split_aforms_def o_def)
  apply (rule exI[where x="Joints (map (λx. fst (split_aform x i)) x)"])
  apply auto
  subgoal by (auto intro!: brI simp: aforms_rel_def)
  apply (rule exI[where x="Joints (map (λx. snd (split_aform x i)) x)"])
  apply (rule conjI)
  subgoal by (auto intro!: brI simp: aforms_rel_def)
  subgoal
    using Joints_map_split_aform[of x i]
    by (auto simp: br_def aforms_rel_def env_len_def Joints_imp_length_eq)
  done

lemma length_summarize_pdevs_list[simp]:
  "length (summarize_pdevs_list a b c d) = length d"
  by (auto simp: summarize_pdevs_list_def)

lemma length_summarize_aforms[simp]:
  "length (summarize_aforms a b e d) = length d"
  by (auto simp: summarize_aforms_def Let_def)

lemma
  split_aform_largest_take_refine:
  "(ni, n) ∈ nat_rel ⟹
       (xi::real aform list, x) ∈ aforms_rel ⟹
       length xi = d ⟹ (RETURN (split_aforms_largest_uncond_take ni xi), split_spec_params d n x) ∈ ⟨aforms_rel ×r aforms_rel⟩nres_rel"
  apply (auto simp: split_spec_params_def nres_rel_def aforms_relp_def mem_lv_rel_set_rel_iff
      split_aforms_largest_uncond_take_def Let_def comps
      split: prod.splits
      intro!: RETURN_SPEC_refine)
  apply (rule split_aforms_lemma)
  by (auto simp add: aforms_rel_def)

lemma aform_val_pdevs_of_real[simp]: "aform_val e (pdevs_of_real x) = x"
  by (auto simp: pdevs_of_real_def)

lemma degree_aform_pdevs_of_real[simp]: "degree_aform (pdevs_of_real x) = 0"
  by (auto simp: pdevs_of_real_def)

lemma interpret_floatarith_inner_eq:
  shows "interpret_floatarith (inner_floatariths xs ys) vs =
    (∑(x, y) ← (zip xs ys). (interpret_floatarith x vs) * (interpret_floatarith y vs))"
  by (induction xs ys rule: inner_floatariths.induct) auto

lemma approx_slp_outer_sing:
  "approx_slp_outer p (Suc 0) [fa] XS = Some R ⟷ (∃Y.
  approx_floatarith p fa (map (λx. (x, 0)) XS) = Some Y ∧
       [aform_err_to_aform Y (max (degree_aforms XS) (degree_aform_err Y))] = R)"
  by (auto simp: approx_slp_outer_def bind_eq_Some_conv degrees_def)

lemma aforms_err_aform_valsI:
  assumes "vs = aform_vals e XS"
  shows "vs ∈ aforms_err e (map (λx. (x, 0)) (XS))"
    by (auto simp: assms aforms_err_def o_def aform_err_def aform_vals_def)


lemma aform_val_degree_cong:
  "b = d ⟹ (⋀i. i < degree_aform d ⟹ a i = c i) ⟹ aform_val a b = aform_val c d"
  by (auto simp: aform_val_def intro!: pdevs_val_degree_cong)

lemma mem_degree_aformD: "x ∈ set XS ⟹ degree_aform x ≤ degree_aforms XS"
  by (auto simp: degrees_def)
lemma degrees_append_leD1: "(degrees xs) ≤ degrees (xs @ ys)"
  unfolding degrees_def
  by (rule Max_mono) (auto simp: degrees_def min_def max_def Max_ge_iff image_Un Max_gr_iff)

lemma inner_aforms':
  assumes "xs ∈ Joints XS"
  assumes "inner_aforms' p XS (map pdevs_of_real rs) = Some R"
  shows "(∑(x, y) ← (zip xs rs). x * y) ∈ Affine (hd R)" (is ?th1) "length R = 1" (is ?th2)
proof -
  from assms obtain e where "e ∈ funcset UNIV {-1 .. 1}" "xs = aform_vals e XS"
    by (auto simp: Joints_def valuate_def aform_vals_def)
  then have e: "xs @ rs = aform_vals e (XS @ map pdevs_of_real rs)"
    "xs @ rs ∈ Joints (XS @ map pdevs_of_real rs)"
    "length xs = length XS"
    "e ∈ funcset UNIV {- 1..1}"
    by (auto simp: aform_vals_def o_def degrees_def Joints_def valuate_def)
  have "approx_slp_outer p (Suc 0)
     ([inner_floatariths (map floatarith.Var [0..<length XS])
         (map floatarith.Var [length XS..<length XS + length rs])])
     (XS @ map pdevs_of_real rs) =
    Some R"
    using assms(2)
    by (auto simp: inner_aforms'_def)
  then obtain Y where Y:
    "approx_floatarith p
       (inner_floatariths (map floatarith.Var [0..<length XS])
         (map floatarith.Var [length XS..<length XS + length rs]))
       (map (λx. (x, 0)) (XS @ map pdevs_of_real rs)) =
      Some Y"
    "R = [aform_err_to_aform Y (max (degree_aforms (XS @ map pdevs_of_real rs)) (degree_aform_err Y))]"
    unfolding approx_slp_outer_sing by auto
  let ?m = "(max (degree_aforms (XS @ map pdevs_of_real rs)) (degree_aform_err Y))"
  from approx_floatarith_Elem[OF Y(1) e(4) aforms_err_aform_valsI[OF e(1)]]
  have "interpret_floatarith
     (inner_floatariths (map floatarith.Var [0..<length XS]) (map floatarith.Var [length XS..<length XS + length rs]))
     (xs @ rs)
    ∈ aform_err e Y" "degree_aform_err Y ≤ max (degree_aforms (XS @ map pdevs_of_real rs)) (degree_aform_err Y)"
    by auto
  from aform_err_to_aformE[OF this] obtain err where err:
    "interpret_floatarith
           (inner_floatariths (map floatarith.Var [0..<length XS])
             (map floatarith.Var [length XS..<length XS + length rs]))
           (xs @ rs) =
          aform_val (e(max (degree_aforms (XS @ map pdevs_of_real rs)) (degree_aform_err Y) := err))
           (aform_err_to_aform Y (max (degree_aforms (XS @ map pdevs_of_real rs)) (degree_aform_err Y)))"
    "- 1 ≤ err" "err ≤ 1"
    by auto
  let ?e' = "(e(max (degrees (map snd XS @ map (snd ∘ pdevs_of_real) rs)) (degree_aform_err Y) := err))"
  from e(1) have e': "xs @ rs = aform_vals ?e' (XS @ map pdevs_of_real rs)"
    apply (auto simp: aform_vals_def intro!: aform_val_degree_cong)
     apply (frule mem_degree_aformD)
     apply (frule le_less_trans[OF degrees_append_leD1])
     apply (auto simp: o_def)
    done
  from err have "interpret_floatarith
           (inner_floatariths (map floatarith.Var [0..<length XS])
             (map floatarith.Var [length XS..<length XS + length rs]))
           (xs @ rs)#xs@rs ∈ Joints (R @ XS @ map pdevs_of_real rs)"
    using e(1,3,4) e'
    apply (auto simp: valuate_def Joints_def
      intro!: nth_equalityI
        image_eqI[where x="?e'"])
     apply (simp add: Y aform_vals_def; fail)
    apply (simp add: Y o_def)
    apply (simp add: nth_append nth_Cons)
    apply (auto split: nat.splits simp: nth_append nth_Cons aform_vals_def)
    done
  then have "(∑(x, y)←zip (map floatarith.Var [0..<length XS])
                (map floatarith.Var
                  [length XS..<
                   length XS + length rs]). interpret_floatarith x (xs @ rs) * interpret_floatarith y (xs @ rs)) #
    xs @ rs
    ∈ Joints (R @ XS @ map pdevs_of_real rs)"
    apply (subst (asm)interpret_floatarith_inner_eq )
    apply (auto simp:  )
    done
  also have "(∑(x, y)←zip (map floatarith.Var [0..<length XS])
                (map floatarith.Var
                  [length XS..<
                   length XS + length rs]). interpret_floatarith x (xs @ rs) * interpret_floatarith y (xs @ rs)) =
    (∑(x, y)←zip xs rs. x * y)"
    by (auto simp: sum_list_sum_nth assms e(3) nth_append intro!: sum.cong)
  finally show ?th1 ?th2
    by (auto simp: Affine_def valuate_def Joints_def Y)
qed

lemma inner_aforms'_inner_lv_rel:
  "(a, a') ∈ aforms_rel ⟹
       inner_aforms' prec a (map pdevs_of_real a'a) = Some R ⟹
       x ∈ a' ⟹ inner_lv_rel x a'a ∈ Affine (hd R)"
  unfolding mem_lv_rel_set_rel_iff
  unfolding lv_rel_def aforms_rel_def
  apply (auto simp: br_def)
  apply (subst arg_cong2[where f="(∈)", OF _ refl])
  defer
   apply (rule inner_aforms')
     apply (auto simp: br_def Joints_imp_length_eq inner_lv_rel_def)
  done

lemma aform_inf_inner_refine:
  "(RETURN o2 aform_inf_inner optns, Inf_inners) ∈ aforms_rel → rl_rel → ⟨rnv_rel⟩nres_rel"
  by (auto simp: aforms_relp_def nres_rel_def Inf_inners_def aform_inf_inner_def[abs_def] comp2_def
      iaN
      intro!: Inf_aform'_Affine_le inner_aforms'_inner_lv_rel
      split: option.splits list.splits)

lemma aform_sup_inner_refine:
  "(RETURN o2 aform_sup_inner optns, Sup_inners) ∈ aforms_rel → rl_rel → ⟨rnv_rel⟩nres_rel"
  by (auto simp: aforms_relp_def nres_rel_def Sup_inners_def aform_sup_inner_def[abs_def] comp2_def
      iaN
      intro!: Sup_aform'_Affine_ge inner_aforms'_inner_lv_rel
      split: option.splits list.splits)

lemma lv_aforms_rel_comp_br_Affine_le: "lv_aforms_rel O br Affine top ⊆ ⟨lv_rel⟩aforms_relp"
  apply (auto simp: lv_aforms_rel_def aforms_relp_def br_def)
  apply (rule relcompI)
   apply (auto simp: aforms_rel_def intro!: brI )
  by (auto simp: mem_lv_rel_set_rel_iff Joints_imp_length_eq
      intro!: eucl_of_list_mem_eucl_of_list_aform intro!: in_image_eucl_of_list_eucl_of_list_aform)

lemma bijective_lv_rel[relator_props]: "bijective lv_rel"
  unfolding lv_rel_def bijective_def
  apply (auto simp: br_def)
  by (metis eucl_of_list_inj)

lemma sv_lv_rel_inverse[relator_props]: "single_valued (lv_rel¯)"
  using bijective_lv_rel
  by (rule bijective_imp_sv)

lemma list_of_eucl_image_lv_rel_inverse:
  "(x, list_of_eucl ` x) ∈ ⟨lv_rel¯⟩set_rel"
  unfolding set_rel_sv[OF sv_lv_rel_inverse]
  apply auto
    apply (rule ImageI)
     apply (rule converseI)
     apply (rule lv_relI)
     apply auto
   apply (rule image_eqI)
    prefer 2 apply assumption
  unfolding lv_rel_def
   apply (auto simp: br_def)
  subgoal for y
    apply (rule exI[where x="list_of_eucl y"])
    apply auto
    done
  done

lemma lv_rel_comp_lv_rel_inverse:
  "((lv_rel::(_×'a::executable_euclidean_space) set) O lv_rel¯) = {(x, y). x = y ∧ length x = DIM('a)}"
  apply (auto simp: intro!: lv_relI)
  unfolding lv_rel_def
  by (auto simp: br_def intro!: eucl_of_list_inj)

lemma inter_aform_plane_refine_aux:
  "d = CARD('n::enum) ⟹
  (xi, x) ∈ aforms_rel ⟹ (si, s) ∈ ⟨rl_rel⟩sctn_rel ⟹
  length xi = d ⟹
  length (normal si) = d ⟹
  (nres_of (inter_aform_plane_lv (length xi) optns xi si), inter_sctn_specs d x s)
  ∈ ⟨aforms_rel⟩nres_rel"
proof (goal_cases)
  case 1
  from 1 have sis: "si = s"
    by (cases si; cases s) (auto simp: sctn_rel_def)
  have Dp: "DIM_precond TYPE('n rvec) CARD('n)"
    by auto
  have a: "(xi, eucl_of_list_aform xi::'n rvec aform) ∈ lv_aforms_rel"
    by (auto simp: lv_aforms_rel_def br_def aforms_relp_def aforms_rel_def mem_lv_rel_set_rel_iff 1
        Joints_imp_length_eq)
  have b: "(si, (Sctn (eucl_of_list (normal si)) (pstn si))::'n rvec sctn) ∈ ⟨lv_rel⟩sctn_rel"
    using 1
    by (cases si)
      (auto simp: lv_aforms_rel_def br_def aforms_relp_def aforms_rel_def mem_lv_rel_set_rel_iff
        Joints_imp_length_eq sctn_rel_def intro!: lv_relI)
  have a: "(nres_of (inter_aform_plane_lv CARD('n) optns xi si),
    inter_aform_plane optns (eucl_of_list_aform xi::'n rvec aform) (Sctn (eucl_of_list (normal si)) (pstn si)))
    ∈ ⟨lv_aforms_rel⟩nres_rel"
    (is "(_, inter_aform_plane _ ?ea ?se) ∈ _")
    using inter_aform_plane_lv.refine[OF Dp IdI a b, of optns]
    by simp
  have b: "(inter_aform_plane optns ?ea ?se, inter_sctn_spec (Affine ?ea) ?se) ∈ ⟨br Affine top⟩nres_rel"
    using 1
    apply (auto simp: inter_sctn_spec_def nres_rel_def inter_aform_plane_def project_ncoord_aform_def
        inter_aform_plane_ortho_nres_def
        split: option.splits
        intro!: RETURN_SPEC_refine
        dest!: inter_inter_aform_plane_ortho)
     apply (auto simp: aform_rel_def br_def nres_rel_def comps bind_eq_Some_conv
        inter_sctn_spec_def inter_aform_plane_def plane_of_def aforms_relp_def aforms_rel_def
        mem_lv_rel_set_rel_iff
        intro!: RETURN_SPEC_refine nres_of_THE_DRES_le mem_Affine_project_coord_aformI
        dest!: inter_inter_aform_plane_ortho
        split: if_splits)
     apply (auto dest!: mem_Affine_project_coord_aformD simp: abs_in_Basis_iff plane_of_def)
    done
  from relcompI[OF a b]
  have "(nres_of (inter_aform_plane_lv CARD('n) optns xi si), inter_sctn_spec (Affine ?ea) ?se) ∈
      ⟨lv_aforms_rel⟩nres_rel O ⟨br Affine top⟩nres_rel"
    by auto
  also have "⟨lv_aforms_rel⟩nres_rel O ⟨br Affine top⟩nres_rel ⊆ ⟨⟨lv_rel⟩aforms_relp⟩nres_rel"
    by (simp add: Autoref_Misc.nres_rel_comp lv_aforms_rel_comp_br_Affine_le nres_rel_mono)
  finally
  have step1: "(nres_of (inter_aform_plane_lv CARD('n) optns xi si),
      inter_sctn_spec (Affine (eucl_of_list_aform xi)::'n rvec set) (Sctn (eucl_of_list (normal si)) (pstn si)))
      ∈ ⟨⟨lv_rel⟩aforms_relp⟩nres_rel"
    by simp
  have step2: "(inter_sctn_spec (Affine (eucl_of_list_aform xi)::'n rvec set) (Sctn (eucl_of_list (normal si)) (pstn si)),
    inter_sctn_specs CARD('n) (Joints xi) si) ∈ ⟨⟨lv_rel¯⟩set_rel⟩nres_rel"
    apply (auto simp: inter_sctn_specs_def inter_sctn_spec_def simp: nres_rel_def intro!: SPEC_refine)
    subgoal for x
      apply (rule exI[where x="list_of_eucl ` x"])
      apply (auto simp: env_len_def plane_ofs_def intro: list_of_eucl_image_lv_rel_inverse)
      subgoal for y
        apply (rule image_eqI[where x="eucl_of_list y"])
          apply (subst list_of_eucl_eucl_of_list)
          apply (auto simp: 1 Joints_imp_length_eq)
        apply (rule subsetD, assumption)
        apply (auto simp: intro!: eucl_of_list_mem_eucl_of_list_aform 1)
        using inner_lv_rel_autoref[where 'a="'n rvec", param_fo, OF lv_relI lv_relI, of y "normal si"]
        by (auto simp: plane_of_def 1 Joints_imp_length_eq)
      subgoal for y
        apply (drule subsetD, assumption)
        using inner_lv_rel_autoref[where 'a="'n rvec", param_fo, OF lv_relI lv_relI, of "list_of_eucl y" "normal si"]
        by (auto simp: plane_of_def 1 Joints_imp_length_eq)
      done
    done
  from relcompI[OF step1 step2]
  have "(nres_of (inter_aform_plane_lv CARD('n) optns xi si), inter_sctn_specs CARD('n) (Joints xi) si)
    ∈ ⟨⟨lv_rel::(real list × 'n rvec)set⟩aforms_relp⟩nres_rel O ⟨⟨lv_rel¯⟩set_rel⟩nres_rel" by simp
  also have "⟨⟨lv_rel::(real list × 'n rvec)set⟩aforms_relp⟩nres_rel O ⟨⟨lv_rel¯⟩set_rel⟩nres_rel ⊆
    ⟨aforms_rel O Id⟩nres_rel"
    unfolding nres_rel_comp O_assoc aforms_relp_def
    apply (rule nres_rel_mono)
    apply (rule relcomp_mono[OF order_refl])
    unfolding set_rel_sv[OF sv_lv_rel_inverse] set_rel_sv[OF lv_rel_sv]
    apply (auto simp: length_lv_rel)
    unfolding relcomp_Image[symmetric] lv_rel_comp_lv_rel_inverse
     apply (auto simp: Basis_not_uminus length_lv_rel)
    unfolding lv_rel_def
    subgoal for a b c d
      apply (auto dest!: brD simp: length_lv_rel)
      using eucl_of_list_inj[where 'a="'n rvec", of d b]
      by auto
    done
  finally show ?case using 1 by (simp add: aforms_rel_def br_def sis)
qed

end

setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum"}, SOME @{typ "'a list"})›
setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum_all"}, SOME @{typ "('a ⇒ bool) ⇒ bool"})›
setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum_ex"}, SOME @{typ "('a ⇒ bool) ⇒ bool"})›

lemmas inter_aform_plane_refine_unoverloaded0 =
  inter_aform_plane_refine_aux[internalize_sort "'n::enum",
    unoverload enum_class.enum,
    unoverload enum_class.enum_all,
    unoverload enum_class.enum_ex]

theorem
  inter_aform_plane_refine_unoverloaded:
    "class.enum (enum::'a list) enum_all enum_ex ∧ d = CARD('a) ⟹
    (xi, x) ∈ aforms_rel ⟹
    (si, s) ∈ ⟨rl_rel⟩sctn_rel ⟹
    length xi = d ⟹
    length (normal si) = d ⟹
    (nres_of (inter_aform_plane_lv (length xi) optns xi si), inter_sctn_specs d x s)
    ∈ ⟨aforms_rel⟩nres_rel"
  by (rule inter_aform_plane_refine_unoverloaded0) auto

setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum"}, SOME @{typ "'a::enum list"})›
setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum_all"}, SOME @{typ "('a::enum ⇒ bool) ⇒ bool"})›
setup ‹Sign.add_const_constraint (@{const_name "enum_class.enum_ex"}, SOME @{typ "('a::enum ⇒ bool) ⇒ bool"})›

context includes autoref_syntax begin

text ‹TODO: this is a special case of ‹Cancel_Card_Constraint› from ‹AFP/Perron_Frobenius›!›
lemma type_impl_card_n_enum:
  assumes "∃(Rep :: 'a ⇒ nat) Abs. type_definition Rep Abs {0 ..< n :: nat}"
  obtains enum enum_all enum_ex
  where "class.enum (enum::'a list) enum_all enum_ex ∧ n = CARD('a)"
proof -
  from assms obtain rep :: "'a ⇒ nat" and abs :: "nat ⇒ 'a"
    where t: "type_definition rep abs {0 ..<n}" by auto
  then interpret type_definition rep abs ‹{0 ..<n}› .
  have "card (UNIV :: 'a set) = card {0 ..< n}"
    by (rule card)
  also have "… = n" by auto
  finally have bn: "CARD ('a) = n" .

  let ?enum = "(map abs [0..<n])"
  have "class.enum ?enum (Ball (set ?enum)) (Bex (set ?enum))"
    by standard (auto simp: distinct_map Abs_inject simp flip: Abs_image intro: inj_onI)
  with bn have "class.enum ?enum (Ball (set ?enum)) (Bex (set ?enum)) ∧ n = CARD('a)"
    by simp
  then show ?thesis ..
qed


lemma inter_aform_plane_refine_ex_typedef:
  "(xi, x) ∈ aforms_rel ⟹
  (si, s) ∈ ⟨rl_rel⟩sctn_rel ⟹
  length xi = d ⟹
  length (normal si) = d ⟹
  (nres_of (inter_aform_plane_lv (length xi) optns xi si), inter_sctn_specs d x s)
  ∈ ⟨aforms_rel⟩nres_rel"
  if "∃(Rep :: 'a ⇒ nat) Abs. type_definition Rep Abs {0 ..< d :: nat}"
  by (rule type_impl_card_n_enum[OF that]) (rule inter_aform_plane_refine_unoverloaded; assumption)

lemma inter_aform_plane_refine:
  "0 < d ⟹
  (xi, x) ∈ aforms_rel ⟹
  (si, s) ∈ ⟨Id⟩sctn_rel ⟹
  length xi = d ⟹
  length (normal si) = d ⟹
  (nres_of (inter_aform_plane_lv (length xi) optns xi si), inter_sctn_specs d x s)
  ∈ ⟨aforms_rel⟩nres_rel"
  by (rule inter_aform_plane_refine_ex_typedef[cancel_type_definition, simplified])

lemma Joints_reduce_aforms: "x ∈ Joints X ⟹ x ∈ Joints (reduce_aforms prec t X)"
proof (auto simp: reduce_aforms_def summarize_threshold_def[abs_def] Joints_def valuate_def aform_val_def, goal_cases)
  case (1 e)
  from summarize_aformsE[OF ‹e ∈ _› order_refl]
  obtain e' where
    "aform_vals e X = aform_vals e' (summarize_aforms prec t (degree_aforms X) X)"
    "⋀i. i < degree_aforms X ⟹ e i = e' i"
    "e' ∈ funcset UNIV {- 1..1}"
    by blast
  thus ?case
    by (auto intro!: image_eqI[where x=e'] simp: aform_vals_def)
qed

lemma length_reduce_aform[simp]: "length (reduce_aforms optns a x) = length x"
  by (auto simp: reduce_aforms_def)

lemma reduce_aform_refine:
  "(xi, x) ∈ aforms_rel ⟹ length xi = d ⟹
    (RETURN (reduce_aforms prec C xi), reduce_specs d r x) ∈ ⟨aforms_rel⟩nres_rel"
  apply (auto simp: reduce_specs_def nres_rel_def comps aforms_relp_def mem_lv_rel_set_rel_iff
      aforms_rel_def env_len_def
      intro!: RETURN_SPEC_refine)
  apply (auto simp: br_def env_len_def)
  by (auto simp: mem_lv_rel_set_rel_iff Joints_imp_length_eq
      intro!: in_image_eucl_of_list_eucl_of_list_aform Joints_reduce_aforms
      eucl_of_list_mem_eucl_of_list_aform)

lemma aform_euclarithform_refine:
  "(nres_of o2 aform_form optns, approx_form_spec) ∈ Id → aforms_rel → ⟨bool_rel⟩nres_rel"
  by (auto simp: approx_form_spec_def nres_rel_def comps aform_form_def
      aforms_rel_def br_def approx_form_aform
      dest!: approx_form_aux
      intro!: ivls_of_aforms)

lemma aform_isFDERIV:
  "(λN xs fas vs. nres_of (aform_isFDERIV optns N xs fas vs), isFDERIV_spec)
    ∈ nat_rel → ⟨nat_rel⟩list_rel → ⟨Id⟩list_rel → aforms_rel → ⟨bool_rel⟩nres_rel"
  by (auto simp: isFDERIV_spec_def nres_rel_def comps aform_isFDERIV_def
      aforms_rel_def br_def
      dest!: approx_form_aux
      intro!: ivls_of_aforms isFDERIV_approx)

lemma approx_slp_lengthD: "approx_slp p slp a = Some xs ⟹
  length xs = length slp + length a"
  by (induction slp arbitrary: xs a) (auto simp: bind_eq_Some_conv)

lemma approx_slp_outer_lengthD: "approx_slp_outer p d slp a = Some xs ⟹
  length xs = min d (length slp + length a)"
  by (auto simp: approx_slp_outer_def Let_def bind_eq_Some_conv o_def
      aforms_err_to_aforms_def aform_err_to_aform_def dest!: approx_slp_lengthD)

lemma approx_slp_refine:
  "(nres_of o3 aform_slp prec, approx_slp_spec fas) ∈ 
    nat_rel → fas_rel → aforms_rel → ⟨⟨aforms_rel⟩option_rel⟩nres_rel"
  apply (auto simp: approx_slp_spec_def comps aform_slp_def nres_rel_def
      intro!: RETURN_SPEC_refine ASSERT_refine_right)
  subgoal for a b
    apply (rule exI[where x = "map_option Joints (approx_slp_outer prec (length fas) (slp_of_fas fas) a)"])
    apply (auto simp: option.splits aforms_rel_def br_def env_len_def Joints_imp_length_eq)
    apply (auto dest!: approx_slp_outer_lengthD)[]
    using length_slp_of_fas_le trans_le_add1 approx_slp_outer_lengthD apply blast
    using approx_slp_outer_plain by blast
  done

lemma fresh_index_aforms_Nil[simp]: "fresh_index_aforms [] = 0"
  by (auto simp: fresh_index_aforms_def)

lemma independent_aforms_Nil[simp]:
  "independent_aforms x [] = [x]"
  by (auto simp: independent_aforms_def)

lemma mem_Joints_zero_iff[simp]: "x # xs ∈ Joints ((0, zero_pdevs) # XS) ⟷ (x = 0 ∧ xs ∈ Joints XS)"
  by (auto simp: Joints_def valuate_def)

lemma Joints_independent_aforms_eq: "Joints (independent_aforms x xs) = set_Cons (Affine x) (Joints xs)"
  by (simp add: independent_aforms_def Joints_msum_aform degree_le_fresh_index set_Cons_def)

lemma independent_aforms_refine: "(independent_aforms, set_Cons) ∈ ⟨rnv_rel⟩aform_rel → aforms_rel → aforms_rel"
  by (auto simp: aforms_rel_def br_def aform_rel_def Joints_independent_aforms_eq)

end


locale aform_approximate_sets = approximate_sets
  aform_ops
  Joints
  aforms_rel
begin

lemma Joints_in_lv_rel_set_relD:
  "(Joints xs, X) ∈ ⟨lv_rel⟩set_rel ⟹ X = Affine (eucl_of_list_aform xs)"
  unfolding lv_rel_def set_rel_br
  by (auto simp: br_def Joints_imp_length_eq eucl_of_list_image_Joints[symmetric])

lemma ncc_precond: "ncc_precond TYPE('a::executable_euclidean_space)"
  unfolding ncc_precond_def ncc_def appr_rel_def
  by (auto simp: aforms_rel_def compact_Affine convex_Affine dest!: Joints_in_lv_rel_set_relD brD)

lemma fst_eucl_of_list_aform_map: "fst (eucl_of_list_aform (map (λx. (fst x, asdf x)) x)) =
  fst (eucl_of_list_aform x)"
  by (auto simp: eucl_of_list_aform_def o_def)

lemma
  Affine_pdevs_of_list:― ‹TODO: move!›
  "Affine (fst x, pdevs_of_list (map snd (list_of_pdevs (snd x)))) = Affine x"
  by (auto simp: Affine_def valuate_def aform_val_def
      elim: pdevs_val_of_list_of_pdevs2[where X = "snd x"]
        pdevs_val_of_list_of_pdevs[where X = "snd x"])

end

lemma aform_approximate_sets: "aform_approximate_sets prec"
  apply (unfold_locales)
  unfolding aform_ops_def approximate_set_ops.simps
  subgoal unfolding relAPP_def aforms_rel_def .
  subgoal by (force simp: aforms_of_ivls_refine)
  subgoal by (rule product_aforms_refine)
  subgoal by (rule msum_aforms'_refine)
  subgoal by (rule inf_aforms_refine)
  subgoal by (rule sup_aforms_refine)
  subgoal by (rule split_aform_largest_take_refine)
  subgoal by (rule aform_inf_inner_refine)
  subgoal by (rule aform_sup_inner_refine)
  subgoal by (rule inter_aform_plane_refine) simp_all
  subgoal by (auto split: option.splits intro!: reduce_aform_refine)
  subgoal by (force simp: width_spec_def nres_rel_def)
  subgoal by (rule approx_slp_refine)
  subgoal by (rule aform_euclarithform_refine)
  subgoal by (rule aform_isFDERIV)
  subgoal by simp
  subgoal by (auto simp: Joints_imp_length_eq)
  subgoal by (force simp: Affine_def Joints_def valuate_def intro!:)
  subgoal by (force simp: Affine_def Joints_def valuate_def intro!:)
  subgoal by (auto simp: Joints_imp_length_eq)
  done

end