Merge branch 'fma_proofs_merge' into 'certificates'

Fma proofs merge

See merge request AVA/daisy!170

coq/ErrorBounds.v

@@ -258,6 +258,75 @@ Proof.
   apply Rabs_pos.
 Qed.
 
+Lemma fma_abs_err_bounded (e1:exp Q) (e1R:R) (e1F:R) (e2:exp Q) (e2R:R) (e2F:R)
+      (e3:exp Q) (e3R:R) (e3F:R)
+      (vR:R) (vF:R) (E1 E2:env) (m m1 m2 m3:mType) defVars:
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e1)) e1R M0 ->
+  eval_exp E2 defVars (toRExp e1) e1F m1->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e2)) e2R M0 ->
+  eval_exp E2 defVars (toRExp e2) e2F m2 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e3)) e3R M0 ->
+  eval_exp E2 defVars (toRExp e3) e3F m3->
+  eval_exp E1 (toRMap defVars) (toREval (Fma (toRExp e1) (toRExp e2) (toRExp e3))) vR M0 ->
+  eval_exp (updEnv 3 e3F (updEnv 2 e2F (updEnv 1 e1F emptyEnv)))
+           (updDefVars 3 m3 (updDefVars 2 m2 (updDefVars 1 m1 defVars)))
+           (Fma (Var R 1) (Var R 2) (Var R 3)) vF m ->
+  (Rabs (vR - vF) <= Rabs ((e1R - e1F) + (e2R * e3R - e2F * e3F)) + Rabs (e1F + e2F * e3F) * (Q2R (mTypeToQ m)))%R.
+Proof.
+  intros e1_real e1_float e2_real e2_float e3_real e3_float fma_real fma_float.
+  inversion fma_real; subst;
+  assert (m0 = M0) by (eapply toRMap_eval_M0; eauto).
+  assert (m4 = M0) by (eapply toRMap_eval_M0; eauto).
+  assert (m5 = M0) by (eapply toRMap_eval_M0; eauto).
+  subst; simpl in H3; rewrite Q2R0_is_0 in H3; auto.
+  rewrite delta_0_deterministic in fma_real; auto.
+  rewrite delta_0_deterministic; auto.
+  unfold evalFma in *; simpl in *.
+  clear delta H3.
+  rewrite (meps_0_deterministic (toRExp e1) H5 e1_real);
+    rewrite (meps_0_deterministic (toRExp e2) H6 e2_real);
+    rewrite (meps_0_deterministic (toRExp e3) H7 e3_real).
+  rewrite (meps_0_deterministic (toRExp e1) H5 e1_real) in fma_real.
+  rewrite (meps_0_deterministic (toRExp e2) H6 e2_real) in fma_real.
+  rewrite (meps_0_deterministic (toRExp e3) H7 e3_real) in fma_real.
+  clear H5 H6 v1 v2 v3 H7 H2.
+  inversion fma_float; subst.
+  unfold evalFma in *.
+  unfold perturb; simpl.
+  inversion H3; subst; inversion H6; subst; inversion H7; subst.
+  unfold updEnv in *; simpl in *.
+  inversion H5; inversion H1; inversion H9; subst.
+  clear fma_float H7 fma_real e1_real e1_float e2_real e2_float e3_real e3_float H6 H1 H5 H9 H3 H0 H4 H8.
+  repeat rewrite Rmult_plus_distr_l.
+  rewrite Rmult_1_r.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite Ropp_plus_distr.
+  rewrite <- Rplus_assoc.
+  setoid_rewrite <- Rsub_eq_Ropp_Rplus at 2.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite Rsub_eq_Ropp_Rplus.
+  rewrite <- Rplus_assoc.
+  setoid_rewrite Rplus_comm at 8.
+  rewrite <- Rplus_assoc.
+  setoid_rewrite Rplus_comm at 9.
+  rewrite Rplus_assoc.
+  setoid_rewrite Rplus_assoc at 2.
+  rewrite <- Rplus_assoc.
+  rewrite <- Rsub_eq_Ropp_Rplus.
+  rewrite <- Rsub_eq_Ropp_Rplus.
+  rewrite <- Ropp_plus_distr.
+  rewrite <- Rsub_eq_Ropp_Rplus.
+  eapply Rle_trans.
+  eapply Rabs_triang.
+  eapply Rplus_le_compat_l.
+  rewrite Rabs_Ropp.
+  repeat rewrite Rabs_mult.
+  eapply Rmult_le_compat_l; auto.
+  apply Rabs_pos.
+Qed.
+  
+
 Lemma err_prop_inversion_pos_real nF nR err elo ehi
       (float_iv_pos : (0 < elo - err)%R)
       (real_iv_pos : (0 < elo)%R)

coq/ErrorValidation.v

@@ -69,6 +69,25 @@ Fixpoint validErrorbound (e:exp Q) (* analyzed expression *)
          | _, _ => false
          end
        else false
+      | Fma e1 e2 e3 =>
+        if ((validErrorbound e1 typeMap A dVars)
+              && (validErrorbound e2 typeMap A dVars)
+              && (validErrorbound e3 typeMap A dVars))
+        then
+          match DaisyMap.find e1 A, DaisyMap.find e2 A, DaisyMap.find e3 A with
+          | Some (ive1, err1), Some (ive2, err2), Some (ive3, err3) =>
+            let errIve1 := widenIntv ive1 err1 in
+            let errIve2 := widenIntv ive2 err2 in
+            let errIve3 := widenIntv ive3 err3 in
+            let upperBoundE1 := maxAbs ive1 in
+            let upperBoundE2 := maxAbs ive2 in
+            let upperBoundE3 := maxAbs ive3 in
+            let errIntv_prod := multIntv errIve2 errIve3 in
+            let mult_error_bound := (upperBoundE2 * err3 + upperBoundE3 * err2 + err2 * err3) in
+            Qleb (err1 + mult_error_bound + (maxAbs (addIntv errIve1 errIntv_prod)) * (mTypeToQ m)) err
+          | _, _, _ => false
+          end
+        else false
       |Downcast m1 e1 =>
        if validErrorbound e1 typeMap A dVars
        then
@@ -407,55 +426,19 @@ Proof.
       repeat rewrite Q2R_minus; lra.
 Qed.
 
-Lemma validErrorboundCorrectMult E1 E2 A
-      (e1:exp Q) (e2:exp Q) (nR nR1 nR2 nF nF1 nF2 :R) (e err1 err2 :error)
-      (alo ahi e1lo e1hi e2lo e2hi:Q) dVars (m m1 m2:mType) Gamma defVars:
-  m = join m1 m2 ->
-  eval_exp E1 (toRMap defVars) (toREval (toRExp e1)) nR1 M0 ->
-  eval_exp E1 (toRMap defVars) (toREval (toRExp e2)) nR2 M0 ->
-  eval_exp E1 (toRMap defVars) (toREval (toRExp (Binop Mult e1 e2))) nR M0 ->
-  eval_exp E2 defVars (toRExp e1) nF1 m1 ->
-  eval_exp E2 defVars (toRExp e2) nF2 m2 ->
-  eval_exp (updEnv 2 nF2 (updEnv 1 nF1 emptyEnv))
-           (updDefVars 2 m2 (updDefVars 1 m1 defVars))
-           (toRExp (Binop Mult (Var Q 1) (Var Q 2))) nF m ->
-  typeCheck (Binop Mult e1 e2) defVars Gamma = true ->
-  validErrorbound (Binop Mult e1 e2) Gamma A dVars = true ->
+Lemma multiplicationErrorBounded e1lo e1hi e2lo e2hi nR1 nF1 nR2 nF2 err1 err2 :
   (Q2R e1lo <= nR1 <= Q2R e1hi)%R ->
   (Q2R e2lo <= nR2 <= Q2R e2hi)%R ->
-  DaisyMap.find e1 A = Some ((e1lo,e1hi),err1) ->
-  DaisyMap.find e2 A = Some ((e2lo, e2hi),err2) ->
-  DaisyMap.find (Binop Mult e1 e2) A = Some ((alo,ahi),e)->
-  (Rabs (nR1 - nF1) <= (Q2R err1))%R ->
-  (Rabs (nR2 - nF2) <= (Q2R err2))%R ->
-  (Rabs (nR - nF) <= (Q2R e))%R.
+  (Rabs (nR1 - nF1) <= Q2R err1)%R ->
+  (Rabs (nR2 - nF2) <= Q2R err2)%R ->
+  (0 <= Q2R err1)%R ->
+  (0 <= Q2R err2)%R ->
+  (Rabs (nR1 * nR2 - nF1 * nF2) <=
+   RmaxAbsFun (Q2R e1lo, Q2R e1hi) * Q2R err2 + RmaxAbsFun (Q2R e2lo, Q2R e2hi) * Q2R err1 +
+   Q2R err1 * Q2R err2)%R.
 Proof.
-  intros mIsJoin e1_real e2_real eval_real e1_float e2_float eval_float
-         subexpr_ok valid_error valid_e1 valid_e2 A_e1 A_e2 A_mult
-         err1_bounded err2_bounded.
-  cbn in *; Daisy_compute; type_conv; subst.
-  eapply Rle_trans.
-  eapply (mult_abs_err_bounded e1 e2); eauto.
-  pose proof (typingSoundnessExp _ _ R2 e1_float).
-  pose proof (typingSoundnessExp _ _ R1 e2_float).
-  rewrite H in Heqo0; rewrite H0 in Heqo1; inversion Heqo0; inversion Heqo1; subst.
-  clear H H0.
-  rename R0 into valid_error.
-  assert (0 <= Q2R err1)%R as err1_pos.
-  { pose proof (err_always_positive e1 Gamma A dVars);
-      eauto. }
-  assert (0 <= Q2R err2)%R as err2_pos.
-  { pose proof (err_always_positive e2 Gamma A dVars); eauto. }
-  clear R2 R1.
-  canonize_hyps.
-  repeat rewrite Q2R_plus in valid_error.
-  repeat rewrite Q2R_mult in valid_error.
-  repeat rewrite <- maxAbs_impl_RmaxAbs in valid_error.
-  eapply Rle_trans.
-  Focus 2.
-  apply valid_error.
-  apply Rplus_le_compat.
-  - unfold Rabs in err1_bounded.
+  intros valid_e1 valid_e2 err1_bounded err2_bounded err1_pos err2_pos.
+  unfold Rabs in err1_bounded.
   unfold Rabs in err2_bounded.
   (* Before doing case distinction, prove bounds that will be used many times: *)
   assert (nR1 <= RmaxAbsFun (Q2R e1lo, Q2R e1hi))%R
@@ -486,8 +469,6 @@ Proof.
     as nR1_to_sum by lra.
   assert (RmaxAbsFun (Q2R e1lo, Q2R e1hi) * Q2R err2 + RmaxAbsFun (Q2R e2lo, Q2R e2hi) * Q2R err1 <=  RmaxAbsFun (Q2R e1lo, Q2R e1hi) * Q2R err2 + RmaxAbsFun (Q2R e2lo, Q2R e2hi) * Q2R err1 + Q2R err1 * Q2R err2)%R
     as sum_to_errsum by lra.
-    clear e1_real e1_float e2_real e2_float eval_real eval_float valid_error
-      A_e1 A_e2.
   (* Large case distinction for
          a) different cases of the value of Rabs (...) and
          b) wether arguments of multiplication in (nf1 * nF2) are < or >= 0 *)
@@ -909,6 +890,57 @@ Proof.
           apply H.
           lra.
       }
+Qed.
+
+Lemma validErrorboundCorrectMult E1 E2 A
+      (e1:exp Q) (e2:exp Q) (nR nR1 nR2 nF nF1 nF2 :R) (e err1 err2 :error)
+      (alo ahi e1lo e1hi e2lo e2hi:Q) dVars (m m1 m2:mType) Gamma defVars:
+  m = join m1 m2 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e1)) nR1 M0 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e2)) nR2 M0 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp (Binop Mult e1 e2))) nR M0 ->
+  eval_exp E2 defVars (toRExp e1) nF1 m1 ->
+  eval_exp E2 defVars (toRExp e2) nF2 m2 ->
+  eval_exp (updEnv 2 nF2 (updEnv 1 nF1 emptyEnv))
+           (updDefVars 2 m2 (updDefVars 1 m1 defVars))
+           (toRExp (Binop Mult (Var Q 1) (Var Q 2))) nF m ->
+  typeCheck (Binop Mult e1 e2) defVars Gamma = true ->
+  validErrorbound (Binop Mult e1 e2) Gamma A dVars = true ->
+  (Q2R e1lo <= nR1 <= Q2R e1hi)%R ->
+  (Q2R e2lo <= nR2 <= Q2R e2hi)%R ->
+  DaisyMap.find e1 A = Some ((e1lo,e1hi),err1) ->
+  DaisyMap.find e2 A = Some ((e2lo, e2hi),err2) ->
+  DaisyMap.find (Binop Mult e1 e2) A = Some ((alo,ahi),e)->
+  (Rabs (nR1 - nF1) <= (Q2R err1))%R ->
+  (Rabs (nR2 - nF2) <= (Q2R err2))%R ->
+  (Rabs (nR - nF) <= (Q2R e))%R.
+Proof.
+  intros mIsJoin e1_real e2_real eval_real e1_float e2_float eval_float
+         subexpr_ok valid_error valid_e1 valid_e2 A_e1 A_e2 A_mult
+         err1_bounded err2_bounded.
+  cbn in *; Daisy_compute; type_conv; subst.
+  eapply Rle_trans.
+  eapply (mult_abs_err_bounded e1 e2); eauto.
+  pose proof (typingSoundnessExp _ _ R2 e1_float).
+  pose proof (typingSoundnessExp _ _ R1 e2_float).
+  rewrite H in Heqo0; rewrite H0 in Heqo1; inversion Heqo0; inversion Heqo1; subst.
+  clear H H0.
+  rename R0 into valid_error.
+  assert (0 <= Q2R err1)%R as err1_pos.
+  { pose proof (err_always_positive e1 Gamma A dVars);
+      eauto. }
+  assert (0 <= Q2R err2)%R as err2_pos.
+  { pose proof (err_always_positive e2 Gamma A dVars); eauto. }
+  clear R2 R1.
+  canonize_hyps.
+  repeat rewrite Q2R_plus in valid_error.
+  repeat rewrite Q2R_mult in valid_error.
+  repeat rewrite <- maxAbs_impl_RmaxAbs in valid_error.
+  eapply Rle_trans.
+  Focus 2.
+  apply valid_error.
+  apply Rplus_le_compat.
+  - eauto using multiplicationErrorBounded.
   - remember (multIntv (widenIntv (e1lo, e1hi) err1) (widenIntv (e2lo, e2hi) err2)) as iv.
     iv_assert iv iv_unf.
     destruct iv_unf as [ivl [ivh iv_unf]].
@@ -1862,6 +1894,108 @@ Proof.
             rewrite (Qmult_inj_r) in H3; auto. }
 Qed.
 
+Lemma validErrorboundCorrectFma E1 E2 A
+      (e1:exp Q) (e2:exp Q) (e3: exp Q) (nR nR1 nR2 nR3 nF nF1 nF2 nF3: R) (e err1 err2 err3 :error)
+      (alo ahi e1lo e1hi e2lo e2hi e3lo e3hi:Q) dVars (m m1 m2 m3:mType) Gamma defVars:
+  m = join3 m1 m2 m3 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e1)) nR1 M0 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e2)) nR2 M0 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp e3)) nR3 M0 ->
+  eval_exp E1 (toRMap defVars) (toREval (toRExp (Fma e1 e2 e3))) nR M0 ->
+  eval_exp E2 defVars (toRExp e1) nF1 m1 ->
+  eval_exp E2 defVars (toRExp e2) nF2 m2 ->
+  eval_exp E2 defVars (toRExp e3) nF3 m3 ->
+  eval_exp (updEnv 3 nF3 (updEnv 2 nF2 (updEnv 1 nF1 emptyEnv)))
+           (updDefVars 3 m3 (updDefVars 2 m2 (updDefVars 1 m1 defVars)))
+           (toRExp (Fma (Var Q 1) (Var Q 2) (Var Q 3))) nF m ->
+  typeCheck (Fma e1 e2 e3) defVars Gamma = true ->
+  validErrorbound (Fma e1 e2 e3) Gamma A dVars = true ->
+  (Q2R e1lo <= nR1 <= Q2R e1hi)%R ->
+  (Q2R e2lo <= nR2 <= Q2R e2hi)%R ->
+  (Q2R e3lo <= nR3 <= Q2R e3hi)%R ->
+  DaisyMap.find e1 A = Some ((e1lo,e1hi),err1) ->
+  DaisyMap.find e2 A = Some ((e2lo, e2hi),err2) ->
+  DaisyMap.find e3 A = Some ((e3lo, e3hi),err3) ->
+  DaisyMap.find (Fma e1 e2 e3) A = Some ((alo,ahi),e)->
+  (Rabs (nR1 - nF1) <= (Q2R err1))%R ->
+  (Rabs (nR2 - nF2) <= (Q2R err2))%R ->
+  (Rabs (nR3 - nF3) <= (Q2R err3))%R ->
+  (Rabs (nR - nF) <= (Q2R e))%R.
+Proof.
+  intros mIsJoin e1_real e2_real e3_real eval_real e1_float e2_float e3_float eval_float
+         subexpr_ok valid_error valid_e1 valid_e2 valid_e3 A_e1 A_e2 A_e3 A_fma
+         err1_bounded err2_bounded err3_bounded.
+  cbn in *; Daisy_compute; type_conv; subst.
+  eapply Rle_trans.
+  eapply (fma_abs_err_bounded e1 e2 e3); eauto.
+  pose proof (typingSoundnessExp _ _ R4 e1_float).
+  pose proof (typingSoundnessExp _ _ R3 e2_float).
+  pose proof (typingSoundnessExp _ _ R2 e3_float).
+  rewrite H in Heqo0; rewrite H0 in Heqo1; rewrite H1 in Heqo2;
+    inversion Heqo0; inversion Heqo1; inversion Heqo2; subst.
+  rename R0 into valid_error.
+  assert (0 <= Q2R err1)%R as err1_pos by (eapply (err_always_positive e1 Gamma A dVars); eauto).
+  assert (0 <= Q2R err2)%R as err2_pos by (eapply (err_always_positive e2 Gamma A dVars); eauto).
+  assert (0 <= Q2R err3)%R as err3_pos by (eapply (err_always_positive e3 Gamma A dVars); eauto).
+  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; eauto.
+  apply Rplus_le_compat.
+  - eauto using Rle_trans, Rabs_triang, Rplus_le_compat, multiplicationErrorBounded.
+  - apply Rmult_le_compat_r; auto using mTypeToQ_pos_R.
+    remember (multIntv (widenIntv (e2lo, e2hi) err2) (widenIntv (e3lo, e3hi) err3)) as iv_prod.
+    remember (addIntv (widenIntv (e1lo, e1hi) err1) iv_prod) as iv_sum.
+    iv_assert iv_sum iv_unf.
+    destruct iv_unf as [ivl [ivh iv_unf]].
+    rewrite iv_unf.
+    rewrite <- maxAbs_impl_RmaxAbs.
+    assert (ivlo iv_sum = ivl) by (rewrite iv_unf; auto).
+    assert (ivhi iv_sum = ivh) by (rewrite iv_unf; auto).
+    (* rewrite <- H, <- H0. *)
+    assert (contained nR1 (Q2R e1lo, Q2R e1hi)) as contained_intv1 by auto.
+    pose proof (distance_gives_iv (a:=nR1) _ (Q2R e1lo, Q2R e1hi) contained_intv1 err1_bounded).
+    assert (contained nR2 (Q2R e2lo, Q2R e2hi)) as contained_intv2 by auto.
+    pose proof (distance_gives_iv (a:=nR2) _ _ contained_intv2 err2_bounded).
+    assert (contained nR3 (Q2R e3lo, Q2R e3hi)) as contained_intv3 by auto.
+    pose proof (distance_gives_iv (a:=nR3) _ _ contained_intv3 err3_bounded).
+    pose proof (IntervalArith.interval_multiplication_valid _ _ H5 H6).
+    pose proof (IntervalArith.interval_addition_valid _ _ H4 H7).
+    destruct H8.
+    unfold RmaxAbsFun.
+    subst.
+    apply RmaxAbs; subst; simpl in *.
+    unfold RmaxAbsFun.
+    + rewrite Q2R_min4.
+      repeat rewrite Q2R_mult;
+        repeat rewrite Q2R_minus;
+        repeat rewrite Q2R_plus;
+        repeat rewrite Q2R_minus.
+      rewrite Q2R_max4.
+      rewrite Q2R_min4.
+      repeat rewrite Q2R_mult;
+        repeat rewrite Q2R_minus;
+        repeat rewrite Q2R_plus;
+        repeat rewrite Q2R_minus.
+      assumption.
+    + rewrite Q2R_max4.
+      repeat rewrite Q2R_mult;
+        repeat rewrite Q2R_minus;
+        repeat rewrite Q2R_plus;
+        repeat rewrite Q2R_minus.
+      rewrite Q2R_max4.
+      rewrite Q2R_min4.
+      repeat rewrite Q2R_mult;
+        repeat rewrite Q2R_minus;
+        repeat rewrite Q2R_plus;
+        repeat rewrite Q2R_minus.
+      assumption.
+Qed.
 
 Lemma validErrorboundCorrectRounding E1 E2 A (e: exp Q) (nR nF nF1: R) (err err':error) (elo ehi alo ahi: Q) dVars (m: mType) (machineEpsilon:mType) Gamma defVars:
   eval_exp E1 (toRMap defVars) (toREval (toRExp e)) nR M0 ->
@@ -2038,6 +2172,80 @@ Proof.
         { cbn; instantiate (1:=dVars); Daisy_compute.
           rewrite L, L2,L4, R1; simpl; auto. }
         { andb_to_prop R; auto. }
+  (*- simpl in valid_error.
+    destruct (absenv e1) as [[ivlo1 ivhi1] err1] eqn:absenv_e1;
+      destruct (absenv e2) as [[ivlo2 ivhi2] err2] eqn:absenv_e2;
+      destruct (absenv e3) as [[ivlo3 ivhi3] err3] eqn:absenv_e3.
+    subst; simpl in *.
+    rewrite absenv_eq, absenv_e1, absenv_e2, absenv_e3 in *.
+    case_eq (Gamma (Fma e1 e2 e3));
+      intros * type_fma; rewrite type_fma in *; [ | inversion valid_error ].
+    case_eq (Gamma e1);
+      intros * type_e1; rewrite type_e1 in typing_ok; [ | inversion typing_ok ].
+    case_eq (Gamma e2);
+      intros * type_e2; rewrite type_e2 in typing_ok; [ | inversion typing_ok ].
+    case_eq (Gamma e3);
+      intros * type_e3; rewrite type_e3 in typing_ok; [ | inversion typing_ok ].
+    repeat match goal with
+           | [H: _ = true |- _] => andb_to_prop H
+           end.
+    type_conv.*)
+  - cbn in *. rewrite A_eq in *.
+    Daisy_compute; try congruence; type_conv; subst; simpl in *.
+    inversion eval_real; subst.
+    assert (m0 = M0 /\ m4 = M0 /\ m5 = M0) as [? [? ?]] by (split; try split; eapply toRMap_eval_M0; eauto); subst.
+    destruct i as [ivlo1 ivhi1]; destruct i2 as [ivlo2 ivhi2]; destruct i1 as [ivlo3 ivhi3];
+      rename e into err1; rename e5 into err2; rename e4 into err3.
+    destruct (IHe1 E1 E2 fVars dVars A v1 err1 P ivlo1 ivhi1 Gamma defVars)
+      as [[vF1 [mF1 eval_float_e1]] bounded_e1];
+      try auto; set_tac.
+    destruct (IHe2 E1 E2 fVars dVars A v2 err2 P ivlo2 ivhi2 Gamma defVars)
+      as [[vF2 [mF2 eval_float_e2]] bounded_e2];
+      try auto; set_tac.
+    destruct (IHe3 E1 E2 fVars dVars A v3 err3 P ivlo3 ivhi3 Gamma defVars)
+      as [[vF3 [mF3 eval_float_e3]] bounded_e3];
+      try auto; set_tac.
+    destruct (validIntervalbounds_sound _ (E:=E1) (Gamma:=defVars) L (fVars := fVars) (dVars:=dVars))
+      as [iv1' [ err1' [v1' [map_e1 [eval_real_e1 bounds_e1]]]]];
+      try auto; set_tac.
+    rewrite map_e1 in Heqo; inversion Heqo; subst.
+    pose proof (meps_0_deterministic _ eval_real_e1 H5); subst; clear H5.
+    destruct (validIntervalbounds_sound _ (E:=E1) (Gamma:=defVars) R1 (fVars:= fVars) (dVars:=dVars))
+      as [iv2' [err2' [v2' [map_e2 [eval_real_e2 bounds_e2]]]]];
+      try auto; set_tac.
+    rewrite map_e2 in Heqo2; inversion Heqo2; subst.
+    pose proof (meps_0_deterministic _ eval_real_e2 H6); subst; clear H6.
+    destruct (validIntervalbounds_sound _ (E:=E1) (Gamma:=defVars) R0 (fVars:= fVars) (dVars:=dVars))
+      as [iv3' [err3' [v3' [map_e3 [eval_real_e3 bounds_e3]]]]];
+      try auto; set_tac.
+    rewrite map_e3 in Heqo5; inversion Heqo5; subst.
+    pose proof (meps_0_deterministic _ eval_real_e3 H7); subst; clear H7.
+    split.
+    + repeat eexists; econstructor; eauto.
+      rewrite Rabs_right; try lra.
+      instantiate (1 := 0%R).
+      apply mTypeToQ_pos_R.
+      apply Rle_ge.
+      hnf; right; reflexivity.
+    + intros * eval_float.
+      clear eval_float_e1 eval_float_e2 eval_float_e3.
+      inversion eval_float; subst.
+      eapply (fma_unfolding H4 H5 H8 H9) in eval_float; try auto.
+      eapply (validErrorboundCorrectFma (e1:=e1) (e2:=e2) (e3:=e3) A); eauto.
+      { simpl.
+        rewrite Heqo0.
+        rewrite Heqo4.
+        rewrite Heqo6.
+        rewrite Heqo7.
+        rewrite mTypeEq_refl, R5, R6, R7; auto. }
+      { simpl.
+        rewrite A_eq.
+        rewrite Heqo0.
+        rewrite R, L1, L2, R4; simpl.
+        rewrite map_e1, map_e2, map_e3.
+        inversion Heqo1.
+        rewrite <- H0.
+        auto. }
   - cbn in *; Daisy_compute; try congruence; type_conv; subst.
     inversion eval_real; subst.
     apply M0_least_precision in H1.

coq/Expressions.v

@@ -98,6 +98,9 @@ Definition evalUnop (o:unop) (v:R):=
   |Inv => (/ v)%R
   end .
 
+Definition evalFma (v1:R) (v2:R) (v3:R):=
+  evalBinop Plus v1 (evalBinop Mult v2 v3).
+
 (**
   Define expressions parametric over some value type V.
   Will ease reasoning about different instantiations later.
@@ -107,6 +110,7 @@ Inductive exp (V:Type): Type :=
 | Const: mType -> V -> exp V
 | Unop: unop -> exp V -> exp V
 | Binop: binop -> exp V -> exp V -> exp V
+| Fma: exp V -> exp V -> exp V -> exp V
 | Downcast: mType -> exp V -> exp V.
 
 (**
@@ -122,6 +126,8 @@ Fixpoint expEq (e1:exp Q) (e2:exp Q) :=
     (unopEq o1 o2) && (expEq e11 e22)
   | Binop o1 e11 e12, Binop o2 e21 e22 =>
     (binopEq o1 o2) && (expEq e11 e21) && (expEq e12 e22)
+  | Fma e11 e12 e13, Fma e21 e22 e23 =>
+    (expEq e11 e21) && (expEq e12 e22) && (expEq e13 e23)
   | Downcast m1 f1, Downcast m2 f2 =>
     (mTypeEq m1 m2) && (expEq f1 f2)
   | _, _ => false
@@ -136,6 +142,7 @@ Proof.
   - apply Qeq_bool_iff; lra.
   - case u; auto.
   - case b; auto.
+  - firstorder.
   - apply mTypeEq_refl.
 Qed.
 
@@ -156,6 +163,9 @@ Proof.
       * destruct b; auto.
       * apply IHe1.
     + apply IHe2.
+  - f_equal.
+    + f_equal; auto.
+    + auto.
   - f_equal.
     + apply mTypeEq_sym; auto.
     + apply IHe.
@@ -180,6 +190,9 @@ Proof.
     rewrite binopEq_refl; simpl.
     apply andb_true_iff.
     split; [eapply IHe1; eauto | eapply IHe2; eauto].
+  - rewrite andb_true_iff.
+    rewrite andb_true_iff.
+    split; [split; [eapply IHe1; eauto | eapply IHe2; eauto] | eapply IHe3; eauto].
   - rewrite mTypeEq_refl; simpl.
     eapply IHe; eauto.
 Qed.
@@ -206,6 +219,7 @@ Module ExpOrderedType (V_ordered:OrderType) <: OrderType.
       if unopEq u1 u2
       then expCompare e1 e2
       else (if unopEq u1 Neg then Lt else Gt)
+    | Unop _ _, Fma _ _ _ => Lt
     | Unop _ _, Binop _ _ _ => Lt
     | Unop _ _, Downcast _ _ => Lt
     | Unop _ _, _ => Gt
@@ -213,6 +227,7 @@ Module ExpOrderedType (V_ordered:OrderType) <: OrderType.
       if mTypeEq m1 m2
       then expCompare e1 e2
       else (if morePrecise m1 m2 then Lt else Gt)
+    | Downcast _ _, Fma _ _ _ => Lt
     | Downcast _ _, Binop _ _ _ => Lt
     | Downcast _ _, _ => Gt
     | Binop b1 e11 e12, Binop b2 e21 e22 =>
@@ -238,7 +253,19 @@ Module ExpOrderedType (V_ordered:OrderType) <: OrderType.
         end
       | _ => res
       end
-    |_ , _ => Gt
+    | Binop _ _ _, Fma _ _ _ => Lt
+    | Binop _ _ _, _ => Gt
+    | Fma e11 e12 e13, Fma e21 e22 e23 =>
+      match expCompare e11 e21 with
+      | Eq => match expCompare e12 e22 with
+             | Eq => expCompare e13 e23
+             | Lt => Lt
+             | Gt => Gt
+             end
+      | Lt => Lt
+      | Gt => Gt
+      end
+    | Fma _ _ _ , _ => Gt
     end.
 
   Lemma expCompare_refl e: expCompare e e = Eq.
@@ -248,6 +275,7 @@ Module ExpOrderedType (V_ordered:OrderType) <: OrderType.
     - rewrite mTypeEq_refl. apply V_orderedFacts.compare_refl.
     - rewrite unopEq_refl; auto.
     - rewrite IHe1, IHe2. destruct b; auto.
+    - now rewrite IHe1, IHe2, IHe3.
     - rewrite mTypeEq_refl; auto.
   Qed.