(**
Formalization of the base expression language for the daisy framework
 **)
Require Import Coq.Reals.Reals Interval.Interval_tactic.
Require Import Daisy.Infra.RealConstruction.
Set Implicit Arguments.
(**
  Expressions will use binary operators.
  Define them first
**)
Inductive binop : Type := Plus | Sub | Mult | Div.
(**
  Next define an evaluation function for binary operators on reals.
  Errors are added on the expression evaluation level later.
 **)
Fixpoint eval_binop (o:binop) (v1:R) (v2:R) :=
  match o with
  | Plus => Rplus v1 v2
  | Sub => Rminus v1 v2
  | Mult => Rmult v1 v2
  | Div => Rdiv v1 v2
  end.
(*
  Define expressions parametric over some value type V.
  Will ease reasoning about different instantiations later.
*)
Inductive exp (V:Type): Type :=
  Var: nat -> exp V
| Const: V -> exp V
| Binop: binop -> exp V -> exp V -> exp V.
(**
  Define the machine epsilon for floating point operations.
  FIXME: Currently set to 1.0 instead of the concrete value!
 **)
Definition machineEpsilon:R := realFromNum 1 1 53.
(* Definition machineEpsilon:R := realFromNum 11102230246251565 16 16. *)
(**
  Define a perturbation function to ease writing of basic definitions
**)
Definition perturb (r:R) (e:R) :=
  Rmult r (Rplus 1 e).
(**
  Abbreviation for the type of an environment
 **)
Definition env_ty := nat -> R.
(**
  Define expression evaluation relation parametric by an "error" delta.
  This value will be used later to express float computations using a perturbation
  of the real valued computation by (1 + d)
**)
Inductive eval_exp (eps:R) (env:env_ty) : (exp R) -> R -> Prop :=
  Var_load x: eval_exp eps env (Var R x) (env x)
| Const_dist n delta:
    Rle (Rabs delta) eps ->
    eval_exp eps env (Const n) (perturb n delta)
|Binop_dist op e1 e2 v1 v2 delta:
   Rle (Rabs delta) eps ->
                eval_exp eps env e1 v1 ->
                eval_exp eps env e2 v2 ->
                eval_exp eps env (Binop op e1 e2) (perturb (eval_binop op v1 v2) delta).
(**
  Define real evaluation as stated above:
 **)
Definition is_real_value (e:exp R) (env:env_ty) (v:R) := eval_exp R0 env e v.
(**
  Prove that using eps = 0 makes the evaluation deterministic
 **)
Lemma Rabs_0_impl_eq (d:R):
  Rle (Rabs d) R0 -> d = R0.
Proof.
  intros abs_leq_0.
  pose proof (Rabs_pos d) as abs_geq_0.
  pose proof (Rle_antisym (Rabs d) R0 abs_leq_0 abs_geq_0) as Rabs_eq.
  rewrite <- Rabs_R0 in Rabs_eq.
  apply Rsqr_eq_asb_1 in Rabs_eq.
  rewrite Rsqr_0 in Rabs_eq.
  apply Rsqr_0_uniq in Rabs_eq; assumption.
Qed.

Lemma eval_det (e:exp R) (env:env_ty):
  forall v1 v2,
  eval_exp R0 env e v1 ->
  eval_exp R0 env e v2 ->
  v1 = v2.
Proof.
  induction e; intros v1 v2 eval_v1 eval_v2; inversion eval_v1; inversion eval_v2; try auto.
  - apply Rabs_0_impl_eq in H0;
      apply Rabs_0_impl_eq in H3.
    rewrite H0, H3; reflexivity.
  - apply Rabs_0_impl_eq in H2;
      apply Rabs_0_impl_eq in H9.
    rewrite H2, H9.
    subst.
    rewrite (IHe1 v0 v4); auto.
    rewrite (IHe2 v3 v5); auto.
Qed.

Lemma perturb_0_val (v:R) (delta:R):
  (Rabs delta <= 0)%R ->
  perturb v delta = v.
Proof.
  intros abs_0; apply Rabs_0_impl_eq in abs_0; subst.
  unfold perturb.
  rewrite Rmult_plus_distr_l.
  rewrite Rmult_0_r.
  rewrite Rmult_1_r.
  rewrite Rplus_0_r; auto.
Qed.

Lemma Rsub_eq_Ropp_Rplus (a:R) (b:R) :
  (a - b = a + (- b))%R.
Proof.
  field_simplify; reflexivity.
Qed.

Lemma const_abs_err_bounded (n:R) (nR:R) (nF:R):
  forall cenv:nat -> R,
  eval_exp 0%R cenv (Const n) nR ->
  eval_exp machineEpsilon cenv (Const n) nF ->
  (Rabs (nR - nF) <= Rabs (n * machineEpsilon))%R.
Proof.
  intros cenv eval_real eval_float.
  inversion eval_real; subst.
  rewrite perturb_0_val; auto.
  inversion eval_float; subst.
  unfold perturb; simpl.
  rewrite Rmult_plus_distr_l.
  rewrite Rmult_1_r.
  rewrite Rsub_eq_Ropp_Rplus.
  rewrite Ropp_plus_distr.
  assert (-n + (- (n * delta0)) = ( -(n * delta0) + - n))%R by (rewrite Rplus_comm; auto).
  rewrite H.
  rewrite <- Rsub_eq_Ropp_Rplus.
  rewrite Rplus_minus.
  rewrite Rabs_Ropp.
  repeat rewrite Rabs_mult.
  apply Rmult_le_compat_l.
  - apply Rabs_pos.
  - assert (Rabs machineEpsilon = machineEpsilon).
    + unfold machineEpsilon.
      unfold realFromNum, negativePower.
      unfold Rabs.
      destruct Rcase_abs; auto.
      exfalso.
      apply (Rlt_not_ge _ _ r).
      interval.
    + rewrite H2; auto.
Qed.

Lemma var_abs_err_bounded (n:nat) (nR:R) (nF:R) (cenv:nat->R) (nlo:R) (nhi:R):
  (nlo <= cenv n <= nhi)%R ->
  eval_exp 0%R cenv (Var R n) nR ->
  eval_exp machineEpsilon cenv (Var R n) nF ->
  (Rabs (nR - nF) <= Rabs ((Rmax (Rabs nlo) (Rabs nhi)) * machineEpsilon))%R.
Proof.
  intros [lo_le_env env_le_hi] eval_real eval_float.
  inversion eval_real; subst.
  (* rewrite perturb_0_val. *)
Admitted.

(**
  Using the parametric expressions, define boolean expressions for conditionals
**)
Inductive bexp (V:Type) : Type :=
  leq: exp V -> exp V -> bexp V
| less: exp V -> exp V -> bexp V.
(**
  Define evaluation of booleans for reals
 **)
Inductive bval (eps:R) (env:env_ty) : (bexp R) -> Prop -> Prop :=
  leq_eval (e1:exp R) (e2:exp R) (v1:R) (v2:R):
    eval_exp eps env e1 v1 ->
    eval_exp eps env e2 v2 ->
    bval eps env (leq e1 e2) (Rle v1 v2)
|less_eval (e1:exp R) (e2:exp R) (v1:R) (v2:R):
    eval_exp eps env e1 v1 ->
    eval_exp eps env e2 v2 ->
    bval eps env (less e1 e2) (Rlt v1 v2).
(**
 Simplify arithmetic later by making > >= only abbreviations
**)
Definition gr := fun (V:Type) (e1: exp V) (e2: exp V) => less e2 e1.
Definition greq := fun (V:Type) (e1:exp V) (e2: exp V) => leq e2 e1.