(**
  Formalization of the base expression language for the daisy framework
 **)
open preamble miscTheory
open realTheory realLib sptreeTheory
open AbbrevsTheory MachineTypeTheory
open DaisyTactics

val _ = ParseExtras.temp_tight_equality()

val _ = new_theory "Expressions";

val _ = temp_overload_on("abs",``real$abs``);
(**
  EXPRESSIONS WILL use binary operators.
  Define them first
**)
val _ = Datatype `
  binop = Plus | Sub | Mult | Div`;

(**
  Next define an evaluation function for binary operators.
  Errors are added on the expression evaluation level later
*)
val evalBinop_def = Define `
  evalBinop Plus v1 v2 = v1 + v2 /\
  evalBinop Sub  v1 v2 = v1 - v2 /\
  evalBinop Mult v1 v2 = v1 * v2 /\
  evalBinop Div  v1 v2 = v1 / (v2:real)`;

(**
Expressions will use unary operators.
Define them first
**)
val _ = Datatype `
  unop = Neg | Inv`

(**
Define evaluation for unary operators on reals.
Errors are added in the expression evaluation level later
**)
val evalUnop_def = Define `
  evalUnop Neg v = - v /\
  evalUnop Inv v = inv(v:real)`

(**
  Define Expressions parametric over some value type 'v.
  Will ease reasoning about different instantiations later.
**)
val _ = Datatype `
  exp = Var num
      | Const mType 'v
      | Unop unop exp
      | Binop binop exp exp
      | Fma exp exp exp
      | Downcast mType exp`

(**
  Define evaluation of FMA (fused multiply-add) on reals
  Errors are added on the expression evaluation level later
**)
val evalFma_def = Define `
  evalFma v1 v2 v3 = evalBinop Plus v1 (evalBinop Mult v2 v3)`

val toREval_def = Define `
  toREval e :real exp =
    case e of
      | (Var n) => Var n
      | (Const m n) => Const M0 n
      | (Unop u e1) => Unop u (toREval e1)
      | (Binop bop e1 e2) => Binop bop (toREval e1) (toREval e2)
      | (Fma e1 e2 e3) => Fma (toREval e1) (toREval e2) (toREval e3)
      | (Downcast m e1) => (toREval e1)`;

(**
  Define a perturbation function to ease writing of basic definitions
**)
val perturb_def = Define `
  perturb (r:real) (e:real) = r * (1 + e)`

(**
Define expression evaluation relation parametric by an "error" epsilon.
The result value expresses float computations according to the IEEE standard,
using a perturbation of the real valued computation by (1 + delta), where
|delta| <= machine epsilon.
**)
val (eval_exp_rules, eval_exp_ind, eval_exp_cases) = Hol_reln `
  (!E defVars m x v.
     defVars x = SOME m /\
     E x = SOME v ==>
     eval_exp E defVars (Var x) v m) /\
  (!E defVars m n delta.
      abs delta <= (mTypeToQ m) ==>
     eval_exp E defVars (Const m n) (perturb n delta) m) /\
  (!E defVars m f1 v1.
     eval_exp E defVars f1 v1 m ==>
     eval_exp E defVars (Unop Neg f1) (evalUnop Neg v1) m) /\
  (!E defVars m f1 v1 delta.
     abs delta <= (mTypeToQ m) /\
     (v1 <> 0) /\
     eval_exp E defVars f1 v1 m ==>
     eval_exp E defVars (Unop Inv f1) (perturb (evalUnop Inv v1) delta) m) /\
  (!E defVars m m1 f1 v1 delta.
     isMorePrecise m1 m /\
     abs delta <= (mTypeToQ m) /\
     eval_exp E defVars f1 v1 m1 ==>
     eval_exp E defVars (Downcast m f1) (perturb v1 delta) m) /\
  (!E defVars m1 m2 b f1 f2 v1 v2 delta.
     abs delta <= (mTypeToQ (join m1 m2)) /\
     eval_exp E defVars f1 v1 m1 /\
     eval_exp E defVars f2 v2 m2 /\
     ((b = Div) ==> (v2 <> 0)) ==>
     eval_exp E defVars (Binop b f1 f2) (perturb (evalBinop b v1 v2) delta) (join m1 m2)) /\
  (!E defVars m1 m2 m3 f1 f2 f3 v1 v2 v3 delta.
     abs delta <= (mTypeToQ (join3 m1 m2 m3)) /\
     eval_exp E defVars f1 v1 m1 /\
     eval_exp E defVars f2 v2 m2 /\
     eval_exp E defVars f3 v3 m3 ==>
     eval_exp E defVars (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) delta) (join3 m1 m2 m3))`;

val eval_exp_cases_old = save_thm ("eval_exp_cases_old", eval_exp_cases);

(**
  Generate a better case lemma
**)
val eval_exp_cases =
  map (GEN_ALL o SIMP_CONV (srw_ss()) [Once eval_exp_cases])
    [``eval_exp E defVars (Var v) res m``,
     ``eval_exp E defVars (Const m1 n) res m2``,
     ``eval_exp E defVars (Unop u e) res m``,
     ``eval_exp E defVars (Binop n e1 e2) res m``,
     ``eval_exp E defVars (Fma e1 e2 e3) res m``,
     ``eval_exp E defVars (Downcast m1 e1) res m2``]
  |> LIST_CONJ |> curry save_thm "eval_exp_cases";

val [Var_load, Const_dist, Unop_neg, Unop_inv, Downcast_dist, Binop_dist, Fma_dist] = CONJ_LIST 7 eval_exp_rules;
save_thm ("Var_load", Var_load);
save_thm ("Const_dist", Const_dist);
save_thm ("Unop_neg", Unop_neg);
save_thm ("Unop_inv", Unop_inv);
save_thm ("Binop_dist", Binop_dist);
save_thm ("Fma_dist", Fma_dist);
save_thm ("Downcast_dist", Downcast_dist);

(**
  Show some simpler (more general) rule lemmata
**)

val Const_dist' = store_thm (
  "Const_dist'",
  ``!m n delta v m' E Gamma.
      (abs delta) <= (mTypeToQ m) /\
      v = perturb n delta /\
      m' = m ==>
      eval_exp E Gamma (Const m n) v m'``,
  fs [Const_dist]);

val Unop_neg' = store_thm (
  "Unop_neg'",
  ``!m f1 v1 v m' E Gamma.
      eval_exp E Gamma f1 v1 m /\
      v = evalUnop Neg v1 /\
      m' = m ==>
      eval_exp E Gamma (Unop Neg f1) v m'``,
  fs [Unop_neg]);

val Unop_inv' = store_thm (
  "Unop_inv'",
  ``!m f1 v1 delta v m' E Gamma.
      (abs delta) <= (mTypeToQ m) /\
      eval_exp E Gamma f1 v1 m /\
      (v1 <> 0) /\
      v = perturb (evalUnop Inv v1) delta /\
      m' = m ==>
      eval_exp E Gamma (Unop Inv f1) v m'``,
  fs [Unop_inv]);

val Downcast_dist' = store_thm ("Downcast_dist'",
  ``!m m1 f1 v1 delta v m' E Gamma.
      isMorePrecise m1 m /\
      (abs delta) <= (mTypeToQ m) /\
      eval_exp E Gamma f1 v1 m1 /\
      v = perturb v1 delta /\
      m' = m ==>
      eval_exp E Gamma (Downcast m f1) v m'``,
  rpt strip_tac
  \\ rw []
  \\ match_mp_tac Downcast_dist
  \\ qexists_tac `m1` \\ fs[]);

val Binop_dist' = store_thm ("Binop_dist'",
  ``!m1 m2 op f1 f2 v1 v2 delta v m' E Gamma.
      (abs delta) <= (mTypeToQ m') /\
      eval_exp E Gamma f1 v1 m1 /\
      eval_exp E Gamma f2 v2 m2 /\
      ((op = Div) ==> (v2 <> 0)) /\
      v = perturb (evalBinop op v1 v2) delta /\
      m' = join m1 m2 ==>
      eval_exp E Gamma (Binop op f1 f2) v m'``,
  fs [Binop_dist]);

val Fma_dist' = store_thm ("Fma_dist'",
  ``!m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m' E Gamma.
      (abs delta) <= (mTypeToQ m') /\
      eval_exp E Gamma f1 v1 m1 /\
      eval_exp E Gamma f2 v2 m2 /\
      eval_exp E Gamma f3 v3 m3 /\
      v = perturb (evalFma v1 v2 v3) delta /\
      m' = join3 m1 m2 m3 ==>
      eval_exp E Gamma (Fma f1 f2 f3) v m'``,
  fs [Fma_dist]);

(**
  Define the set of "used" variables of an expression to be the set of variables
  occuring in it
**)
val usedVars_def = Define `
  usedVars (e: 'a exp) :num_set =
    case e of
      | Var x => insert x () (LN)
      | Unop u e1 => usedVars e1
      | Binop b e1 e2 => union (usedVars e1) (usedVars e2)
      | Fma e1 e2 e3 => union (usedVars e1) (union (usedVars e2) (usedVars e3))
      | Downcast m e1 => usedVars e1
      | _ => LN`;

(**
  If |delta| <= 0 then perturb v delta is exactly v.
**)
val delta_0_deterministic = store_thm("delta_0_deterministic",
  ``!(v:real) (delta:real). abs delta <= 0 ==> perturb v delta = v``,
  fs [perturb_def,ABS_BOUNDS,REAL_LE_ANTISYM]);

val delta_M0_deterministic = store_thm("delta_M0_deterministic",
  ``!(v:real) (delta:real). abs delta <= mTypeToQ M0 ==> perturb v delta = v``,
  fs [mTypeToQ_def,perturb_def,ABS_BOUNDS,REAL_LE_ANTISYM]);

val toRMap_def = Define `
  toRMap (d:num -> mType option) (n:num) : mType option =
    case d n of
      | SOME m => SOME M0
      | NONE => NONE`;

val toRMap_eval_M0 = store_thm (
  "toRMap_eval_M0",
  ``!f v E Gamma m.
      eval_exp E (toRMap Gamma) (toREval f) v m ==> m = M0``,
  Induct \\ simp[Once toREval_def] \\ fs[eval_exp_cases, toRMap_def]
  \\ rpt strip_tac \\ fs[]
  >- (every_case_tac \\ fs[])
  >- (rveq \\ first_x_assum drule \\ strip_tac \\ fs[])
  >- (rveq \\ first_x_assum drule \\ strip_tac \\ fs[])
  >- (`m1 = M0` by (rpt (first_x_assum drule \\ strip_tac) \\ fs[])
      \\ `m2 = M0` by (rpt (first_x_assum drule \\ strip_tac) \\ fs[])
      \\ rveq \\ fs[join_def])
  >- (`m1 = M0` by (rpt (first_x_assum drule \\ strip_tac) \\ fs[])
      \\ `m2 = M0` by (rpt (first_x_assum drule \\ strip_tac) \\ fs[])
      \\ `m3 = M0` by (rpt (first_x_assum drule \\ strip_tac) \\ fs[])
      \\ rveq \\ fs[join3_def] \\ fs[join_def]));

(**
Evaluation with 0 as machine epsilon is deterministic
**)
val meps_0_deterministic = store_thm("meps_0_deterministic",
  ``!(f: real exp) v1:real v2:real E defVars.
      eval_exp E (toRMap defVars) (toREval f) v1 M0 /\
      eval_exp E (toRMap defVars) (toREval f) v2 M0 ==>
      v1 = v2``,
  Induct_on `f`
  >- (rw [toREval_def] \\ fs [eval_exp_cases])
  >- (rw [toREval_def]
      \\ fs [eval_exp_cases, mTypeToQ_def, delta_0_deterministic])
  >- (rw []
      \\ rpt (
           qpat_x_assum `eval_exp _ _ (toREval _) _ _`
             (fn thm => assume_tac (ONCE_REWRITE_RULE [toREval_def] thm)))
      \\ Cases_on `u` \\ fs [eval_exp_cases] \\ rw []
      \\ fs [eval_exp_cases]
      >- (res_tac \\ fs [REAL_NEG_EQ])
      >- (res_tac \\ fs [mTypeToQ_def, delta_0_deterministic]))
  >- (rw[]
      \\ rename1 `Binop b f1 f2`
      \\ rpt (
           qpat_x_assum `eval_exp _ _ (toREval _) _ _`
            (fn thm => assume_tac (ONCE_REWRITE_RULE [toREval_def] thm)))
      \\ Cases_on `b` \\ fs [eval_exp_cases]
      \\ `m1 = M0 /\ m2 = M0` by (conj_tac \\ irule  toRMap_eval_M0 \\ asm_exists_tac \\ fs [])
      \\ rw[]
      \\ rename1 `eval_exp E _ (toREval f1) vf11 M0`
      \\ rename1 `eval_exp E _ (toREval f1) vf12 m1`
      \\ rename1 `eval_exp E _ (toREval f2) vf21 M0`
      \\ rename1 `eval_exp _ _ (toREval f2) vf22 m2`
      \\ `m1 = M0 /\ m2 = M0` by (conj_tac \\ irule toRMap_eval_M0 \\ asm_exists_tac \\  fs [])
      \\ rw []
      \\ fs [join_def, mTypeToQ_def, delta_0_deterministic]
      \\ qpat_x_assum `!v1 v2 E defVars. _ /\ _ ==> v1 = v2` (fn thm =>qspecl_then [`vf21`,`vf22`] ASSUME_TAC thm)
      \\ qpat_x_assum `!v1 v2 E defVars. _ /\ _ ==> v1 = v2` (fn thm =>qspecl_then [`vf11`,`vf12`] ASSUME_TAC thm)
      \\ res_tac
      \\ rw[])
  >- (rw[]
      \\ rpt (
           qpat_x_assum `eval_exp _ _ (toREval _) _ _`
            (fn thm => assume_tac (ONCE_REWRITE_RULE [toREval_def] thm)))
      \\ fs [eval_exp_cases]
      \\ `m1 = M0 /\ m2 = M0` by (conj_tac \\ irule  toRMap_eval_M0 \\ asm_exists_tac \\ fs [])
      \\ `m3 = M0` by (irule  toRMap_eval_M0 \\ asm_exists_tac \\ fs [])
      \\ `m1' = M0 /\ m2' = M0` by (conj_tac \\ irule  toRMap_eval_M0 \\ asm_exists_tac \\ fs [])
      \\ `m3' = M0` by (irule  toRMap_eval_M0 \\ asm_exists_tac \\ fs [])
      \\ rw[]
      \\ qpat_x_assum `!v1 v2 E defVars. _ /\ _ ==> v1 = v2` (fn thm =>qspecl_then [`v3`,`v3'`, `E`, `defVars`] ASSUME_TAC thm)
      \\ qpat_x_assum `!v1 v2 E defVars. _ /\ _ ==> v1 = v2` (fn thm =>qspecl_then [`v2'`,`v2''`, `E`, `defVars`] ASSUME_TAC thm)
      \\ qpat_x_assum `!v1 v2 E defVars. _ /\ _ ==> v1 = v2` (fn thm =>qspecl_then [`v1'`,`v1''`, `E`, `defVars`] ASSUME_TAC thm)
      \\ fs [join3_def, join_def, mTypeToQ_def, delta_0_deterministic])
  >- (rw[]
      \\ rpt (
           qpat_x_assum `eval_exp _ _ (toREval (Downcast _ _)) _ _`
            (fn thm => assume_tac (ONCE_REWRITE_RULE [toREval_def] thm)))
      \\ fs [eval_exp_cases]
      \\ res_tac));

(**
Helping lemmas. Needed in soundness proof.
For each evaluation of using an arbitrary epsilon, we can replace it by
evaluating the subExpressions and then binding the result values to different
variables in the Environment.
**)
val binary_unfolding = store_thm("binary_unfolding",
  ``!(b:binop) (f1:(real)exp) (f2:(real)exp) E Gamma (v:real) v1 v2 m1 m2 delta.
      (b = Div ==> (v2 <> 0)) /\
      (abs delta) <= (mTypeToQ (join m1 m2)) /\
      eval_exp E Gamma f1 v1 m1 /\
      eval_exp E Gamma f2 v2 m2 /\
      eval_exp E Gamma (Binop b f1 f2) (perturb (evalBinop b v1 v2) delta) (join m1 m2) ==>
      eval_exp (updEnv 2 v2 (updEnv 1 v1 emptyEnv))
        (updDefVars 2 m2 (updDefVars 1 m1 Gamma))
        (Binop b (Var 1) (Var 2)) (perturb (evalBinop b v1 v2) delta) (join m1 m2)``,
  fs [updEnv_def,updDefVars_def,join_def,eval_exp_cases,APPLY_UPDATE_THM,PULL_EXISTS]
  \\ metis_tac []);

val fma_unfolding = store_thm("fma_unfolding",
  ``!(f1:(real)exp) (f2:(real)exp) (f3:(real)exp) E Gamma (v:real) v1 v2 v3 m1 m2 m3 delta.
      (abs delta) <= (mTypeToQ (join3 m1 m2 m3)) /\
      eval_exp E Gamma f1 v1 m1 /\
      eval_exp E Gamma f2 v2 m2 /\
      eval_exp E Gamma f3 v3 m3 /\
      eval_exp E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) delta) (join3 m1 m2 m3) ==>
      eval_exp (updEnv 3 v3 (updEnv 2 v2 (updEnv 1 v1 emptyEnv)))
        (updDefVars 3 m3 (updDefVars 2 m2 (updDefVars 1 m1 Gamma)))
        (Fma (Var 1) (Var 2) (Var 3)) (perturb (evalFma v1 v2 v3) delta) (join3 m1 m2 m3)``,
  fs [updEnv_def,updDefVars_def,join3_def,join_def,eval_exp_cases,APPLY_UPDATE_THM,PULL_EXISTS]
  \\ rpt strip_tac
  \\ qexists_tac `delta'`
  \\ conj_tac \\ fs[]);

val eval_eq_env = store_thm (
  "eval_eq_env",
  ``!e E1 E2 Gamma v m.
      (!x. E1 x = E2 x) /\
      eval_exp E1 Gamma e v m ==>
      eval_exp E2 Gamma e v m``,
  Induct \\ rpt strip_tac \\ fs[eval_exp_cases]
  >- (`E1 n = E2 n` by (first_x_assum irule)
      \\ fs[])
  >- (qexists_tac `delta'` \\ fs[])
  >- (rveq \\ qexists_tac `v1` \\ fs[]
      \\ first_x_assum drule \\ disch_then irule \\ fs[])
  >- (rveq \\ qexists_tac `v1` \\ fs[]
      \\ qexists_tac `delta'` \\ fs[]
      \\ first_x_assum drule \\ disch_then irule \\ fs[])
  >- (rveq \\ qexistsl_tac [`m1`, `m2`, `v1`, `v2`, `delta'`]
      \\ fs[] \\ conj_tac \\ first_x_assum irule \\ asm_exists_tac \\ fs[])
  >- (rveq \\ qexistsl_tac [`m1`, `m2`, `m3`, `v1`, `v2`, `v3`, `delta'`]
      \\ fs[] \\ prove_tac [])
  >- (rveq \\ qexistsl_tac [`m1'`, `v1`, `delta'`] \\ fs[]
      \\ first_x_assum drule \\ disch_then irule \\ fs[]));

(* (** *)
(*   Analogous lemma for unary expressions *)
(* **) *)
(* val unary_unfolding = store_thm("unary_unfolding", *)
(* ``!(u:unop) (e1:(real)exp) (m:mType) E defVars (v:real). *)
(*        (eval_exp E defVars (Unop Inv e1) v m <=> *)
(*        (?(v1:real). *)
(*           eval_exp E defVars e1 v1 m /\ *)
(*           eval_exp (updEnv 1 v1 emptyEnv) (updDefVars 1 m defVars) (Unop Inv (Var 1)) v m))``, *)
(*   fs [updEnv_def, updDefVars_def, eval_exp_cases,APPLY_UPDATE_THM,PULL_EXISTS] *)
(*   \\ metis_tac []); *)

(*
  Using the parametric Expressions, define boolean Expressions for conditionals
*)
(* val _ = Datatype ` *)
(*   bexp = Leq ('v exp) ('v exp) *)
(*        | Less ('v exp) ('v exp)` *)

(*
  Define evaluation of boolean expressions
*)
(* val (bval_exp_rules, bval_exp_ind, bval_exp_cases) = Hol_reln ` *)
(*   (!E defVars e1 e2 v1 v2 m. *)
(*       eval_exp E defVars e1 v1 m /\ *)
(*       eval_exp E defVars e2 v2 m ==> *)
(*       bval E defVars m (Leq e1 e2) (v1 <= v2))/\ *)
(*   (!E defVars e1 e2 v1 v2 m. *)
(*       eval_exp E defVars e1 v1 m /\ *)
(*       eval_exp E defVars e2 v2 m ==> *)
(*       bval E defVars m (Less e1 e2) (v1 < v2))`; *)

(* val bval_exp_cases = *)
(*   map (GEN_ALL o SIMP_CONV (srw_ss()) [Once bval_exp_cases]) *)
(*     [``bval E defVars m (Leq e1 e2) res``, *)
(*      ``bval E defVars m (Less e1 e2) res``] *)
(*   |> LIST_CONJ |> curry save_thm "bval_exp_cases"; *)

(* (** *)
(*   Simplify arithmetic later by making > >= only abbreviations *)
(* **) *)
(* val _ = overload_on("Gr",``\(e1:('v)exp). \(e2:('v)exp). Less e2 e1``); *)
(* val _ = overload_on("Greq",``\(e1:('v)exp). \(e2:('v)exp). Leq e2 e1``); *)

val _ = export_theory();