Prove error bounds for FMA

hol4/ErrorBoundsScript.sml

@@ -201,6 +201,52 @@ val div_abs_err_bounded = store_thm ("div_abs_err_bounded",
       \\ once_rewrite_tac[REAL_ABS_MUL]
       \\ match_mp_tac REAL_LE_MUL2 \\ fs[REAL_ABS_POS]));
 
+val fma_abs_err_bounded = store_thm ("fma_abs_err_bounded",
+  ``!(e1:real exp) (e1R:real) (e1F:real) (e2:real exp) (e2R:real) (e2F:real) (e3:real exp) (e3R:real) (e3F:real) (err1:real) (err2:real) (err3:real)
+     (vR:real) (vF:real) (E1 E2 :env) (m m1 m2 m3:mType) (defVars: num -> mType option).
+       eval_exp E1 (toRMap defVars) (toREval e1) e1R M0 /\
+       eval_exp E2 defVars e1 e1F m1 /\
+       eval_exp E1 (toRMap defVars) (toREval e2) e2R M0 /\
+       eval_exp E2 defVars e2 e2F m2 /\
+       eval_exp E1 (toRMap defVars) (toREval e3) e3R M0 /\
+       eval_exp E2 defVars e3 e3F m2 /\
+       eval_exp E1 (toRMap defVars) (toREval (Fma e1 e2 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 1) (Var 2) (Var 3)) vF m /\
+       abs (e1R - e1F) <= err1 /\
+       abs (e2R - e2F) <= err2 /\
+       abs (e3R - e3F) <= err3 ==>
+       abs (vR - vF) <= abs ((e1R - e1F) + (e2R * e3R - e2F * e3F)) + abs (e1F + e2F * e3F) * (mTypeToQ m)``,
+  rpt strip_tac
+  \\ qpat_x_assum `eval_exp E1 _ (toREval (Fma e1 e2 e3)) _ _` (fn thm => assume_tac (ONCE_REWRITE_RULE [toREval_def] thm))
+  \\ fs []
+  \\ inversion `eval_exp E1 _ (Fma _ _ _) _ _` eval_exp_cases
+  \\ rename1 `vR = perturb (evalFma v1R v2R v3R) deltaR`
+  \\ inversion `eval_exp _ _ (Fma (Var 1) (Var 2) (Var 3)) _ _` eval_exp_cases
+  \\ rename1 `vF = perturb (evalFma v1F v2F v3F) deltaF`
+  \\ `(m1' = M0) /\ (m2' = M0) /\ (m3' = M0)` by (rpt conj_tac \\ irule toRMap_eval_M0\\ asm_exists_tac \\ fs[]) \\ fs []
+  \\ rpt (qpat_x_assum `M0 = _` (fn thm => fs [GSYM thm]))
+  \\ `perturb (evalFma v1R v2R v3R) deltaR = evalFma v1R v2R v3R` by (match_mp_tac delta_M0_deterministic \\ fs[])
+  \\ `vR = evalFma v1R v2R v3R` by simp[]
+  \\ `v1R = e1R` by metis_tac[meps_0_deterministic]
+  \\ `v2R = e2R` by metis_tac[meps_0_deterministic]
+  \\ `v3R = e3R` by metis_tac[meps_0_deterministic]
+  \\ fs[evalFma_def, evalBinop_def, perturb_def]
+  \\ rpt (inversion `eval_exp (updEnv 3 e3F (updEnv 2 e2F (updEnv 1 e1F emptyEnv))) _ _ _ _` eval_exp_cases)
+  \\ fs [updEnv_def] \\ rveq
+  \\ fs [updDefVars_def] \\ rveq
+  \\ rewrite_tac [real_sub]
+  \\ `e1R + e2R * e3R + -((e1F + e2F * e3F) * (1 + deltaF)) = (e1R + e2R * e3R + -(e1F + e2F * e3F)) + (- (e1F + e2F * e3F) * deltaF)` by REAL_ASM_ARITH_TAC
+  \\ simp[]
+  \\ qspecl_then [`abs (e1R + e2R * e3R + -(e1F + e2F * e3F)) + abs (-(e1F + e2F * e3F) * deltaF)`] match_mp_tac real_le_trans2
+  \\ conj_tac
+  >- (REAL_ASM_ARITH_TAC)
+  >- (match_mp_tac REAL_LE_ADD2
+      \\ conj_tac \\ TRY (REAL_ASM_ARITH_TAC)
+      \\ once_rewrite_tac[REAL_ABS_MUL]
+      \\ match_mp_tac REAL_LE_MUL2 \\ fs[REAL_ABS_POS]
+      \\ once_rewrite_tac[GSYM REAL_NEG_LMUL, REAL_ABS_MUL]
+      \\ once_rewrite_tac[ABS_NEG] \\ fs[]));
+
 val round_abs_err_bounded = store_thm ("round_abs_err_bounded",
   ``!(e:real exp) (nR:real) (nF1:real) (nF:real) (E1:env) (E2:env) (err:real) (machineEpsilon:mType) (m:mType) (defVars: num -> mType option).
        eval_exp E1 (toRMap defVars) (toREval e) nR M0 /\