Theory Affine_Arithmetic.Floatarith_Expression

section ‹Operations on Expressions›
theory Floatarith_Expression
imports
  "HOL-Decision_Procs.Approximation"
  Affine_Arithmetic_Auxiliarities
  Executable_Euclidean_Space
begin

text ‹Much of this could move to the distribution...›

subsection ‹Approximating Expression*s*›

unbundle floatarith_notation

text ‹\label{sec:affineexpr}›

primrec interpret_floatariths :: "floatarith list  real list  real list"
where
    "interpret_floatariths [] vs = []"
  | "interpret_floatariths (a#bs) vs = interpret_floatarith a vs#interpret_floatariths bs vs"

lemma length_interpret_floatariths[simp]: "length (interpret_floatariths fas xs) = length fas"
  by (induction fas) auto

lemma interpret_floatariths_nth[simp]:
  "interpret_floatariths fas xs ! n = interpret_floatarith (fas ! n) xs"
  if "n < length fas"
  using that
  by (induction fas arbitrary: n) (auto simp: nth_Cons split: nat.splits)

abbreviation "einterpret  λfas vs. eucl_of_list (interpret_floatariths fas vs)"

subsection ‹Syntax›

syntax interpret_floatarith::"floatarith  real list  real"

instantiation floatarith :: "{plus, minus, uminus, times, inverse, zero, one}"
begin

definition "- f = Minus f"
lemma interpret_floatarith_uminus[simp]:
  "interpret_floatarith (- f) xs = - interpret_floatarith f xs"
  by (auto simp: uminus_floatarith_def)

definition "f + g = Add f g"
lemma interpret_floatarith_plus[simp]:
  "interpret_floatarith (f + g) xs = interpret_floatarith f xs + interpret_floatarith g xs"
  by (auto simp: plus_floatarith_def)

definition "f - g = Add f (Minus g)"
lemma interpret_floatarith_minus[simp]:
  "interpret_floatarith (f - g) xs = interpret_floatarith f xs - interpret_floatarith g xs"
  by (auto simp: minus_floatarith_def)

definition "inverse f = Inverse f"
lemma interpret_floatarith_inverse[simp]:
  "interpret_floatarith (inverse f) xs = inverse (interpret_floatarith f xs)"
  by (auto simp: inverse_floatarith_def)

definition "f * g = Mult f g"
lemma interpret_floatarith_times[simp]:
  "interpret_floatarith (f * g) xs = interpret_floatarith f xs * interpret_floatarith g xs"
  by (auto simp: times_floatarith_def)

definition "f div g = f * Inverse g"
lemma interpret_floatarith_divide[simp]:
  "interpret_floatarith (f div g) xs = interpret_floatarith f xs / interpret_floatarith g xs"
  by (auto simp: divide_floatarith_def inverse_eq_divide)

definition "1 = Num 1"
lemma interpret_floatarith_one[simp]:
  "interpret_floatarith 1 xs = 1"
  by (auto simp: one_floatarith_def)

definition "0 = Num 0"
lemma interpret_floatarith_zero[simp]:
  "interpret_floatarith 0 xs = 0"
  by (auto simp: zero_floatarith_def)

instance proof qed
end


subsection ‹Derived symbols›

definition "Re r = (case quotient_of r of (n, d)  Num (of_int n) / Num (of_int d))"
declare [[coercion Re ]]

lemma interpret_Re[simp]: "interpret_floatarith (Re x) xs = real_of_rat x"
  by (auto simp: Re_def of_rat_divide dest!: quotient_of_div split: prod.splits)

definition "Sin x = Cos ((Pi * (Num (Float 1 (-1)))) - x)"

lemma interpret_floatarith_Sin[simp]:
  "interpret_floatarith (Sin x) vs = sin (interpret_floatarith x vs)"
  by (auto simp: Sin_def approximation_preproc_floatarith(11))

definition "Half x = Mult (Num (Float 1 (-1))) x"
lemma interpret_Half[simp]: "interpret_floatarith (Half x) xs = interpret_floatarith x xs / 2"
  by (auto simp: Half_def)

definition "Tan x = (Sin x) / (Cos x)"

lemma interpret_floatarith_Tan[simp]:
  "interpret_floatarith (Tan x) vs = tan (interpret_floatarith x vs)"
  by (auto simp: Tan_def approximation_preproc_floatarith(12) inverse_eq_divide)

primrec Sume where
  "Sume f [] = 0"
| "Sume f (x#xs) = f x + Sume f xs" 

lemma interpret_floatarith_Sume[simp]:
  "interpret_floatarith (Sume f x) vs = (ix. interpret_floatarith (f i) vs)"
  by (induction x) auto

definition Norm where "Norm is = Sqrt (Sume (λi. i * i) is)"

lemma interpret_floatarith_norm[simp]:
  assumes [simp]: "length x = DIM('a)"
  shows "interpret_floatarith (Norm x) vs = norm (einterpret x vs::'a::executable_euclidean_space)"
  apply (auto simp: Norm_def norm_eq_sqrt_inner)
  apply (subst euclidean_inner[where 'a='a])
  apply (auto simp: power2_eq_square[symmetric] )
  apply (subst sum_list_Basis_list[symmetric])
  apply (rule sum_list_nth_eqI)
  by (auto simp: in_set_zip eucl_of_list_inner)

notation floatarith.Power (infixr "^e" 80)

subsection ‹Constant Folding›

fun dest_Num_fa where
  "dest_Num_fa (floatarith.Num x) = Some x"
| "dest_Num_fa _ = None"

fun_cases dest_Num_fa_None: "dest_Num_fa fa = None"
  and dest_Num_fa_Some: "dest_Num_fa fa = Some x"

fun fold_const_fa where
  "fold_const_fa (Add fa1 fa2) =
    (let (ffa1, ffa2) = (fold_const_fa fa1, fold_const_fa fa2)
    in case (dest_Num_fa ffa1, dest_Num_fa (ffa2)) of
      (Some a, Some b)  Num (a + b)
    | (Some a, None)  (if a = 0 then ffa2 else Add (Num a) ffa2)
    | (None, Some a)  (if a = 0 then ffa1 else Add ffa1 (Num a))
    | (None, None)  Add ffa1 ffa2)"
| "fold_const_fa (Minus a) =
    (case (fold_const_fa a) of
      (Num x)  Num (-x)
    | x  Minus x)"
| "fold_const_fa (Mult fa1 fa2) =
    (let (ffa1, ffa2) = (fold_const_fa fa1, fold_const_fa fa2)
  in case (dest_Num_fa ffa1, dest_Num_fa (ffa2)) of
    (Some a, Some b)  Num (a * b)
  | (Some a, None)  (if a = 0 then Num 0 else if a = 1 then ffa2 else Mult (Num a) ffa2)
  | (None, Some a)  (if a = 0 then Num 0 else if a = 1 then ffa1 else Mult ffa1 (Num a))
  | (None, None)  Mult ffa1 ffa2)"
| "fold_const_fa (Inverse a) = Inverse (fold_const_fa a)"
| "fold_const_fa (Abs a) =
    (case (fold_const_fa a) of
      (Num x)  Num (abs x)
    | x  Abs x)"
| "fold_const_fa (Max a b) =
    (case (fold_const_fa a, fold_const_fa b) of
      (Num x, Num y)  Num (max x y)
    | (x, y)  Max x y)"
| "fold_const_fa (Min a b) =
    (case (fold_const_fa a, fold_const_fa b) of
      (Num x, Num y)  Num (min x y)
    | (x, y)  Min x y)"
| "fold_const_fa (Floor a) =
    (case (fold_const_fa a) of
      (Num x)  Num (floor_fl x)
    | x  Floor x)"
| "fold_const_fa (Power a b) =
    (case (fold_const_fa a) of
      (Num x)  Num (x ^ b)
    | x  Power x b)"
| "fold_const_fa (Cos a) = Cos (fold_const_fa a)"
| "fold_const_fa (Arctan a) = Arctan (fold_const_fa a)"
| "fold_const_fa (Sqrt a) = Sqrt (fold_const_fa a)"
| "fold_const_fa (Exp a) = Exp (fold_const_fa a)"
| "fold_const_fa (Ln a) = Ln (fold_const_fa a)"
| "fold_const_fa (Powr a b) = Powr (fold_const_fa a) (fold_const_fa b)"
| "fold_const_fa Pi = Pi"
| "fold_const_fa (Var v) = Var v"
| "fold_const_fa (Num x) = Num x"

fun_cases fold_const_fa_Num: "fold_const_fa fa = Num y"
  and fold_const_fa_Add: "fold_const_fa fa = Add x y"
  and fold_const_fa_Minus: "fold_const_fa fa = Minus y"

lemma fold_const_fa[simp]: "interpret_floatarith (fold_const_fa fa) xs = interpret_floatarith fa xs"
  by (induction fa) (auto split!: prod.splits floatarith.splits option.splits
      elim!: dest_Num_fa_None dest_Num_fa_Some
      simp: max_def min_def floor_fl_def)


subsection ‹Free Variables›

primrec max_Var_floatarith where― ‹TODO: include bound in predicate›
  "max_Var_floatarith (Add a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
| "max_Var_floatarith (Mult a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
| "max_Var_floatarith (Inverse a) = max_Var_floatarith a"
| "max_Var_floatarith (Minus a) = max_Var_floatarith a"
| "max_Var_floatarith (Num a) = 0"
| "max_Var_floatarith (Var i) = Suc i"
| "max_Var_floatarith (Cos a) = max_Var_floatarith a"
| "max_Var_floatarith (floatarith.Arctan a) = max_Var_floatarith a"
| "max_Var_floatarith (Abs a) = max_Var_floatarith a"
| "max_Var_floatarith (floatarith.Max a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
| "max_Var_floatarith (floatarith.Min a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
| "max_Var_floatarith (floatarith.Pi) = 0"
| "max_Var_floatarith (Sqrt a) = max_Var_floatarith a"
| "max_Var_floatarith (Exp a) = max_Var_floatarith a"
| "max_Var_floatarith (Powr a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
| "max_Var_floatarith (floatarith.Ln a) = max_Var_floatarith a"
| "max_Var_floatarith (Power a n) = max_Var_floatarith a"
| "max_Var_floatarith (Floor a) = max_Var_floatarith a"       
  
primrec max_Var_floatariths where
  "max_Var_floatariths [] = 0"
| "max_Var_floatariths (x#xs) = max (max_Var_floatarith x) (max_Var_floatariths xs)"

primrec max_Var_form where
  "max_Var_form (Conj a b) = max (max_Var_form a) (max_Var_form b)"
|  "max_Var_form (Disj a b) = max (max_Var_form a) (max_Var_form b)"
|  "max_Var_form (Less a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
|  "max_Var_form (LessEqual a b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
|  "max_Var_form (Bound a b c d) = linorder_class.Max {max_Var_floatarith a,max_Var_floatarith b, max_Var_floatarith c, max_Var_form d}"
|  "max_Var_form (AtLeastAtMost a b c) = linorder_class.Max {max_Var_floatarith a,max_Var_floatarith b, max_Var_floatarith c}"
|  "max_Var_form (Assign a b c) = linorder_class.Max {max_Var_floatarith a,max_Var_floatarith b, max_Var_form c}"

lemma
  interpret_floatarith_eq_take_max_VarI:
  assumes "take (max_Var_floatarith ra) ys = take (max_Var_floatarith ra) zs"
  shows "interpret_floatarith ra ys = interpret_floatarith ra zs"
  using assms
  by (induct ra) (auto dest!: take_max_eqD simp: take_Suc_eq split: if_split_asm)

lemma
  interpret_floatariths_eq_take_max_VarI:
  assumes "take (max_Var_floatariths ea) ys = take (max_Var_floatariths ea) zs"
  shows "interpret_floatariths ea ys = interpret_floatariths ea zs"
  using assms
  apply (induction ea)
  subgoal by simp
  subgoal by (clarsimp) (metis interpret_floatarith_eq_take_max_VarI take_map take_max_eqD)
  done


lemma Max_Image_distrib:
  includes no_floatarith_notation
  assumes "finite X" "X  {}"
  shows "Max ((λx. max (f1 x) (f2 x)) ` X) = max (Max (f1 ` X)) (Max (f2 ` X))"
  apply (rule Max_eqI)
  subgoal using assms by simp
  subgoal for y
    using assms
    by (force intro: max.coboundedI1 max.coboundedI2 Max_ge)
  subgoal
  proof -
    have "Max (f1 ` X)  f1 ` X" using assms by auto
    then obtain x1 where x1: "x1  X" "Max (f1 ` X) = f1 x1" by auto
    have "Max (f2 ` X)  f2 ` X" using assms by auto
    then obtain x2 where x2: "x2  X" "Max (f2 ` X) = f2 x2" by auto
    show ?thesis
      apply (rule image_eqI[where x="if f1 x1  f2 x2 then x2 else x1"])
      using x1 x2 assms
       apply (auto simp: max_def)
       apply (metis Max_ge dual_order.trans finite_imageI image_eqI assms(1))
      apply (metis Max_ge dual_order.trans finite_imageI image_eqI assms(1))
      done
  qed
  done

lemma max_Var_floatarith_simps[simp]:
  "max_Var_floatarith (a / b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
  "max_Var_floatarith (a * b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
  "max_Var_floatarith (a + b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
  "max_Var_floatarith (a - b) = max (max_Var_floatarith a) (max_Var_floatarith b)"
  "max_Var_floatarith (- b) = (max_Var_floatarith b)"
  by (auto simp: divide_floatarith_def times_floatarith_def plus_floatarith_def minus_floatarith_def
      uminus_floatarith_def)

lemma max_Var_floatariths_Max:
  "max_Var_floatariths xs = (if set xs = {} then 0 else linorder_class.Max (max_Var_floatarith ` set xs))"
  by (induct xs) auto


lemma max_Var_floatariths_map_plus[simp]:
  "max_Var_floatariths (map (λi. fa1 i + fa2 i) xs) = max (max_Var_floatariths (map fa1 xs)) (max_Var_floatariths (map fa2 xs))"
  by (auto simp: max_Var_floatariths_Max image_image Max_Image_distrib)

lemma max_Var_floatariths_map_times[simp]:
  "max_Var_floatariths (map (λi. fa1 i * fa2 i) xs) = max (max_Var_floatariths (map fa1 xs)) (max_Var_floatariths (map fa2 xs))"
  by (auto simp: max_Var_floatariths_Max image_image Max_Image_distrib)

lemma max_Var_floatariths_map_divide[simp]:
  "max_Var_floatariths (map (λi. fa1 i / fa2 i) xs) = max (max_Var_floatariths (map fa1 xs)) (max_Var_floatariths (map fa2 xs))"
  by (auto simp: max_Var_floatariths_Max image_image Max_Image_distrib)

lemma max_Var_floatariths_map_uminus[simp]:
  "max_Var_floatariths (map (λi. - fa1 i) xs) = max_Var_floatariths (map fa1 xs)"
  by (auto simp: max_Var_floatariths_Max image_image Max_Image_distrib)

lemma max_Var_floatariths_map_const[simp]:
  "max_Var_floatariths (map (λi. fa) xs) = (if xs = [] then 0 else max_Var_floatarith fa)"
  by (auto simp: max_Var_floatariths_Max image_image image_constant_conv)

lemma max_Var_floatariths_map_minus[simp]:
  "max_Var_floatariths (map (λi. fa1 i - fa2 i) xs) = max (max_Var_floatariths (map fa1 xs)) (max_Var_floatariths (map fa2 xs))"
  by (auto simp: max_Var_floatariths_Max image_image Max_Image_distrib)


primrec fresh_floatarith where
  "fresh_floatarith (Var y) x  (x  y)"
| "fresh_floatarith (Num a) x  True"
| "fresh_floatarith Pi x  True"
| "fresh_floatarith (Cos a) x  fresh_floatarith a x"
| "fresh_floatarith (Abs a) x  fresh_floatarith a x"
| "fresh_floatarith (Arctan a) x  fresh_floatarith a x"
| "fresh_floatarith (Sqrt a) x  fresh_floatarith a x"
| "fresh_floatarith (Exp a) x  fresh_floatarith a x"
| "fresh_floatarith (Floor a) x  fresh_floatarith a x"
| "fresh_floatarith (Power a n) x  fresh_floatarith a x"
| "fresh_floatarith (Minus a) x  fresh_floatarith a x"
| "fresh_floatarith (Ln a) x  fresh_floatarith a x"
| "fresh_floatarith (Inverse a) x  fresh_floatarith a x"
| "fresh_floatarith (Add a b) x  fresh_floatarith a x  fresh_floatarith b x"
| "fresh_floatarith (Mult a b) x  fresh_floatarith a x  fresh_floatarith b x"
| "fresh_floatarith (Max a b) x  fresh_floatarith a x  fresh_floatarith b x"
| "fresh_floatarith (Min a b) x  fresh_floatarith a x  fresh_floatarith b x"
| "fresh_floatarith (Powr a b) x  fresh_floatarith a x  fresh_floatarith b x"

lemma fresh_floatarith_subst:
  fixes v::float
  assumes "fresh_floatarith e x"
  assumes "x < length vs"
  shows "interpret_floatarith e (vs[x:=v]) = interpret_floatarith e vs"
  using assms
  by (induction e) (auto simp: map_update)

lemma fresh_floatarith_max_Var:
  assumes "max_Var_floatarith ea  i"
  shows "fresh_floatarith ea i"
  using assms
  by (induction ea) auto

primrec fresh_floatariths where
  "fresh_floatariths [] x  True"
| "fresh_floatariths (a#as) x  fresh_floatarith a x  fresh_floatariths as x"

lemma fresh_floatariths_max_Var:
  assumes "max_Var_floatariths ea  i"
  shows "fresh_floatariths ea i"
  using assms
  by (induction ea) (auto simp: fresh_floatarith_max_Var)

lemma
  interpret_floatariths_take_eqI:
  assumes "take n ys = take n zs"
  assumes "max_Var_floatariths ea  n"
  shows "interpret_floatariths ea ys = interpret_floatariths ea zs"
  by (rule interpret_floatariths_eq_take_max_VarI) (rule take_greater_eqI[OF assms])

lemma
  interpret_floatarith_fresh_eqI:
  assumes "i. fresh_floatarith ea i  (i < length ys  i < length zs  ys ! i = zs ! i)"
  shows "interpret_floatarith ea ys = interpret_floatarith ea zs"
  using assms
  by (induction ea) force+

lemma
  interpret_floatariths_fresh_eqI:
  assumes "i. fresh_floatariths ea i  (i < length ys  i < length zs  ys ! i = zs ! i)"
  shows "interpret_floatariths ea ys = interpret_floatariths ea zs"
  using assms
  apply (induction ea)
  subgoal by (force simp: interpret_floatarith_fresh_eqI intro: interpret_floatarith_fresh_eqI)
  subgoal for e ea
    apply clarsimp
    apply (auto simp: list_eq_iff_nth_eq)
    using interpret_floatarith_fresh_eqI by blast
  done

lemma
  interpret_floatarith_max_Var_cong:
  assumes "i. i < max_Var_floatarith f  xs ! i = ys ! i"
  shows "interpret_floatarith f ys = interpret_floatarith f xs"
  using assms
  by (induction f) auto

lemma
  interpret_floatarith_fresh_cong:
  assumes "i. ¬fresh_floatarith f i  xs ! i = ys ! i"
  shows "interpret_floatarith f ys = interpret_floatarith f xs"
  using assms
  by (induction f) auto

lemma max_Var_floatarith_le_max_Var_floatariths:
  "fa  set fas  max_Var_floatarith fa  max_Var_floatariths fas"
  by (induction fas) (auto simp: nth_Cons max_def split: nat.splits)

lemma max_Var_floatarith_le_max_Var_floatariths_nth:
  "n < length fas  max_Var_floatarith (fas ! n)  max_Var_floatariths fas"
  by (rule max_Var_floatarith_le_max_Var_floatariths) auto

lemma max_Var_floatariths_leI:
  assumes "i. i < length xs  max_Var_floatarith (xs ! i)  F"
  shows "max_Var_floatariths xs  F"
  using assms
  by (auto simp: max_Var_floatariths_Max in_set_conv_nth)

lemma fresh_floatariths_map_Var[simp]:
  "fresh_floatariths (map floatarith.Var xs) i  i  set xs"
  by (induction xs) auto


lemma max_Var_floatarith_fold_const_fa:
  "max_Var_floatarith (fold_const_fa fa)  max_Var_floatarith fa"
  by (induction fa) (auto simp: fold_const_fa.simps split!: option.splits floatarith.splits)

lemma max_Var_floatariths_fold_const_fa:
  "max_Var_floatariths (map fold_const_fa xs)  max_Var_floatariths xs"
  by (auto simp: intro!: max_Var_floatariths_leI max_Var_floatarith_le_max_Var_floatariths_nth
      max_Var_floatarith_fold_const_fa[THEN order_trans])

lemma interpret_form_max_Var_cong:
  assumes "i. i < max_Var_form f  xs ! i = ys ! i"
  shows "interpret_form f xs = interpret_form f ys"
  using assms
  by (induction f) (auto simp: interpret_floatarith_max_Var_cong[where xs=xs and ys=ys])

lemma max_Var_floatariths_lessI: "i < max_Var_floatarith (fas ! j)  j < length fas  i < max_Var_floatariths fas"
  by (metis leD le_trans less_le max_Var_floatarith_le_max_Var_floatariths nth_mem)

lemma interpret_floatariths_max_Var_cong:
  assumes "i. i < max_Var_floatariths f  xs ! i = ys ! i"
  shows "interpret_floatariths f ys = interpret_floatariths f xs"
  by (auto intro!: nth_equalityI interpret_floatarith_max_Var_cong assms max_Var_floatariths_lessI)


lemma max_Var_floatarithimage_Var[simp]: "max_Var_floatarith ` Var ` X = Suc ` X" by force

lemma max_Var_floatariths_map_Var[simp]:
  "max_Var_floatariths (map Var xs) = (if xs = [] then 0 else Suc (linorder_class.Max (set xs)))"
  by (auto simp: max_Var_floatariths_Max hom_Max_commute split: if_splits)

lemma Max_atLeastLessThan_nat[simp]: "a < b  linorder_class.Max {a..<b} = b - 1" for a b::nat
  by (auto intro!: Max_eqI)


subsection ‹Derivatives›

lemma isDERIV_Power_iff: "isDERIV j (Power fa n) xs = (if n = 0 then True else isDERIV j fa xs)"
  by (cases n) auto

lemma isDERIV_max_Var_floatarithI:
  assumes "isDERIV n f ys"
  assumes "i. i < max_Var_floatarith f  xs ! i = ys ! i"
  shows "isDERIV n f xs"
  using assms
proof (induction f)
  case (Power f n) then show ?case by (cases n) auto
qed (auto simp: max_def interpret_floatarith_max_Var_cong[of _ xs ys] split: if_splits)

definition isFDERIV where "isFDERIV n xs fas vs 
  (i<n. j<n. isDERIV (xs ! i) (fas ! j) vs)  length fas = n  length xs = n"

lemma isFDERIV_I: "(i j. i < n  j < n  isDERIV (xs ! i) (fas ! j) vs) 
  length fas = n  length xs = n  isFDERIV n xs fas vs"
  by (auto simp: isFDERIV_def)

lemma isFDERIV_isDERIV_D: "isFDERIV n xs fas vs  i < n  j < n  isDERIV (xs ! i) (fas ! j) vs"
  by (auto simp: isFDERIV_def)

lemma isFDERIV_lengthD: "length fas = n" "length xs = n" if "isFDERIV n xs fas vs"
  using that by (auto simp: isFDERIV_def)

lemma isFDERIV_uptD: "isFDERIV n [0..<n] fas vs  i < n  j < n  isDERIV i (fas ! j) vs"
  by (auto simp: isFDERIV_def)

lemma isFDERIV_max_Var_congI: "isFDERIV n xs fas ws"
  if f: "isFDERIV n xs fas vs" and c: "(i. i < max_Var_floatariths fas  vs ! i = ws ! i)"
  using c f
  by (auto simp: isFDERIV_def max_Var_floatariths_lessI
      intro!: isFDERIV_I isDERIV_max_Var_floatarithI[OF isFDERIV_isDERIV_D[OF f]])

lemma isFDERIV_max_Var_cong: "isFDERIV n xs fas ws  isFDERIV n xs fas vs"
  if c: "(i. i < max_Var_floatariths fas  vs ! i = ws ! i)"
  using c by (auto intro: isFDERIV_max_Var_congI)

lemma isDERIV_max_VarI:
  "i  max_Var_floatarith fa  isDERIV j fa xs  isDERIV i fa xs"
  by (induction fa) (auto simp: isDERIV_Power_iff)

lemmas max_Var_floatarith_le_max_Var_floatariths_nthI =
  max_Var_floatarith_le_max_Var_floatariths_nth[THEN order_trans]


lemma
  isFDERIV_appendD1:
  assumes "isFDERIV (J + K) [0..<J + K] (es @ rs) xs"
  assumes "length es = J"
  assumes "length rs = K"
  assumes "max_Var_floatariths es  J"
  shows "isFDERIV J [0..<J] (es) xs"
  unfolding isFDERIV_def
  apply (safe)
  subgoal for i j
    using assms
    apply (cases "i < length es")
    subgoal by (auto simp: nth_append isFDERIV_def) (metis add.commute trans_less_add2)
    subgoal
      apply (rule isDERIV_max_VarI[where j=0])
       apply (rule max_Var_floatarith_le_max_Var_floatariths_nthI)
         apply force
        apply force
       apply force
      done
    done
  subgoal by (auto simp: assms)
  subgoal by (auto simp: assms)
  done

lemma interpret_floatariths_Var[simp]:
  "interpret_floatariths (map floatarith.Var xs) vs = (map (nth vs) xs)"
  by (induction xs) auto

lemma max_Var_floatariths_append[simp]: "max_Var_floatariths (xs @ ys) = max (max_Var_floatariths xs) (max_Var_floatariths ys)"
  by (induction xs) (auto)

lemma map_nth_append_upt[simp]:
  assumes "a  length xs"
  shows "map ((!) (xs @ ys)) [a..<b] = map ((!) ys) [a - length xs..<b - length xs]"
  using assms
  by (auto intro!: nth_equalityI simp: nth_append)

lemma map_nth_Cons_upt[simp]:
  assumes "a > 0"
  shows "map ((!) (x # ys)) [a..<b] = map ((!) ys) [a - Suc 0..<b - Suc 0]"
  using assms
  by (auto intro!: nth_equalityI simp: nth_append)

lemma map_nth_eq_self[simp]:
  shows "length fas = l  (map ((!) fas) [0..<l]) = fas"
  by (auto simp: intro!: nth_equalityI)


lemma
  isFDERIV_appendI1:
  assumes "isFDERIV J [0..<J] (es) xs"
  assumes "i j. i < J + K  j < K  isDERIV i (rs ! j) xs"
  assumes "length es = J"
  assumes "length rs = K"
  assumes "max_Var_floatariths es  J"
  shows "isFDERIV (J + K) [0..<J + K] (es @ rs) xs"
  unfolding isFDERIV_def
  apply safe
  subgoal for i j
    using assms
    apply (cases "j < length es")
    subgoal
      apply (auto simp: nth_append isFDERIV_def)
      by (metis (no_types, opaque_lifting) isDERIV_max_VarI le_trans less_le
          max_Var_floatarith_le_max_Var_floatariths_nthI nat_le_linear)
    subgoal by (auto simp: nth_append)
    done
  subgoal by (auto simp: assms)
  subgoal by (auto simp: assms)
  done


lemma matrix_matrix_mult_zero[simp]:
  "a ** 0 = 0" "0 ** a = 0"
  by (vector matrix_matrix_mult_def)+

lemma scaleR_blinfun_compose_left: "i *R (A oL B) = i *R A oL B"
  and scaleR_blinfun_compose_right: "i *R (A oL B) = A oL i *R B"
  by (auto intro!: blinfun_eqI simp: blinfun.bilinear_simps)

lemma
  matrix_blinfun_compose:
  fixes A B::"(real ^ 'n) L (real ^ 'n)"
  shows "matrix (A oL B) = (matrix A) ** (matrix B)"
  by transfer (auto simp: matrix_compose linear_linear)

lemma matrix_add_rdistrib: "((B + C) ** A) = (B ** A) + (C ** A)"
  by (vector matrix_matrix_mult_def sum.distrib[symmetric] field_simps)

lemma matrix_scaleR_right: "r *R (a::'a::real_algebra_1^'n^'m) ** b = r *R (a ** b)"
  by (vector matrix_matrix_mult_def algebra_simps scaleR_sum_right)

lemma matrix_scaleR_left: "(a::'a::real_algebra_1^'n^'m) ** r *R b = r *R (a ** b)"
  by (vector matrix_matrix_mult_def algebra_simps scaleR_sum_right)

lemma bounded_bilinear_matrix_matrix_mult[bounded_bilinear]:
   "bounded_bilinear ((**)::
    ('a::{euclidean_space, real_normed_algebra_1}^'n^'m) 
    ('a::{euclidean_space, real_normed_algebra_1}^'p^'n) 
    ('a::{euclidean_space, real_normed_algebra_1}^'p^'m))"
  unfolding bilinear_conv_bounded_bilinear[symmetric]
  unfolding bilinear_def
  apply safe
  by unfold_locales (auto simp: matrix_add_ldistrib matrix_add_rdistrib matrix_scaleR_right matrix_scaleR_left)

lemma norm_axis: "norm (axis ia 1::'a::{real_normed_algebra_1}^'n) = 1"
  by (auto simp: axis_def norm_vec_def L2_set_def if_distrib if_distribR sum.delta
      cong: if_cong)

lemma abs_vec_nth_blinfun_apply_lemma:
  fixes x::"(real^'n) L (real^'m)"
  shows "abs (vec_nth (blinfun_apply x (axis ia 1)) i)  norm x"
  apply (rule component_le_norm_cart[THEN order_trans])
  apply (rule norm_blinfun[THEN order_trans])
  by (auto simp: norm_axis)

lemma bounded_linear_matrix_blinfun_apply: "bounded_linear (λx::(real^'n) L (real^'m). matrix (blinfun_apply x))"
  apply standard
  subgoal by (vector blinfun.bilinear_simps matrix_def)
  subgoal by (vector blinfun.bilinear_simps matrix_def)
  apply (rule exI[where x="real (CARD('n) * CARD('m))"])
  apply (auto simp: matrix_def)
  apply (subst norm_vec_def)
  apply (rule L2_set_le_sum[THEN order_trans])
  apply simp
  apply auto
  apply (rule sum_mono[THEN order_trans])
  apply (subst norm_vec_def)
   apply (rule L2_set_le_sum)
   apply simp
  apply (rule sum_mono[THEN order_trans])
   apply (rule sum_mono)
    apply simp
    apply (rule abs_vec_nth_blinfun_apply_lemma)
  apply (simp add: abs_vec_nth_blinfun_apply_lemma)
  done

lemma matrix_has_derivative:
  shows "((λx::(real^'n)L(real^'n). matrix (blinfun_apply x)) has_derivative (λh. matrix (blinfun_apply h))) (at x)"
  apply (auto simp: has_derivative_at2)
  unfolding linear_linear
  subgoal by (rule bounded_linear_matrix_blinfun_apply)
  subgoal
    by (auto simp: blinfun.bilinear_simps matrix_def) vector
  done

lemma
  matrix_comp_has_derivative[derivative_intros]:
  fixes f::"'a::real_normed_vector  ((real^'n)L(real^'n))"
  assumes "(f has_derivative f') (at x within S)"
  shows "((λx. matrix (blinfun_apply (f x))) has_derivative (λx. matrix (f' x))) (at x within S)"
  using has_derivative_compose[OF assms matrix_has_derivative]
  by auto

fun inner_floatariths where
  "inner_floatariths [] _ = Num 0"
| "inner_floatariths _ [] = Num 0"
| "inner_floatariths (x#xs) (y#ys) = Add (Mult x y) (inner_floatariths xs ys)"

lemma interpret_floatarith_inner_eq:
  assumes "length xs = length ys"
  shows "interpret_floatarith (inner_floatariths xs ys) vs =
    (i<length ys. (interpret_floatariths xs vs ! i) * (interpret_floatariths ys vs ! i))"
  using assms
proof (induction rule: list_induct2)
  case Nil
  then show ?case by simp
next
  case (Cons x xs y ys)
  then show ?case
    unfolding length_Cons sum.lessThan_Suc_shift
    by simp
qed

lemma
  interpret_floatarith_inner_floatariths:
  assumes "length xs = DIM('a::executable_euclidean_space)"
  assumes "length ys = DIM('a)"
  assumes "eucl_of_list (interpret_floatariths xs vs) = (x::'a)"
  assumes "eucl_of_list (interpret_floatariths ys vs) = y"
  shows "interpret_floatarith (inner_floatariths xs ys) vs = x  y"
  using assms
  by (subst euclidean_inner)
    (auto simp: interpret_floatarith_inner_eq sum_Basis_sum_nth_Basis_list eucl_of_list_inner
      index_nth_id
      intro!: euclidean_eqI[where 'a='a] sum.cong)

lemma max_Var_floatarith_inner_floatariths[simp]:
  assumes "length f = length g"
  shows "max_Var_floatarith (inner_floatariths f g) = max (max_Var_floatariths f) (max_Var_floatariths g)"
  using assms
  by (induction f g rule: list_induct2) auto


definition FDERIV_floatarith where
  "FDERIV_floatarith fa xs d = inner_floatariths (map (λx. fold_const_fa (DERIV_floatarith x fa)) xs) d"
― ‹TODO: specialize to FDERIV_floatarith fa [0..<n] [m..<m + n]› and do the rest with @{term subst_floatarith}?
   TODO: introduce approximation on type @{typ "real^'i^'j"} and use @{term jacobian}?›

lemma interpret_floatariths_map: "interpret_floatariths (map f xs) vs = map (λx. interpret_floatarith (f x) vs) xs"
  by (induct xs) auto

lemma max_Var_floatarith_DERIV_floatarith:
  "max_Var_floatarith (DERIV_floatarith x fa)  max_Var_floatarith fa"
  by (induction x fa rule: DERIV_floatarith.induct) (auto)

lemma max_Var_floatarith_FDERIV_floatarith:
  "length xs = length d 
    max_Var_floatarith (FDERIV_floatarith fa xs d)  max (max_Var_floatarith fa) (max_Var_floatariths d)"
  unfolding FDERIV_floatarith_def
  by (auto simp: max_Var_floatariths_Max intro!: max_Var_floatarith_DERIV_floatarith[THEN order_trans]
      max_Var_floatarith_fold_const_fa[THEN order_trans])

definition FDERIV_floatariths where
  "FDERIV_floatariths fas xs d = map (λfa. FDERIV_floatarith fa xs d) fas"

lemma max_Var_floatarith_FDERIV_floatariths:
  "length xs = length d  max_Var_floatariths (FDERIV_floatariths fa xs d)  max (max_Var_floatariths fa) (max_Var_floatariths d)"
  by (auto simp: FDERIV_floatariths_def max_Var_floatariths_Max
      intro!: max_Var_floatarith_FDERIV_floatarith[THEN order_trans])
    (auto simp: max_def)

lemma length_FDERIV_floatariths[simp]:
  "length (FDERIV_floatariths fas xs ds) = length fas"
  by (auto simp: FDERIV_floatariths_def)

lemma FDERIV_floatariths_nth[simp]:
  "i < length fas  FDERIV_floatariths fas xs ds ! i  = FDERIV_floatarith (fas ! i) xs ds"
  by (auto simp: FDERIV_floatariths_def)

definition "FDERIV_n_floatariths fas xs ds n = ((λx. FDERIV_floatariths x xs ds)^^n) fas"

lemma FDERIV_n_floatariths_Suc[simp]:
  "FDERIV_n_floatariths fa xs ds 0 = fa"
  "FDERIV_n_floatariths fa xs ds (Suc n) = FDERIV_floatariths (FDERIV_n_floatariths fa xs ds n) xs ds"
  by (auto simp: FDERIV_n_floatariths_def)

lemma length_FDERIV_n_floatariths[simp]: "length (FDERIV_n_floatariths fa xs ds n) = length fa"
  by (induction n) (auto simp: FDERIV_n_floatariths_def)

lemma max_Var_floatarith_FDERIV_n_floatariths:
  "length xs = length d  max_Var_floatariths (FDERIV_n_floatariths fa xs d n)  max (max_Var_floatariths fa) (max_Var_floatariths d)"
  by (induction n)
    (auto intro!: max_Var_floatarith_FDERIV_floatariths[THEN order_trans] simp: FDERIV_n_floatariths_def)

lemma interpret_floatarith_FDERIV_floatarith_cong:
  assumes rq: "i. i < max_Var_floatarith f  rs ! i = qs ! i"
  assumes [simp]: "length ds = length xs" "length es = length xs"
  assumes "interpret_floatariths ds qs = interpret_floatariths es rs"
  shows "interpret_floatarith (FDERIV_floatarith f xs ds) qs =
   interpret_floatarith (FDERIV_floatarith f xs es) rs"
  apply (auto simp: FDERIV_floatarith_def interpret_floatarith_inner_eq)
  apply (rule sum.cong[OF refl])
  subgoal premises prems for i
  proof -
    have "interpret_floatarith (DERIV_floatarith (xs ! i) f) qs = interpret_floatarith (DERIV_floatarith (xs ! i) f) rs"
      apply (rule interpret_floatarith_max_Var_cong)
      apply (auto simp: intro!: rq)
      by (metis leD le_trans max_Var_floatarith_DERIV_floatarith nat_less_le)
    moreover
    have "interpret_floatarith (ds ! i) qs = interpret_floatarith (es ! i) rs"
      using assms
      by (metis i  {..<length xs} interpret_floatariths_nth lessThan_iff)
    ultimately show ?thesis by auto
  qed
  done

theorem
  matrix_vector_mult_eq_list_of_eucl_nth:
  "(M::real^'n::enum^'m::enum) *v v =
    (i<CARD('m).
      (j<CARD('n). list_of_eucl M ! (i * CARD('n) + j) * list_of_eucl v ! j) *R Basis_list ! i)"
  using eucl_of_list_matrix_vector_mult_eq_sum_nth_Basis_list[of "list_of_eucl M" "list_of_eucl v",
      where 'n='n and 'm = 'm]
  by auto

definition "mmult_fa l m n AS BS =
  concat (map (λi. map (λk. inner_floatariths
    (map (λj. AS ! (i * m + j)) [0..<m]) (map (λj. BS ! (j * n + k)) [0..<m])) [0..<n]) [0..<l])"

lemma length_mmult_fa[simp]: "length (mmult_fa l m n AS BS) = l * n"
  by (auto simp: mmult_fa_def length_concat o_def sum_list_distinct_conv_sum_set)

lemma einterpret_mmult_fa:
  assumes [simp]: "Dn = CARD('n::enum)" "Dm = CARD('m::enum)" "Dl = CARD('l::enum)"
    "length A = CARD('l)*CARD('m)" "length B = CARD('m)*CARD('n)"
  shows "einterpret (mmult_fa Dl Dm Dn A B) vs = (einterpret A vs::((real, 'm::enum) vec, 'l) vec) ** (einterpret B vs::((real, 'n::enum) vec, 'm) vec)"
  apply (vector matrix_matrix_mult_def)
  apply (auto simp: mmult_fa_def vec_nth_eucl_of_list_eq2 index_Basis_list_axis2
      concat_map_map_index length_concat o_def sum_list_distinct_conv_sum_set
      interpret_floatarith_inner_eq)
  apply (subst sum_index_enum_eq)
  apply simp
  done

lemma max_Var_floatariths_mmult_fa:
  assumes [simp]: "length A = D * E" "length B = E * F"
  shows "max_Var_floatariths (mmult_fa D E F A B)  max (max_Var_floatariths A) (max_Var_floatariths B)"
  apply (auto simp: mmult_fa_def concat_map_map_index intro!: max_Var_floatariths_leI)
   apply (rule max.coboundedI1)
   apply (auto intro!: max_Var_floatarith_le_max_Var_floatariths_nth max.coboundedI2)
  apply (cases "F = 0")
   apply simp_all
  done

lemma isDERIV_inner_iff:
  assumes "length xs = length ys"
  shows "isDERIV i (inner_floatariths xs ys) vs 
    (k < length xs. isDERIV i (xs ! k) vs)  (k < length ys. isDERIV i (ys ! k) vs)"
  using assms
  by (induction xs ys rule: list_induct2) (auto simp: nth_Cons split: nat.splits)

lemma isDERIV_Power: "isDERIV x (fa) vs  isDERIV x (fa ^e n) vs"
  by (induction n) (auto simp: FDERIV_floatarith_def isDERIV_inner_iff)

lemma isDERIV_mmult_fa_nth:
  assumes "j. j < D * E  isDERIV i (A ! j) xs"
  assumes "j. j < E * F  isDERIV i (B ! j) xs"
  assumes [simp]: "length A = D * E" "length B = E * F" "j < D * F"
  shows "isDERIV i (mmult_fa D E F A B ! j) xs"
  using assms
  apply (cases "F = 0")
  apply (auto simp: mmult_fa_def concat_map_map_index isDERIV_inner_iff ac_simps)
  apply (metis add.commute assms(5) in_square_lemma less_square_imp_div_less mult.commute)
  done

definition "mvmult_fa n m AS B =
  map (λi. inner_floatariths (map (λj. AS ! (i * m + j)) [0..<m]) (map (λj. B ! j) [0..<m])) [0..<n]"

lemma einterpret_mvmult_fa:
  assumes [simp]: "Dn = CARD('n::enum)" "Dm = CARD('m::enum)"
    "length A = CARD('n)*CARD('m)" "length B = CARD('m)"
  shows "einterpret (mvmult_fa Dn Dm A B) vs = (einterpret A vs::((real, 'm::enum) vec, 'n) vec) *v (einterpret B vs::(real, 'm) vec)"
  apply (vector matrix_vector_mult_def)
  apply (auto simp: mvmult_fa_def vec_nth_eucl_of_list_eq2 index_Basis_list_axis2 index_Basis_list_axis1 vec_nth_eucl_of_list_eq
      concat_map_map_index length_concat o_def sum_list_distinct_conv_sum_set
      interpret_floatarith_inner_eq)
  apply (subst sum_index_enum_eq)
  apply simp
  done


lemma max_Var_floatariths_mvult_fa:
  assumes [simp]: "length A = D * E" "length B = E"
  shows "max_Var_floatariths (mvmult_fa D E A B)  max (max_Var_floatariths A) (max_Var_floatariths B)"
  apply (auto simp: mvmult_fa_def concat_map_map_index intro!: max_Var_floatariths_leI)
   apply (rule max.coboundedI1)
  by (auto intro!: max_Var_floatarith_le_max_Var_floatariths_nth max.coboundedI2)

lemma isDERIV_mvmult_fa_nth:
  assumes "j. j < D * E  isDERIV i (A ! j) xs"
  assumes "j. j < E  isDERIV i (B ! j) xs"
  assumes [simp]: "length A = D * E" "length B = E" "j < D"
  shows "isDERIV i (mvmult_fa D E A B ! j) xs"
  using assms
  apply (auto simp: mvmult_fa_def concat_map_map_index isDERIV_inner_iff ac_simps)
  by (metis assms(5) in_square_lemma semiring_normalization_rules(24) semiring_normalization_rules(7))

lemma max_Var_floatariths_mapI:
  assumes "x. x  set xs  max_Var_floatarith (f x)  m"
  shows "max_Var_floatariths (map f xs)  m"
  using assms
  by (force intro!: max_Var_floatariths_leI simp: in_set_conv_nth)

lemma
  max_Var_floatariths_list_updateI:
  assumes "max_Var_floatariths xs  m"
  assumes "max_Var_floatarith v  m"
  assumes "i < length xs"
  shows "max_Var_floatariths (xs[i := v])  m"
  using assms
  apply (auto simp: nth_list_update intro!: max_Var_floatariths_leI )
  using max_Var_floatarith_le_max_Var_floatariths_nthI by blast

lemma
  max_Var_floatariths_replicateI:
  assumes "max_Var_floatarith v  m"
  shows "max_Var_floatariths (replicate n v)  m"
  using assms
  by (auto intro!: max_Var_floatariths_leI )

definition "FDERIV_n_floatarith fa xs ds n = ((λx. FDERIV_floatarith x xs ds)^^n) fa"
lemma FDERIV_n_floatariths_nth: "i < length fas  FDERIV_n_floatariths fas xs ds n ! i = FDERIV_n_floatarith (fas ! i) xs ds n"
  by (induction n)
    (auto simp: FDERIV_n_floatarith_def FDERIV_floatariths_nth)


lemma einterpret_fold_const_fa[simp]:
  "(einterpret (map (λi. fold_const_fa (fa i)) xs) vs::'a::executable_euclidean_space) =
    einterpret (map fa xs) vs" if "length xs = DIM('a)"
  using that
  by (auto intro!: euclidean_eqI[where 'a='a] simp: algebra_simps eucl_of_list_inner)

lemma einterpret_plus[simp]:
  shows "(einterpret (map (λi. fa1 i + fa2 i) [0..<DIM('a)]) vs::'a) =
    einterpret (map fa1 [0..<DIM('a::executable_euclidean_space)]) vs + einterpret (map fa2 [0..<DIM('a)]) vs"
  by (auto intro!: euclidean_eqI[where 'a='a] simp: algebra_simps eucl_of_list_inner)

lemma einterpret_uminus[simp]:
  shows "(einterpret (map (λi. - fa1 i) [0..<DIM('a)]) vs::'a::executable_euclidean_space) =
    - einterpret (map fa1 [0..<DIM('a)]) vs"
  by (auto intro!: euclidean_eqI[where 'a='a] simp: algebra_simps eucl_of_list_inner)

lemma diff_floatarith_conv_add_uminus: "a - b = a + - b" for a b::floatarith
  by (auto simp: minus_floatarith_def plus_floatarith_def uminus_floatarith_def)

lemma einterpret_minus[simp]:
  shows "(einterpret (map (λi. fa1 i - fa2 i) [0..<DIM('a)]) vs::'a::executable_euclidean_space) =
    einterpret (map fa1 [0..<DIM('a)]) vs - einterpret (map fa2 [0..<DIM('a)]) vs"
  by (simp add: diff_floatarith_conv_add_uminus)

lemma einterpret_scaleR[simp]:
  shows "(einterpret (map (λi. fa1 * fa2 i) [0..<DIM('a)]) vs::'a::executable_euclidean_space) =
    interpret_floatarith (fa1) vs *R einterpret (map fa2 [0..<DIM('a)]) vs"
  by (auto intro!: euclidean_eqI[where 'a='a] simp: algebra_simps eucl_of_list_inner)

lemma einterpret_nth[simp]:
  assumes [simp]: "length xs = DIM('a)"
  shows "(einterpret (map ((!) xs) [0..<DIM('a)]) vs::'a::executable_euclidean_space) = einterpret xs vs"
  by (auto intro!: euclidean_eqI[where 'a='a] simp: algebra_simps eucl_of_list_inner)

type_synonym 'n rvec = "(real, 'n) vec"

lemma length_mvmult_fa[simp]: "length (mvmult_fa D E xs ys) = D"
  by (auto simp: mvmult_fa_def)

lemma interpret_mvmult_nth:
  assumes "D = CARD('n::enum)"
  assumes "E = CARD('m::enum)"
  assumes "length xs = D * E"
  assumes "length ys = E"
  assumes "n < CARD('n)"
  shows "interpret_floatarith (mvmult_fa D E xs ys ! n) vs =
    ((einterpret xs vs::((real, 'm) vec, 'n) vec) *v einterpret ys vs)  (Basis_list ! n)"
proof -
  have "interpret_floatarith (mvmult_fa D E xs ys ! n) vs = einterpret (mvmult_fa D E xs ys) vs  (Basis_list ! n::'n rvec)"
    using assms
    by (auto simp: eucl_of_list_inner)
  also
  from einterpret_mvmult_fa[OF assms(1,2), of xs ys vs]
  have "einterpret (mvmult_fa D E xs ys) vs = (einterpret xs vs::((real, 'm) vec, 'n) vec) *v einterpret ys vs"
    using assms by simp
  finally show ?thesis by simp
qed


lemmas [simp del] = fold_const_fa.simps

lemma take_eq_map_nth: "n < length xs  take n xs = map ((!) xs) [0..<n]"
  by (induction xs) (auto intro!: nth_equalityI)

lemmas [simp del] = upt_rec_numeral
lemmas map_nth_eq_take = take_eq_map_nth[symmetric]


subsection ‹Definition of Approximating Function using Affine Arithmetic›

lemma interpret_Floatreal: "interpret_floatarith (floatarith.Num f) vs = (real_of_float f)"
  by simp

ML (* Make a congruence rule out of a defining equation for the interpretation

   th is one defining equation of f,
     i.e. th is "f (Cp ?t1 ... ?tn) = P(f ?t1, .., f ?tn)" 
   Cp is a constructor pattern and P is a pattern 

   The result is:
     [|?A1 = f ?t1 ; .. ; ?An= f ?tn |] ==> P (?A1, .., ?An) = f (Cp ?t1 .. ?tn)
       + the a list of names of the A1 .. An, Those are fresh in the ctxt *)

fun mk_congeq ctxt fs th =
  let
    val Const (fN, _) = th |> Thm.prop_of |> HOLogic.dest_Trueprop |> HOLogic.dest_eq
      |> fst |> strip_comb |> fst;
    val ((_, [th']), ctxt') = Variable.import true [th] ctxt;
    val (lhs, rhs) = HOLogic.dest_eq (HOLogic.dest_Trueprop (Thm.prop_of th'));
    fun add_fterms (t as t1 $ t2) =
          if exists (fn f => Term.could_unify (t |> strip_comb |> fst, f)) fs
          then insert (op aconv) t
          else add_fterms t1 #> add_fterms t2
      | add_fterms (t as Abs _) =
          if exists_Const (fn (c, _) => c = fN) t
          then K [t]
          else K []
      | add_fterms _ = I;
    val fterms = add_fterms rhs [];
    val (xs, ctxt'') = Variable.variant_fixes (replicate (length fterms) "x") ctxt';
    val tys = map fastype_of fterms;
    val vs = map Free (xs ~~ tys);
    val env = fterms ~~ vs; (*FIXME*)
    fun replace_fterms (t as t1 $ t2) =
        (case AList.lookup (op aconv) env t of
            SOME v => v
          | NONE => replace_fterms t1 $ replace_fterms t2)
      | replace_fterms t =
        (case AList.lookup (op aconv) env t of
            SOME v => v
          | NONE => t);
    fun mk_def (Abs (x, xT, t), v) =
          HOLogic.mk_Trueprop (HOLogic.all_const xT $ Abs (x, xT, HOLogic.mk_eq (v $ Bound 0, t)))
      | mk_def (t, v) = HOLogic.mk_Trueprop (HOLogic.mk_eq (v, t));
    fun tryext x =
      (x RS @{lemma "(x. f x = g x)  f = g" by blast} handle THM _ => x);
    val cong =
      (Goal.prove ctxt'' [] (map mk_def env)
        (HOLogic.mk_Trueprop (HOLogic.mk_eq (lhs, replace_fterms rhs)))
        (fn {context = goal_ctxt, prems} =>
          Local_Defs.unfold0_tac goal_ctxt (map tryext prems) THEN resolve_tac goal_ctxt [th'] 1))
        RS sym;
    val (cong' :: vars') =
      Variable.export ctxt'' ctxt (cong :: map (Drule.mk_term o Thm.cterm_of ctxt'') vs);
    val vs' = map (fst o fst o Term.dest_Var o Thm.term_of o Drule.dest_term) vars';

  in (vs', cong') end;

fun mk_congs ctxt eqs =
  let
    val fs = fold_rev (fn eq => insert (op =) (eq |> Thm.prop_of |> HOLogic.dest_Trueprop
      |> HOLogic.dest_eq |> fst |> strip_comb
      |> fst)) eqs [];
    val tys = fold_rev (fn f => fold (insert (op =)) (f |> fastype_of |> binder_types |> tl)) fs [];
    val (vs, ctxt') = Variable.variant_fixes (replicate (length tys) "vs") ctxt;
    val subst =
      the o AList.lookup (op =)
        (map2 (fn T => fn v => (T, Thm.cterm_of ctxt' (Free (v, T)))) tys vs);
    fun prep_eq eq =
      let
        val (_, _ :: vs) = eq |> Thm.prop_of |> HOLogic.dest_Trueprop
          |> HOLogic.dest_eq |> fst |> strip_comb;
        val subst = map_filter (fn Var v => SOME (v, subst (#2 v)) | _ => NONE) vs;
      in Thm.instantiate (TVars.empty, Vars.make subst) eq end;
    val (ps, congs) = map_split (mk_congeq ctxt' fs o prep_eq) eqs;
    val bds = AList.make (K ([], [])) tys;
  in (ps ~~ Variable.export ctxt' ctxt congs, bds) end

ML fun interpret_floatariths_congs ctxt =
  mk_congs ctxt @{thms interpret_floatarith.simps interpret_floatariths.simps}
  |> fst
  |> map snd

ML fun preproc_form_conv ctxt =
  Simplifier.rewrite
   (put_simpset HOL_basic_ss ctxt addsimps
     (Named_Theorems.get ctxt @{named_theorems approximation_preproc}))

ML fun reify_floatariths_tac ctxt i =
  CONVERSION (preproc_form_conv ctxt) i
  THEN REPEAT_ALL_NEW (fn i => resolve_tac ctxt (interpret_floatariths_congs ctxt) i) i

method_setup reify_floatariths = Scan.succeed (fn ctxt => SIMPLE_METHOD' (reify_floatariths_tac ctxt)) "reification of floatariths expression"

schematic_goal reify_example:
  "[xs!i * xs!j, xs!i + xs!j powr (sin (xs!0)), xs!k + (2 / 3 * xs!i * xs!j)] = interpret_floatariths ?fas xs"
  by (reify_floatariths)

ML fun interpret_floatariths_step_tac ctxt i = resolve_tac ctxt (interpret_floatariths_congs ctxt) i

method_setup reify_floatariths_step = Scan.succeed (fn ctxt => SIMPLE_METHOD' (interpret_floatariths_step_tac ctxt)) "reification of floatariths expression (step)"

lemma eucl_of_list_interpret_floatariths_cong:
  fixes y::"'a::executable_euclidean_space"
  assumes "b. b  Basis  interpret_floatarith (fa (index Basis_list b)) vs = y  b"
  assumes "length xs = DIM('a)"
  shows "eucl_of_list (interpret_floatariths (map fa [0..<DIM('a)]) vs) = y"
  apply (rule euclidean_eqI)
  apply (subst eucl_of_list_inner)
  by (auto simp: assms)

lemma interpret_floatariths_fold_const_fa[simp]:
  "interpret_floatariths (map fold_const_fa ds) = interpret_floatariths ds"
  by (auto intro!: nth_equalityI)

fun subst_floatarith where
"subst_floatarith s (Add a b)         = Add (subst_floatarith s a) (subst_floatarith s b)" |
"subst_floatarith s (Mult a b)        = Mult (subst_floatarith s a) (subst_floatarith s b)" |
"subst_floatarith s (Minus a)         = Minus (subst_floatarith s a)" |
"subst_floatarith s (Inverse a)       = Inverse (subst_floatarith s a)" |
"subst_floatarith s (Cos a)           = Cos (subst_floatarith s a)" |
"subst_floatarith s (Arctan a)        = Arctan (subst_floatarith s a)" |
"subst_floatarith s (Min a b)         = Min (subst_floatarith s a) (subst_floatarith s b)" |
"subst_floatarith s (Max a b)         = Max (subst_floatarith s a) (subst_floatarith s b)" |
"subst_floatarith s (Abs a)           = Abs (subst_floatarith s a)" |
"subst_floatarith s Pi                = Pi" |
"subst_floatarith s (Sqrt a)          = Sqrt (subst_floatarith s a)" |
"subst_floatarith s (Exp a)           = Exp (subst_floatarith s a)" |
"subst_floatarith s (Powr a b)        = Powr (subst_floatarith s a) (subst_floatarith s b)" |
"subst_floatarith s (Ln a)            = Ln (subst_floatarith s a)" |
"subst_floatarith s (Power a i)       = Power (subst_floatarith s a) i" |
"subst_floatarith s (Floor a)         = Floor (subst_floatarith s a)" |
"subst_floatarith s (Num f)           = Num f" |
"subst_floatarith s (Var n)           = s n"

lemma interpret_floatarith_subst_floatarith:
  assumes "max_Var_floatarith fa  D"
  shows "interpret_floatarith (subst_floatarith s fa) vs =
    interpret_floatarith fa (map (λi. interpret_floatarith (s i) vs) [0..<D])"
  using assms
  by (induction fa) auto

lemma max_Var_floatarith_subst_floatarith_le[THEN order_trans]:
  assumes "length xs  max_Var_floatarith fa"
  shows "max_Var_floatarith (subst_floatarith ((!) xs) fa)  max_Var_floatariths xs"
  using assms
  by (induction fa) (auto intro!: max_Var_floatarith_le_max_Var_floatariths_nth)

lemma max_Var_floatariths_subst_floatarith_le[THEN order_trans]:
  assumes "length xs  max_Var_floatariths fas"
  shows "max_Var_floatariths (map (subst_floatarith ((!) xs)) fas)  max_Var_floatariths xs"
  using assms
  by (induction fas) (auto simp: max_Var_floatarith_subst_floatarith_le)

fun continuous_on_floatarith :: "floatarith  bool" where
  "continuous_on_floatarith (Add a b)         = (continuous_on_floatarith a  continuous_on_floatarith b)" |
"continuous_on_floatarith (Mult a b)        = (continuous_on_floatarith a  continuous_on_floatarith b)" |
"continuous_on_floatarith (Minus a)         = continuous_on_floatarith a" |
"continuous_on_floatarith (Inverse a)       = False" |
"continuous_on_floatarith (Cos a)           = continuous_on_floatarith a" |
"continuous_on_floatarith (Arctan a)        = continuous_on_floatarith a" |
"continuous_on_floatarith (Min a b)         = (continuous_on_floatarith a  continuous_on_floatarith b)" |
"continuous_on_floatarith (Max a b) = (continuous_on_floatarith a  continuous_on_floatarith b)" |
"continuous_on_floatarith (Abs a)           = continuous_on_floatarith a" |
"continuous_on_floatarith Pi                = True" |
"continuous_on_floatarith (Sqrt a)          = False" |
"continuous_on_floatarith (Exp a)           = continuous_on_floatarith a" |
"continuous_on_floatarith (Powr a b)        = False" |
"continuous_on_floatarith (Ln a)            = False" |
"continuous_on_floatarith (Floor a)         = False" |
"continuous_on_floatarith (Power a n)       = (if n = 0 then True else continuous_on_floatarith a)" |
"continuous_on_floatarith (Num f)           = True" |
"continuous_on_floatarith (Var n)           = True"

definition "Maxse xs = fold (λa b. floatarith.Max a b) xs"
definition "norm2e n = Maxse (map (λj. Norm (map (λi. Var (Suc j * n + i)) [0..<n])) [0..<n]) (Num 0)"

definition "Nr l = Num (float_of l)"

lemma interpret_floatarith_Norm:
  "interpret_floatarith (Norm xs) vs = L2_set (λi. interpret_floatarith (xs ! i) vs) {0..<length xs}"
  by (auto simp: Norm_def L2_set_def sum_list_sum_nth power2_eq_square)

lemma interpret_floatarith_Nr[simp]: "interpret_floatarith (Nr U) vs = real_of_float (float_of U)"
  by (auto simp: Nr_def)


fun list_updates where
  "list_updates [] _ xs = xs"
| "list_updates _ [] xs = xs"
| "list_updates (i#is) (u#us) xs = list_updates is us (xs[i:=u])"


lemma list_updates_nth_notmem:
  assumes "length xs = length ys"
  assumes "i  set xs"
  shows "list_updates xs ys vs ! i = vs ! i"
  using assms
  by (induction xs ys arbitrary: i vs rule: list_induct2) auto

lemma list_updates_nth_less:
  assumes "length xs = length ys" "distinct xs"
  assumes "i < length vs"
  shows "list_updates xs ys vs ! i = (if i  set xs then ys ! (index xs i) else vs ! i)"
  using assms
  by (induction xs ys arbitrary: i vs rule: list_induct2) (auto simp: list_updates_nth_notmem)

lemma length_list_updates[simp]: "length (list_updates xs ys vs) = length vs"
  by (induction xs ys vs rule: list_updates.induct) simp_all

lemma list_updates_nth_ge[simp]:
  "x  length vs  length xs = length ys  list_updates xs ys vs ! x = vs ! x"
  apply (induction xs ys vs rule: list_updates.induct)
  apply (auto simp: nth_list_update)
  by (metis list_update_beyond nth_list_update_neq)

lemma
  list_updates_nth:
  assumes [simp]: "length xs = length ys" "distinct xs"
  shows "list_updates xs ys vs ! i = (if i < length vs  i  set xs then ys ! index xs i else vs ! i)" 
  by (auto simp: list_updates_nth_less list_updates_nth_notmem)

lemma list_of_eucl_coord_update:
  assumes [simp]: "length xs = DIM('a::executable_euclidean_space)"
  assumes [simp]: "distinct xs"
  assumes [simp]: "i  Basis"
  assumes [simp]: "n. n  set xs  n < length vs"
  shows "list_updates xs (list_of_eucl (x + (p - x  i) *R i::'a)) vs =
   (list_updates xs (list_of_eucl x) vs)[xs ! index Basis_list i := p]"
  apply (auto intro!: nth_equalityI simp: list_updates_nth nth_list_update)
   apply (simp add: algebra_simps inner_Basis index_nth_id)
  apply (auto simp add: algebra_simps inner_Basis index_nth_id)
  done

definition "eucl_of_env is vs = eucl_of_list (map (nth vs) is)"

lemma list_updates_list_of_eucl_of_env[simp]:
  assumes [simp]: "length xs = DIM('a::executable_euclidean_space)" "distinct xs"
  shows "list_updates xs (list_of_eucl (eucl_of_env xs vs::'a)) vs = vs"
  by (auto intro!: nth_equalityI simp: list_updates_nth nth_list_update eucl_of_env_def)

lemma nth_nth_eucl_of_env_inner:
  "b  Basis  length is = DIM('a)  vs ! (is ! index Basis_list b) = eucl_of_env is vs  b"
  for b::"'a::executable_euclidean_space"
  by (auto simp: eucl_of_env_def eucl_of_list_inner)

lemma list_updates_idem[simp]:
  assumes "(i. i  set X0  i < length vs)"
  shows "(list_updates X0 (map ((!) vs) X0) vs) = vs"
  using assms
  by (induction X0) auto


lemma pairwise_orthogonal_Basis[intro, simp]: "pairwise orthogonal Basis"
  by (auto simp: pairwise_alt orthogonal_def inner_Basis)

primrec freshs_floatarith where
  "freshs_floatarith (Var y) x  (y  set x)"
| "freshs_floatarith (Num a) x  True"
| "freshs_floatarith Pi x  True"
| "freshs_floatarith (Cos a) x  freshs_floatarith a x"
| "freshs_floatarith (Abs a) x  freshs_floatarith a x"
| "freshs_floatarith (Arctan a) x  freshs_floatarith a x"
| "freshs_floatarith (Sqrt a) x  freshs_floatarith a x"
| "freshs_floatarith (Exp a) x  freshs_floatarith a x"
| "freshs_floatarith (Floor a) x  freshs_floatarith a x"
| "freshs_floatarith (Power a n) x  freshs_floatarith a x"
| "freshs_floatarith (Minus a) x  freshs_floatarith a x"
| "freshs_floatarith (Ln a) x  freshs_floatarith a x"
| "freshs_floatarith (Inverse a) x  freshs_floatarith a x"
| "freshs_floatarith (Add a b) x  freshs_floatarith a x  freshs_floatarith b x"
| "freshs_floatarith (Mult a b) x  freshs_floatarith a x  freshs_floatarith b x"
| "freshs_floatarith (floatarith.Max a b) x  freshs_floatarith a x  freshs_floatarith b x"
| "freshs_floatarith (floatarith.Min a b) x  freshs_floatarith a x  freshs_floatarith b x"
| "freshs_floatarith (Powr a b) x  freshs_floatarith a x  freshs_floatarith b x"

lemma freshs_floatarith[simp]:
  assumes "freshs_floatarith fa ds" "length ds = length xs"
  shows "interpret_floatarith fa (list_updates ds xs vs) = interpret_floatarith fa vs"
  using assms
  by (induction fa) (auto simp: list_updates_nth_notmem)

lemma freshs_floatarith_max_Var_floatarithI:
  assumes "x. x  set xs  max_Var_floatarith f  x"
  shows "freshs_floatarith f xs"
  using assms Suc_n_not_le_n
  by (induction f; force)

definition "freshs_floatariths fas xs = (faset fas. freshs_floatarith fa xs)"

lemma freshs_floatariths_max_Var_floatarithsI:
  assumes "x. x  set xs  max_Var_floatariths f  x"
  shows "freshs_floatariths f xs"
  using assms le_trans max_Var_floatarith_le_max_Var_floatariths
  by (force simp: freshs_floatariths_def intro!: freshs_floatarith_max_Var_floatarithI)

lemma freshs_floatariths_freshs_floatarithI:
  assumes "fa. fa  set fas  freshs_floatarith fa xs"
  shows "freshs_floatariths fas xs"
  by (auto simp: freshs_floatariths_def assms)

lemma fresh_floatariths_fresh_floatarithI:
  assumes "freshs_floatariths fas xs"
  assumes "fa  set fas"
  shows "freshs_floatarith fa xs"
  using assms
  by (auto simp: freshs_floatariths_def)

lemma fresh_floatariths_fresh_floatarith[simp]:
  "fresh_floatariths (fas) i  fa  set fas  fresh_floatarith fa i"
  by (induction fas) auto

lemma interpret_floatariths_fresh_cong:
  assumes "i. ¬fresh_floatariths f i  xs ! i = ys ! i"
  shows "interpret_floatariths f ys = interpret_floatariths f xs"
  by (auto intro!: nth_equalityI assms interpret_floatarith_fresh_cong simp: )

fun subterms :: "floatarith  floatarith set" where
"subterms (Add a b) = insert (Add a b) (subterms a  subterms b)" |
"subterms (Mult a b) = insert (Mult a b) (subterms a  subterms b)" |
"subterms (Min a b) = insert (Min a b) (subterms a  subterms b)" |
"subterms (floatarith.Max a b) = insert (floatarith.Max a b) (subterms a  subterms b)" |
"subterms (Powr a b) = insert (Powr a b) (subterms a  subterms b)" |
"subterms (Inverse a) = insert (Inverse a) (subterms a)" |
"subterms (Cos a) = insert (Cos a) (subterms a)" |
"subterms (Arctan a) = insert (Arctan a) (subterms a)" |
"subterms (Abs a) = insert (Abs a) (subterms a)" |
"subterms (Sqrt a) = insert (Sqrt a) (subterms a)" |
"subterms (Exp a) = insert (Exp a) (subterms a)" |
"subterms (Ln a) = insert (Ln a) (subterms a)" |
"subterms (Power a n) = insert (Power a n) (subterms a)" |
"subterms (Floor a) = insert (Floor a) (subterms a)" |
"subterms (Minus a) = insert (Minus a) (subterms a)" |
"subterms Pi = {Pi}" |
"subterms (