Prove error for multiplication, still needs testing wether computed bound is...

Prove error for multiplication, still needs testing wether computed bound is low enough for actual validation

coq/ErrorBounds.v

@@ -163,4 +163,88 @@ Proof.
        apply Rabs_triang.
        rewrite Rabs_Ropp.
        apply Req_le; auto.
-Qed.
\ No newline at end of file
+Qed.
+
+Lemma mult_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 Mult e1 e2) vR ->
+  eval_exp machineEpsilon (updEnv 2 e2F (updEnv 1 e1F cenv)) (Binop Mult (Var R 1) (Var R 2)) vF ->
+  (Rabs (e1R - e1F) <= err1)%R ->
+  (Rabs (e2R - e2F) <= err2)%R ->
+  (Rabs (vR - vF) <= Rabs e1R * Rabs e2R + (Rabs e1F * Rabs e2F + Rabs e1F * Rabs e2F * machineEpsilon))%R.
+Proof.
+  intros e1_real e1_float e2_real e2_float mult_real mult_float bound_e1 bound_e2.
+  (* Prove that e1R and e2R are the correct values and that vR is e1R + e2R *)
+  inversion mult_real; subst.
+  rewrite perturb_0_val in mult_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 mult_real.
+  rewrite (eval_det H5 e2_real) in mult_real.
+  clear H4 H5 v1 v2.
+  (* Now unfold the float valued evaluation to get the deltas we need for the inequality *)
+  inversion mult_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 mult_float H4 H5 mult_real e1_real e1_float e2_real e2_float.
+  repeat rewrite Rmult_plus_distr_l.
+  rewrite Rmult_1_r.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite Ropp_plus_distr.
+  eapply Rle_trans.
+  eapply Rabs_triang.
+  eapply Rle_trans.
+  eapply Rplus_le_compat_l.
+  eapply Rabs_triang.
+  repeat rewrite Rabs_Ropp.
+  repeat rewrite Rabs_mult.
+  eapply Rle_trans.
+  eapply Rplus_le_compat_l.
+  eapply Rplus_le_compat_l.
+  eapply Rmult_le_compat_l.
+  rewrite <- Rabs_mult.
+  apply Rabs_pos.
+  apply H2.
+  apply Req_le; auto.
+Qed.
+(*
+  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 <- Rmult_eq_Ropp_Rplus.
+  rewrite Rabs_Ropp.
+  assert (Rabs (- e2R - - e2F)%R = Rabs (e2R - e2F)).
+  - rewrite Rmult_eq_Ropp_Rplus.
+    rewrite <- Ropp_plus_distr.
+    rewrite Rabs_Ropp.
+    rewrite <- Rmult_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.
+       rewrite Rmult_eq_Ropp_Rplus.
+       eapply Rle_trans.
+       eapply Rmult_le_compat_r.
+       unfold machineEpsilon, RealConstruction.realFromNum, RealConstruction.negativePower; interval.
+       apply Rabs_triang.
+       rewrite Rabs_Ropp.
+       apply Req_le; auto.
+*)
\ No newline at end of file

coq/ErrorValidation.v

@@ -34,6 +34,19 @@ Section ComputableErrors.
     unfold Qle_bool; auto.
   Qed.
 
+(*
+Multiplication :
+  let eterm := (Rabs (e1R * err2) + Rabs(e2R * err1) + Rabs (err1 * err2) +
+               Rabs (e1R * e2R * machineEpsilon) + Rabs (e1R * err1 * machineEpsilon) +
+               Rabs (e2R * err2 * machineEpsilon) + (err1 * err2 * machineEpsilon))%R
+ WAS:
+Qleb ((upperBoundE1 * upperBoundE2) + (e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2) + (((e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2))*machineEpsilon)) err
+WAS:
+Qabs (upperBoundE1 * err2) + Qabs(upperBoundE2 * err1) + Qabs (err1 * err2) +
+                           Qabs (upperBoundE1 * upperBoundE2 * machineEpsilon) + Qabs (upperBoundE1 * err1 * machineEpsilon) +
+                           Qabs (upperBoundE2 * err2 * machineEpsilon) + (err1 * err2 * machineEpsilon)
+*)
+
   Fixpoint validErrorbound (e:exp Q) (env:analysisResult) :=
     let (intv, err) := (env e) in
       match e with
@@ -57,7 +70,7 @@ Section ComputableErrors.
             | Plus => Qleb (err1 + err2 + (Qabs e1F + Qabs e2F) * machineEpsilon) err
             (* TODO:Validity of next two computations *)
             | Sub => Qleb (err1 + err2 + ((Qabs e1F + Qabs e2F) * machineEpsilon)) err
-            | Mult => Qleb ((upperBoundE1 * upperBoundE2) + (e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2) + (((e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2))*machineEpsilon)) err
+            | Mult => Qleb (Qabs upperBoundE1 * Qabs upperBoundE2 + (Qabs e1F * Qabs e2F + Qabs e1F * Qabs e2F * machineEpsilon)) err
             | Div => false
             end
               in andb rec theVal
@@ -397,6 +410,110 @@ Section ComputableErrors.
           apply Rle_abs.
   Qed.
 
+  Lemma validErrorboundCorrectMult cenv absenv (e1:exp Q) (e2:exp Q) (nR nR1 nR2 nF nF1 nF2 :R) (e err1 err2 :error) (alo ahi e1lo e1hi e2lo e2hi:Q) P :
+        envApproximatesPrecond P absenv ->
+        precondValidForExec P cenv ->
+        eval_exp 0%R cenv (toRExp e1) nR1 ->
+        eval_exp 0%R cenv (toRExp e2) nR2 ->
+        eval_exp 0%R cenv (toRExp (Binop Mult e1 e2)) nR ->
+        eval_exp (Q2R machineEpsilon) cenv (toRExp e1) nF1 ->
+        eval_exp (Q2R machineEpsilon) cenv (toRExp e2) nF2 ->
+        eval_exp (Q2R machineEpsilon) (updEnv 2 nF2 (updEnv 1 nF1 cenv)) (toRExp (Binop Mult (Var Q 1) (Var Q 2))) nF ->
+        validErrorbound (Binop Mult e1 e2) absenv = true ->
+        validIntervalbounds (Binop Mult e1 e2) absenv P = true ->
+        absenv e1 = ((e1lo,e1hi),err1) ->
+        absenv e2 = ((e2lo, e2hi),err2) ->
+        absenv (Binop Mult e1 e2) = ((alo,ahi),e)->
+        (Rabs (nR1 - nF1) <= (Q2R err1))%R ->
+        (Rabs (nR2 - nF2) <= (Q2R err2))%R ->
+        (Rabs (nR - nF) <= (Q2R e))%R.
+  Proof.
+    intros env_approx_p p_valid e1_real e2_real eval_real e1_float e2_float eval_float
+           valid_error valid_intv absenv_e1 absenv_e2 absenv_mult
+           err1_bounded err2_bounded.
+    eapply Rle_trans.
+    eapply (mult_abs_err_bounded (toRExp e1) _ _ (toRExp e2)); eauto.
+    rewrite <- mEps_eq_Rmeps.
+    apply e1_float.
+    rewrite <- mEps_eq_Rmeps.
+    apply e2_float.
+    rewrite <- mEps_eq_Rmeps.
+    apply eval_float.
+    rewrite <- mEps_eq_Rmeps.
+    unfold validErrorbound in valid_error at 1.
+    rewrite absenv_mult, absenv_e1, absenv_e2 in valid_error.
+    apply Is_true_eq_left in valid_error.
+    apply andb_prop_elim in valid_error.
+    destruct valid_error as [valid_rec valid_error].
+    apply andb_prop_elim in valid_rec.
+    (** FIXME: Remove if trouble ahead *)
+    clear valid_rec.
+    apply Is_true_eq_true in valid_error.
+    apply Qle_bool_iff in valid_error.
+    apply Qle_Rle in valid_error.
+    repeat rewrite Q2R_plus in valid_error.
+    repeat rewrite Q2R_mult in valid_error.
+    repeat rewrite Q2R_plus in valid_error.
+    repeat rewrite <- Rabs_eq_Qabs in valid_error.
+    repeat rewrite Q2R_plus in valid_error.
+    repeat rewrite <- maxAbs_impl_RmaxAbs in valid_error.
+    eapply Rle_trans.
+    Focus 2.
+    apply valid_error.
+    (* Simplify Interval correctness assumption *)
+    simpl in valid_intv.
+    rewrite absenv_mult, absenv_e1, absenv_e2 in valid_intv.
+    apply Is_true_eq_left in valid_intv.
+    apply andb_prop_elim in valid_intv.
+    destruct valid_intv as [valid_rec valid_intv].
+    apply andb_prop_elim in valid_rec.
+    destruct valid_rec as [valid_iv_e1 valid_iv_e2].
+    apply Is_true_eq_true in valid_iv_e1; apply Is_true_eq_true in valid_iv_e2.
+    assert (Rabs nR1 <= RmaxAbsFun (e1lo,e1hi))%R.
+    - pose proof (validIntervalbounds_correct _ _ _ _ _ env_approx_p p_valid valid_iv_e1 e1_real).
+      unfold RmaxAbsFun.
+      rewrite absenv_e1 in H; simpl in H; destruct H.
+      eapply RmaxAbs; auto.
+    - assert (Rabs nR2 <= RmaxAbsFun (e2lo, e2hi))%R.
+      + pose proof (validIntervalbounds_correct _ _ _ _ _ env_approx_p p_valid valid_iv_e2 e2_real).
+        unfold RmaxAbsFun.
+        rewrite absenv_e2 in H0; simpl in H0; destruct H0.
+        eapply RmaxAbs; auto.
+      + rewrite Rabs_minus_sym in err1_bounded; rewrite Rabs_minus_sym in err2_bounded.
+        assert (Rabs nF1 - Rabs nR1 <= Q2R err1)%R.
+        * eapply Rle_trans; [apply Rabs_triang_inv | auto].
+        * assert (Rabs nF1 <= Rabs nR1 + Q2R err1)%R by lra.
+          assert (Rabs nF1 <= RmaxAbsFun (e1lo, e1hi) + Q2R err1)%R by (eapply Rle_trans; [ apply H2 | eapply Rplus_le_compat_r; auto]).
+          assert (Rabs nF1 <= Rabs (RmaxAbsFun (e1lo, e1hi) + Q2R err1))%R by (eapply Rle_trans; [apply H3 | apply Rle_abs]).
+          clear H1 H2 H3.
+          assert (Rabs nF2 - Rabs nR2 <= Q2R err2)%R by (eapply Rle_trans; [apply Rabs_triang_inv | auto]).
+          assert (Rabs nF2 <= Rabs nR2 + Q2R err2)%R by lra.
+          assert (Rabs nF2 <= RmaxAbsFun (e2lo, e2hi) + Q2R err2)%R by (eapply Rle_trans; [ apply H2 | eapply Rplus_le_compat_r; auto]).
+          assert (Rabs nF2 <= Rabs (RmaxAbsFun (e2lo, e2hi) + Q2R err2))%R by (eapply Rle_trans; [apply H3 | apply Rle_abs]).
+          clear H1 H2 H3.
+          rewrite <- Rplus_assoc.
+          eapply Rle_trans.
+          eapply Rplus_le_compat.
+          eapply Rplus_le_compat.
+          eapply Rmult_le_compat.
+          apply Rabs_pos.
+          apply Rabs_pos.
+          eapply Rle_trans.
+          apply H. apply Rle_abs.
+          eapply Rle_trans.
+          apply H0. apply Rle_abs.
+          apply Rmult_le_compat; [apply Rabs_pos | apply Rabs_pos | apply H4 | apply H5].
+          apply Rmult_le_compat.
+          rewrite <- Rabs_mult; apply Rabs_pos.
+          apply mEps_geq_zero.
+          apply Rmult_le_compat; [apply Rabs_pos | apply Rabs_pos | apply H4 | apply H5].
+          apply Req_le; auto.
+          apply Req_le.
+          unfold RmaxAbsFun.
+          simpl.
+          rewrite Rplus_assoc; auto.
+  Qed.
+
   Ltac iv_assert iv name:=
     assert (exists ivlo ivhi, iv = (ivlo, ivhi)) as name by (destruct iv; repeat eexists; auto).
 
@@ -481,8 +598,22 @@ Section ComputableErrors.
           { apply andb_prop_intro.
             split; apply Is_true_eq_left; auto. }
           { apply Is_true_eq_left; auto. }
-      + admit.
+      + eapply (validErrorboundCorrectMult cenv absenv e1 e2); eauto.
+        simpl.
+        rewrite absenv_eq, absenv_e1, absenv_e2; simpl.
+        apply Is_true_eq_true.
+        apply andb_prop_intro; split.
+        * apply andb_prop_intro.
+          split; apply Is_true_eq_left; auto.
+        * apply Is_true_eq_left; auto.
+        * unfold validIntervalbounds.
+          rewrite absenv_eq, absenv_e1, absenv_e2; simpl.
+          apply Is_true_eq_true.
+          apply andb_prop_intro; split.
+          { apply andb_prop_intro.
+            split; apply Is_true_eq_left; auto. }
+          { apply Is_true_eq_left; auto. }
       + inversion valid_error.
-  Admitted.
+  Qed.
 
 End ComputableErrors.

coq/Infra/RealSimps.v

@@ -81,4 +81,10 @@ Qed.
   Define the machine epsilon for floating point operations.
   Current value: 2^(-53)
  **)
-Definition machineEpsilon:R := realFromNum 1 1 53.
\ No newline at end of file
+Definition machineEpsilon:R := realFromNum 1 1 53.
+
+Lemma RmaxAbs_peel_Rabs a b:
+  Rmax (Rabs a) (Rabs b) = Rabs (Rmax (Rabs a) (Rabs b)).
+Proof.
+  apply Rmax_case_strong; intros; rewrite Rabs_Rabsolu; auto.
+Qed.