Refactoring and prove error bounds for addition and subtraction

coq/ErrorBounds.v

+Require Import Coq.Reals.Reals Coq.micromega.Psatz Interval.Interval_tactic.
+Require Import Daisy.Infra.Abbrevs Daisy.Infra.RealSimps Daisy.Expressions.
+(**
+  TODO: Check wether we need Rabs (n * machineEpsilon) instead
+**)
+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 Rabs_err_simpl.
+  rewrite Rabs_mult.
+  apply Rmult_le_compat_l; [apply Rabs_pos | auto].
+Qed.
+
+(**
+  TODO: Maybe improve bound by adding interval constraint and proving that it is leq than maxAbs of bounds
+  (nlo <= cenv n <= nhi)%R ->   (Rabs (nR - nF) <= Rabs ((Rmax (Rabs nlo) (Rabs nhi)) * machineEpsilon))%R.
+**)
+Lemma param_abs_err_bounded (n:nat) (nR:R) (nF:R) (cenv:nat->R):
+  eval_exp 0%R cenv (Param R n) nR ->
+  eval_exp machineEpsilon cenv (Param R n) nF ->
+  (Rabs (nR - nF) <= Rabs (cenv n) * machineEpsilon)%R.
+Proof.
+  intros eval_real eval_float.
+  inversion eval_real; subst.
+  rewrite perturb_0_val; auto.
+  inversion eval_float; subst.
+  unfold perturb; simpl.
+  rewrite Rabs_err_simpl.
+  repeat rewrite Rabs_mult.
+  apply Rmult_le_compat_l; [ apply Rabs_pos | auto].
+Qed.
+
+Lemma add_abs_err_bounded (e1:exp R) (e1R:R) (e1F:R) (e2:exp R) (e2R:R) (e2F:R) (vR:R) (vF:R) (cenv:nat->R) (err1:R) (err2:R):
+  eval_exp 0%R cenv e1 e1R ->
+  eval_exp machineEpsilon cenv e1 e1F ->
+  eval_exp 0%R cenv e2 e2R ->
+  eval_exp machineEpsilon cenv e2 e2F ->
+  eval_exp 0%R cenv (Binop Plus e1 e2) vR ->
+  eval_exp machineEpsilon (updEnv 2 e2F (updEnv 1 e1F cenv))  (Binop Plus (Var R 1) (Var R 2)) vF ->
+  (Rabs (e1R - e1F) <= err1)%R ->
+  (Rabs (e2R - e2F) <= err2)%R ->
+  (Rabs (vR - vF) <= err1 + err2 + (Rabs (e1F + e2F) * machineEpsilon))%R.
+Proof.
+  intros e1_real e1_float e2_real e2_float plus_real plus_float bound_e1 bound_e2.
+  (* Prove that e1R and e2R are the correct values and that vR is e1R + e2R *)
+  inversion plus_real; subst.
+  rewrite perturb_0_val in plus_real; auto.
+  rewrite (perturb_0_val (eval_binop Plus v1 v2) delta); auto.
+  unfold eval_binop in *; simpl in *.
+  clear delta H2.
+  rewrite (eval_det H4 e1_real);
+    rewrite (eval_det H5 e2_real).
+  rewrite (eval_det H4 e1_real) in plus_real.
+  rewrite (eval_det H5 e2_real) in plus_real.
+  clear H4 H5 v1 v2.
+  (* Now unfold the float valued evaluation to get the deltas we need for the inequality *)
+  inversion plus_float; subst.
+  unfold perturb; simpl.
+  inversion H4; subst; inversion H5; subst.
+  unfold updEnv; simpl.
+  (* We have now obtained all necessary values from the evaluations --> remove them for readability *)
+  clear plus_float H4 H5 plus_real e1_real e1_float e2_real e2_float.
+  repeat rewrite Rmult_plus_distr_l.
+  rewrite Rmult_1_r.
+  rewrite Rsub_eq_Ropp_Rplus.
+  repeat rewrite Ropp_plus_distr.
+  rewrite plus_bounds_simplify.
+  pose proof (Rabs_triang (e1R + - e1F) ((e2R + - e2F) + - ((e1F + e2F) * delta))).
+  rewrite Rplus_assoc.
+  eapply Rle_trans.
+  apply H.
+  pose proof (Rabs_triang (e2R + - e2F) (- ((e1F + e2F) * delta))).
+  pose proof (Rplus_le_compat_l (Rabs (e1R + - e1F)) _ _ H0).
+  eapply Rle_trans.
+  apply H1.
+  rewrite <- Rplus_assoc.
+  repeat rewrite <- Rsub_eq_Ropp_Rplus.
+  rewrite Rabs_Ropp.
+  eapply Rplus_le_compat.
+  - eapply Rplus_le_compat; auto.
+  - rewrite Rabs_mult.
+    eapply Rle_trans.
+    eapply Rmult_le_compat_l.
+    apply Rabs_pos.
+    apply H2.
+    apply Req_le; auto.
+Qed.
+
+(**
+  Copy-Paste proof with minor differences, was easier then manipulating the evaluations and then applying the lemma
+**)
+Lemma substract_abs_err_bounded (e1:exp R) (e1R:R) (e1F:R) (e2:exp R) (e2R:R) (e2F:R) (vR:R) (vF:R) (cenv:nat->R) (err1:R) (err2:R):
+  eval_exp 0%R cenv e1 e1R ->
+  eval_exp machineEpsilon cenv e1 e1F ->
+  eval_exp 0%R cenv e2 e2R ->
+  eval_exp machineEpsilon cenv e2 e2F ->
+  eval_exp 0%R cenv (Binop Sub e1 e2) vR ->
+  eval_exp machineEpsilon (updEnv 2 e2F (updEnv 1 e1F cenv))  (Binop Sub (Var R 1) (Var R 2)) vF ->
+  (Rabs (e1R - e1F) <= err1)%R ->
+  (Rabs (e2R - e2F) <= err2)%R ->
+  (Rabs (vR - vF) <= err1 + err2 + (Rabs (e1F - e2F) * machineEpsilon))%R.
+Proof.
+  intros e1_real e1_float e2_real e2_float sub_real sub_float bound_e1 bound_e2.
+  (* Prove that e1R and e2R are the correct values and that vR is e1R + e2R *)
+  inversion sub_real; subst.
+  rewrite perturb_0_val in sub_real; auto.
+  rewrite perturb_0_val; auto.
+  unfold eval_binop in *; simpl in *.
+  clear delta H2.
+  rewrite (eval_det H4 e1_real);
+    rewrite (eval_det H5 e2_real).
+  rewrite (eval_det H4 e1_real) in sub_real.
+  rewrite (eval_det H5 e2_real) in sub_real.
+  clear H4 H5 v1 v2.
+  (* Now unfold the float valued evaluation to get the deltas we need for the inequality *)
+  inversion sub_float; subst.
+  unfold perturb; simpl.
+  inversion H4; subst; inversion H5; subst.
+  unfold updEnv; simpl.
+  (* We have now obtained all necessary values from the evaluations --> remove them for readability *)
+  clear sub_float H4 H5 sub_real e1_real e1_float e2_real e2_float.
+  repeat rewrite Rmult_plus_distr_l.
+  rewrite Rmult_1_r.
+  repeat rewrite Rsub_eq_Ropp_Rplus.
+  repeat rewrite Ropp_plus_distr.
+  rewrite plus_bounds_simplify.
+  pose proof (Rabs_triang (e1R + - e1F) ((- e2R + - - e2F) + - ((e1F + - e2F) * delta))).
+  rewrite Rplus_assoc.
+  eapply Rle_trans.
+  apply H.
+  pose proof (Rabs_triang (- e2R + - - e2F) (- ((e1F + - e2F) * delta))).
+  pose proof (Rplus_le_compat_l (Rabs (e1R + - e1F)) _ _ H0).
+  eapply Rle_trans.
+  apply H1.
+  rewrite <- Rplus_assoc.
+  repeat rewrite <- Rsub_eq_Ropp_Rplus.
+  rewrite Rabs_Ropp.
+  assert (Rabs (- e2R - - e2F)%R = Rabs (e2R - e2F)).
+  - rewrite Rsub_eq_Ropp_Rplus.
+    rewrite <- Ropp_plus_distr.
+    rewrite Rabs_Ropp.
+    rewrite <- Rsub_eq_Ropp_Rplus; auto.
+  - rewrite H3.
+    eapply Rplus_le_compat.
+    + eapply Rplus_le_compat; auto.
+    +  rewrite Rabs_mult.
+       eapply Rle_trans.
+       eapply Rmult_le_compat_l.
+       apply Rabs_pos.
+       apply H2.
+       apply Req_le; auto.
+Qed.
\ No newline at end of file

coq/Infra/Abbrevs.v

@@ -23,4 +23,14 @@ Definition getIntv (val:ann) := fst val.
 Definition getErr (val:ann) := snd val.
 
 Arguments getIntv _ /.
-Arguments getErr _ /.
\ No newline at end of file
+Arguments getErr _ /.
+
+(**
+  Define environment update function for arbitrary type as abbreviation.
+  This must be defined here, to prove some lemmas about expression evaluation.
+  See below.
+**)
+Definition updEnv (x:nat) (v: R) (E: nat -> R) (y:nat) :R :=
+  if y =? x
+  then v
+  else E y.
\ No newline at end of file

coq/Infra/RealSimps.v

+Require Import Coq.Reals.Reals Coq.micromega.Psatz Coq.Setoids.Setoid.
+
+(**
+  TODO REORDER AND DOCUMENT
+**)
+
+Lemma Rsub_eq_Ropp_Rplus (a:R) (b:R) :
+  (a - b = a + (- b))%R.
+Proof.
+  field_simplify; reflexivity.
+Qed.
+
+Lemma Rabs_simplify (a b: R) :
+  (Rabs a <= b <-> ((a <= 0) /\ a >= - b) \/ (a >= 0 /\ a <= b))%R.
+Proof.
+  split.
+  - intros abs.
+    unfold Rabs in abs.
+    destruct Rcase_abs in abs; lra.
+  - intros cases_abs.
+    unfold Rabs.
+    destruct Rcase_abs; lra.
+Qed.
+
+
+Lemma bound_flip_sub:
+  forall a b e,
+    (Rabs (a - b) <= e ->
+     Rabs ( - a - (- b)) <= e)%R.
+Proof.
+  intros a b e abs_less.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite <- Ropp_plus_distr.
+  rewrite Rabs_Ropp.
+  rewrite <- Rsub_eq_Ropp_Rplus; auto.
+Qed.
+
+
+Lemma plus_bounds_simplify:
+  forall a b c d e, (a + b + ( c + d + e) = (a + c) + (b + d) + e)%R.
+Proof.
+  intros.
+  repeat rewrite <- Rplus_assoc.
+  setoid_rewrite Rplus_comm at 4.
+  setoid_rewrite Rplus_assoc at 3.
+  setoid_rewrite Rplus_comm at 3; auto.
+Qed.
+
+Lemma Rabs_err_simpl:
+  forall a b,
+    (Rabs (a - (a * (1 + b))) = Rabs (a * b))%R.
+Proof.
+  intros a b.
+  rewrite Rmult_plus_distr_l.
+  rewrite Rmult_1_r.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite Ropp_plus_distr.
+  assert (- a + (- (a * b)) = ( - (a * b) + - a))%R by (rewrite Rplus_comm; auto).
+  rewrite H.
+  rewrite <- Rsub_eq_Ropp_Rplus.
+  rewrite Rplus_minus.
+  rewrite Rabs_Ropp; reflexivity.
+Qed.
+
+(**
+  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.
\ No newline at end of file