diff --git a/hol4/ErrorBoundsScript.sml b/hol4/ErrorBoundsScript.sml
index 5817a2582118f38391b1a1bfa8273e9fe88fc88e..3f87939a93e880c62f6d90459bf85e0ddfe9cf47 100644
--- a/hol4/ErrorBoundsScript.sml
+++ b/hol4/ErrorBoundsScript.sml
@@ -17,21 +17,32 @@ val triangle_tac =
   irule triangle_trans \\ rpt conj_tac \\ TRY (fs[REAL_ABS_TRIANGLE] \\ NO_TAC);
 
 Theorem const_abs_err_bounded:
-  !(n:real) (nR:real) (nF:real) (E1 E2:env) (m:mType) Gamma.
-    eval_expr E1 (toRTMap Gamma) (Const REAL n) nR REAL /\
-    eval_expr E2 Gamma (Const m n) nF m ==>
+  ∀ (n:real) (nR:real) (nF:real) (E1 E2:env) (m:mType) Gamma.
+    eval_expr E1 (toRTMap Gamma) (Const REAL n) nR REAL ∧
+    eval_expr E2 Gamma (Const m n) nF m ⇒
     abs (nR - nF) <= computeError n m
 Proof
   rpt strip_tac
-  \\ fs[eval_expr_cases]
-  \\ Cases_on `m` \\ fs[perturb_def, Rabs_err_simpl, REAL_ABS_MUL, computeError_def]
-  \\ rveq
+  \\ fs[eval_expr_cases, mTypeToR_def, perturb_def] \\ rveq
+  \\ Cases_on `m` \\ fs[Excl "RMUL_LEQNORM", perturb_def, computeError_def]
+  \\ TRY COND_CASES_TAC \\ fs[Excl "RMUL_LEQNORM", Rabs_err_simpl, REAL_ABS_MUL, mTypeToR_def]
   \\ TRY (REAL_ASM_ARITH_TAC)
-  \\ metis_tac [REAL_ABS_POS, REAL_LE_LMUL_IMP, REAL_MUL_COMM]
+  \\ irule REAL_LE_RMUL_IMP \\ rewrite_tac[GSYM REAL_INV_1OVER] \\ fs[REAL_ABS_POS]
 QED
 
 val float_add_tac =
-  rewrite_tac [GSYM REAL_ADD_ASSOC]
+  rewrite_tac [GSYM REAL_ADD_ASSOC, computeError_def]
+  \\ COND_CASES_TAC \\ simp[]
+  >- (
+    ‘e1R + (e2R + - (e1F + (e2F + deltaF))) = (e1R - e1F) + ((e2R - e2F) - deltaF)’
+    by REAL_ASM_ARITH_TAC
+    \\ pop_assum (rewrite_tac o single)
+    \\ triangle_tac
+    \\ irule REAL_LE_ADD2 \\ conj_tac \\ fs[]
+    \\ ‘e2R - e2F - deltaF = e2R - e2F + - deltaF’ by REAL_ASM_ARITH_TAC
+    \\ pop_assum (rewrite_tac o single)
+    \\ triangle_tac
+    \\ irule REAL_LE_ADD2 \\ conj_tac \\ fs[])
   \\ `e1R + (e2R + -((e1F + e2F) * (1 + deltaF))) =
       (e1R + - e1F) + ((e2R + - e2F) + - (e1F + e2F) * deltaF)`
         by REAL_ASM_ARITH_TAC
@@ -105,8 +116,17 @@ Proof
 QED
 
 val float_sub_tac =
-  rewrite_tac [GSYM REAL_ADD_ASSOC]
-  \\ `e1R + (-e2R + -((1 + deltaF) * (e1F + -e2F))) =
+  rewrite_tac [GSYM REAL_ADD_ASSOC, computeError_def]
+  \\ COND_CASES_TAC \\ simp[]
+  >- (‘e1R + (-e2R + -(e1F + (-e2F + deltaF))) = (e1R - e1F) + (-(e2R - e2F) + - deltaF)’
+      by REAL_ASM_ARITH_TAC
+      \\ pop_assum (rewrite_tac o single)
+      \\ triangle_tac \\ irule REAL_LE_ADD2 \\ conj_tac
+      >- fs[GSYM real_sub]
+      \\ triangle_tac \\ irule REAL_LE_ADD2 \\ conj_tac
+      >- fs[GSYM real_sub]
+      \\ fs[ABS_NEG, mTypeToR_def, real_sub])
+  \\ `e1R + (-e2R + -((e1F + -e2F) * (1 + deltaF))) =
     (e1R + - e1F) + ((- e2R + e2F) + - (e1F + - e2F) * deltaF)`
       by REAL_ASM_ARITH_TAC
    \\ pop_assum (once_rewrite_tac o single)
@@ -117,7 +137,7 @@ val float_sub_tac =
    \\ irule REAL_LE_ADD2 \\ conj_tac
    >- REAL_ASM_ARITH_TAC
    \\ rewrite_tac [REAL_ABS_MUL, ABS_NEG]
-   \\ metis_tac [REAL_ABS_POS, REAL_LE_LMUL_IMP, REAL_MUL_COMM];
+   \\ irule REAL_LE_LMUL_IMP \\ fs[mTypeToR_def, real_sub];
 
 Theorem subtract_abs_err_bounded:
   !(e1:real expr) (e1R:real) (e1F:real) (e2:real expr) (e2R:real) (e2F:real)
@@ -178,17 +198,25 @@ Proof
 QED
 
 val float_mul_tac =
-  `e1R * e2R + -(e1F * e2F * (1 + deltaF)) =
+  COND_CASES_TAC \\ simp[]
+  >- (‘e1R * e2R + - (e1F * e2F + deltaF) = (e1R * e2R + - (e1F * e2F)) + - deltaF’
+      by REAL_ASM_ARITH_TAC
+      \\ pop_assum (rewrite_tac o single)
+      \\ triangle_tac \\ irule REAL_LE_ADD2 \\ fs[computeError_def])
+  \\ `e1R * e2R + -(e1F * e2F * (1 + deltaF)) =
     (e1R * e2R + - (e1F * e2F)) + - (e1F * e2F * deltaF)`
       by REAL_ASM_ARITH_TAC
   \\ pop_assum (once_rewrite_tac o single)
-  \\ triangle_tac
-  \\ rewrite_tac [ABS_NEG, REAL_ABS_MUL] \\ rename1 ‘abs errF ≤ mTypeToR _’
-  \\ simp[computeError_def]
+  \\ triangle_tac \\ irule REAL_LE_ADD2 \\ conj_tac >- fs[]
+  \\ rewrite_tac [ABS_NEG, computeError_def, REAL_ABS_MUL]
+  \\ fs[]
+  \\ ‘abs deltaF * abs e1F * abs e2F = abs e1F * abs e2F * abs deltaF’
+    by fs[]
+  \\ pop_assum (rewrite_tac o single)
   \\ metis_tac [REAL_LE_LMUL_IMP, REAL_ABS_MUL, REAL_ABS_POS];
 
 Theorem mult_abs_err_bounded:
-  ! (e1:real expr) (e1R:real) (e1F:real) (e2:real expr) (e2R:real) (e2F:real)
+  ∀ (e1:real expr) (e1R:real) (e1F:real) (e2:real expr) (e2R:real) (e2F:real)
     (err1:real) (err2:real) (vR:real) (vF:real) (E1 E2 :env) (m m1 m2:mType)
     (Gamma: real expr -> mType option).
       eval_expr E1 (toRTMap Gamma) (toREval e1) e1R REAL /\
@@ -242,15 +270,21 @@ Proof
 QED
 
 val float_div_tac =
-  `e1R / e2R + -(e1F / e2F * (1 + deltaF)) =
-    (e1R / e2R + - (e1F / e2F)) + - (e1F / e2F * deltaF)`
+  COND_CASES_TAC \\ simp[]
+  >- (‘e1R / e2R + - (e1F / e2F + deltaF) = (e1R / e2R + - (e1F / e2F)) + - deltaF’
       by REAL_ASM_ARITH_TAC
-  \\ once_rewrite_tac [GSYM real_div]
-  \\ pop_assum (fn thm => once_rewrite_tac [thm])
-  \\ triangle_tac
-  \\ once_rewrite_tac [ABS_NEG]
-  \\ simp[computeError_def, REAL_ABS_MUL]
-  \\ metis_tac [REAL_LE_LMUL_IMP, REAL_ABS_POS, REAL_LE_MUL_SWAP];
+      \\ pop_assum (rewrite_tac o single)
+      \\ triangle_tac \\ irule REAL_LE_ADD2 \\ fs[computeError_def])
+  \\ `e1R / e2R + -(e1F * inv e2F * (1 + deltaF)) =
+    (e1R / e2R + - (e1F * inv e2F)) + - (e1F * inv e2F * deltaF)`
+      by REAL_ASM_ARITH_TAC
+  \\ pop_assum (once_rewrite_tac o single)
+  \\ triangle_tac \\ irule REAL_LE_ADD2 \\ conj_tac >- fs[real_div]
+  \\ fs[real_div, ABS_NEG, computeError_def, REAL_ABS_MUL]
+  \\ ‘abs deltaF * abs e1F = abs e1F * abs deltaF’
+    by fs[]
+  \\ pop_assum (rewrite_tac o single)
+  \\ metis_tac [REAL_LE_LMUL_IMP, REAL_ABS_MUL, REAL_ABS_POS];
 
 Theorem div_abs_err_bounded:
   !(e1:real expr) (e1R:real) (e1F:real) (e2:real expr) (e2R:real) (e2F:real)
@@ -300,7 +334,12 @@ Proof
 QED
 
 val float_fma_tac =
-  `e1R * e2R + e3R + -((1 + deltaF) * (e1F * e2F + e3F)) =
+  COND_CASES_TAC \\ simp[]
+>- (‘e1R * e2R + e3R + - (e1F * e2F + e3F + deltaF) =
+      (e1R * e2R + - (e1F * e2F)) + (e3R + - e3F) + - deltaF’ by REAL_ASM_ARITH_TAC
+    \\ pop_assum (rewrite_tac o single)
+    \\ triangle_tac \\ irule REAL_LE_ADD2 \\ fs[])
+  \\ `e1R * e2R + e3R + -((1 + deltaF) * (e1F * e2F + e3F)) =
     (e1R * e2R + e3R + -(e1F * e2F + e3F)) + (- (e1F * e2F + e3F) * deltaF)`
       by REAL_ASM_ARITH_TAC
   \\ simp[]
@@ -396,7 +435,9 @@ Proof
   \\ Cases_on `m1` \\ fs[perturb_def, computeError_def]
   \\ once_rewrite_tac [REAL_LDISTRIB]
   \\ simp[real_sub, REAL_NEG_ADD, REAL_ADD_ASSOC, ABS_NEG, ABS_MUL]
-  \\ metis_tac [REAL_LE_LMUL_IMP, ABS_POS, REAL_MUL_COMM, REAL_LE_MUL_SWAP]
+  \\ COND_CASES_TAC \\ fs[REAL_NEG_ADD, REAL_ADD_ASSOC]
+  \\ ‘delta * nF1 = nF1 * delta’ by fs[] \\ pop_assum (rewrite_tac o single)
+  \\ rewrite_tac [ABS_MUL] \\ irule REAL_LE_LMUL_IMP \\ fs[]
 QED
 
 Theorem nonzerop_EQ1_I'[simp]:
diff --git a/hol4/ErrorValidationScript.sml b/hol4/ErrorValidationScript.sml
index 5cef7fa12a1824a563f56433ac1ce0c5e635cd43..88b79df79180b367f440f3987f1e61cf40b6425b 100644
--- a/hol4/ErrorValidationScript.sml
+++ b/hol4/ErrorValidationScript.sml
@@ -29,7 +29,7 @@ val triangle_tac =
 fun real_rewrite t =
   rewrite_tac [REAL_ARITH (Parse.Term t)];
 
-val validErrorbound_def = Define `
+Definition validErrorbound_def:
   validErrorbound e (typeMap: typeMap) (A:analysisResult) (dVars:num_set)=
     case FloverMapTree_find e A, FloverMapTree_find e typeMap of
     | NONE, _ => F
@@ -114,13 +114,14 @@ val validErrorbound_def = Define `
                     mAbs = maxAbs errIve1 in
                 err1 + computeError mAbs m1 <= err
            else F
-        else F)`;
+       else F)
+End
 
 add_unevaluated_function ``validErrorbound``;
 add_unevaluated_function ``minAbsFun``;
 add_unevaluated_function ``noDivzero``;
 
-val validErrorboundCmd_def = Define `
+Definition validErrorboundCmd_def:
   validErrorboundCmd (f:real cmd) (typeMap: (real expr # mType) binTree)
                     (A:analysisResult) (dVars:num_set) =
     case f of
@@ -132,41 +133,46 @@ val validErrorboundCmd_def = Define `
              else F
            | _ , _ => F)
       | Ret e =>
-        validErrorbound e typeMap A dVars`;
+        validErrorbound e typeMap A dVars
+End
 
 add_unevaluated_function ``validErrorboundCmd``;
 
-val eval_Real_def = Define `
+Definition eval_Real_def:
   eval_Real E1 Gamma e nR =
-    eval_expr E1 (toRTMap (toRExpMap Gamma)) (toREval e) nR REAL`;
+    eval_expr E1 (toRTMap (toRExpMap Gamma)) (toREval e) nR REAL
+End
 
-val eval_Fin_def = Define `
+Definition eval_Fin_def:
   eval_Fin E2 Gamma e nF m =
-    eval_expr E2 (toRExpMap Gamma) e nF m`;
+    eval_expr E2 (toRExpMap Gamma) e nF m
+End
 
 val _ = export_rewrites ["eval_Real_def", "eval_Fin_def"]
 
-val err_always_positive = store_thm (
-  "err_always_positive",
-  ``! (e:real expr) (A:analysisResult) (iv:interval) (err:real) dVars tmap.
-      (validErrorbound e tmap A dVars) /\
-      (FloverMapTree_find e A = SOME (iv,err)) ==>
-      0 <= err``,
+Theorem err_always_positive:
+  ∀ (e:real expr) (A:analysisResult) (iv:interval) (err:real) dVars tmap.
+    (validErrorbound e tmap A dVars) /\
+    (FloverMapTree_find e A = SOME (iv,err)) ==>
+    0 <= err
+Proof
   rpt strip_tac \\ Cases_on `e`
-  \\ Flover_compute ``validErrorbound`` \\ rveq \\ fs[])
+  \\ qpat_x_assum `validErrorbound _ _ _ _` mp_tac
+  \\ simp[Once validErrorbound_def] \\ TOP_CASE_TAC \\ strip_tac
+QED
 
-val validErrorboundCorrectVariable_eval = store_thm (
-  "validErrorboundCorrectVariable_eval",
-  ``! E1 E2 A v e nR nlo nhi P fVars dVars Gamma.
-      eval_Real E1 Gamma (Var v) nR /\
-      approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
-      validTypes (Var v) Gamma /\
-      validRanges (Var v) A E1 (toRTMap (toRExpMap Gamma)) /\
-      validErrorbound (Var v) Gamma A dVars /\
-      (domain (usedVars ((Var v):real expr)) DIFF (domain dVars)) SUBSET (domain fVars) /\
-      FloverMapTree_find (Var v) A = SOME ((nlo,nhi),e) ==>
-      ? nF m.
-        eval_Fin E2 Gamma (Var v) nF m``,
+Theorem validErrorboundCorrectVariable_eval:
+  ∀ E1 E2 A v e nR nlo nhi P fVars dVars Gamma.
+    eval_Real E1 Gamma (Var v) nR /\
+    approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
+    validTypes (Var v) Gamma /\
+    validRanges (Var v) A E1 (toRTMap (toRExpMap Gamma)) /\
+    validErrorbound (Var v) Gamma A dVars /\
+    (domain (usedVars ((Var v):real expr)) DIFF (domain dVars)) SUBSET (domain fVars) /\
+    FloverMapTree_find (Var v) A = SOME ((nlo,nhi),e) ==>
+    ∃ nF m.
+      eval_Fin E2 Gamma (Var v) nF m
+Proof
   rpt strip_tac
   \\ IMP_RES_TAC validRanges_single \\ fs[]
   \\ IMP_RES_TAC meps_0_deterministic \\ rveq
@@ -177,30 +183,34 @@ val validErrorboundCorrectVariable_eval = store_thm (
       destruct approxEnv_gives_value
   \\ fs[domain_union, domain_lookup, usedVars_def]
   >- (Cases_on `lookup v dVars` \\ fs[domain_lookup])
-  \\ fs[eval_expr_cases, toRExpMap_def, toRTMap_def, option_case_eq]);
+  \\ fs[eval_expr_cases, toRExpMap_def, toRTMap_def, option_case_eq]
+QED
 
 add_unevaluated_function ``computeError``;
 add_unevaluated_function ``maxAbs``;
+add_unevaluated_function “normal”;
 
-val validErrorboundCorrectVariable = store_thm (
-  "validErrorboundCorrectVariable",
-  ``! (E1 E2:env) A fVars dVars  (v:num) (nR nF err nlo nhi:real) (P:precond) m
-      Gamma.
-      approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
-      eval_Real E1 Gamma (Var v) nR /\
-      eval_Fin E2 Gamma (Var v) nF m /\
-      validTypes (Var v) Gamma /\
-      validRanges (Var v) A E1 (toRTMap (toRExpMap Gamma)) /\
-      validErrorbound (Var v) Gamma A dVars /\
-      (domain (usedVars ((Var v):real expr)) DIFF (domain dVars)) SUBSET (domain fVars) /\
-      FloverMapTree_find (Var v) A = SOME ((nlo,nhi),err) ==>
-      abs (nR - nF) <= err``,
+Theorem validErrorboundCorrectVariable:
+  ∀ (E1 E2:env) A fVars dVars  (v:num) (nR nF err nlo nhi:real) (P:precond) m
+    Gamma.
+    approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
+    eval_Real E1 Gamma (Var v) nR /\
+    eval_Fin E2 Gamma (Var v) nF m /\
+    validTypes (Var v) Gamma /\
+    validRanges (Var v) A E1 (toRTMap (toRExpMap Gamma)) /\
+    validErrorbound (Var v) Gamma A dVars /\
+    (domain (usedVars ((Var v):real expr)) DIFF (domain dVars)) SUBSET (domain fVars) /\
+    FloverMapTree_find (Var v) A = SOME ((nlo,nhi),err) ==>
+    abs (nR - nF) <= err
+Proof
   rpt strip_tac
   \\ IMP_RES_TAC validRanges_single \\ fs[]
   \\ IMP_RES_TAC meps_0_deterministic \\ rveq
   \\ fs[toREval_def]
   \\ rpt (inversion `eval_expr _ _ _ _ _` eval_expr_cases)
-  \\ Flover_compute ``validErrorbound`` \\ rveq \\ fs[]
+  \\ qpat_x_assum `validErrorbound _ _ _ _` mp_tac
+  \\ simp[validErrorbound_def] \\ TOP_CASE_TAC \\ fs[]
+  \\ rpt strip_tac \\ fs[]
   >- (drule approxEnv_dVar_bounded
       \\ rpt (disch_then drule)
       \\ disch_then (qspecl_then [`m`, `(nlo,nhi)`, `e`] irule)
@@ -214,33 +224,32 @@ val validErrorboundCorrectVariable = store_thm (
   \\ qexists_tac `computeError nR m`
   \\ conj_tac \\ TRY (first_x_assum irule \\ fs[domain_lookup])
   \\ fs[toRExpMap_def] \\ rveq
-  \\ irule computeError_up
-  \\ irule contained_leq_maxAbs \\ fs[contained_def, IVlo_def, IVhi_def]);
+  \\ irule computeError_up \\ fs[]
+  \\ irule contained_leq_maxAbs \\ fs[contained_def, IVlo_def, IVhi_def]
+QED
 
 Theorem validErrorboundCorrectConstant_eval:
   !(E:env) (n nR:real) m Gamma.
-      validTypes (Const m n) Gamma ==>
-      ? nF m1.
-        eval_Fin E Gamma (Const m n) nF m1
+    validTypes (Const m n) Gamma ==>
+    ? nF m1.
+      eval_Fin E Gamma (Const m n) nF m1
 Proof
   rpt strip_tac \\ fs[validTypes_def]
-  \\ qexistsl_tac [`perturb n m (mTypeToR m)`,`m`] \\ irule Const_dist'
+  \\ qexistsl_tac [`perturb n m (mTypeToR m n)`,`m`] \\ irule Const_dist'
   \\ fs[]
-  \\ qexists_tac `mTypeToR m`
-  \\ fs[mTypeToR_def, realTheory.abs]
-  \\ Cases_on `m` \\ fs[mTypeToR_pos]
+  \\ qexists_tac `mTypeToR m n` \\ conj_tac \\ fs[realTheory.abs, mTypeToR_pos]
 QED
 
-val validErrorboundCorrectConstant = store_thm (
-  "validErrorboundCorrectConstant",
-  ``!(E1 E2:env) (A:analysisResult) (n nR nF e nlo nhi:real) dVars m Gamma.
-      eval_Real E1 Gamma (Const m n) nR /\
-      eval_Fin E2 Gamma (Const m n) nF m /\
-      validTypes (Const m n) Gamma /\
-      validErrorbound (Const m n) Gamma A dVars /\
-      (nlo <= nR /\ nR <= nhi) /\
-      FloverMapTree_find (Const m n) A = SOME ((nlo,nhi),e) ==>
-      (abs (nR - nF)) <= e``,
+Theorem validErrorboundCorrectConstant:
+  ∀(E1 E2:env) (A:analysisResult) (n nR nF e nlo nhi:real) dVars m Gamma.
+    eval_Real E1 Gamma (Const m n) nR /\
+    eval_Fin E2 Gamma (Const m n) nF m /\
+    validTypes (Const m n) Gamma /\
+    validErrorbound (Const m n) Gamma A dVars /\
+    (nlo <= nR /\ nR <= nhi) /\
+    FloverMapTree_find (Const m n) A = SOME ((nlo,nhi),e) ==>
+    (abs (nR - nF)) <= e
+Proof
   rpt strip_tac \\ fs[toREval_def]
   \\ Flover_compute ``validErrorbound``
   \\ rveq \\ fs[validTypes_def] \\ rveq
@@ -257,40 +266,41 @@ val validErrorboundCorrectConstant = store_thm (
       \\ conj_tac \\ irule Const_dist' \\ fs[]
       \\ find_exists_tac \\ fs[])
   \\ irule computeError_up
-  \\ fs[maxAbs_def] \\  irule maxAbs \\ fs[]);
+  \\ fs[maxAbs_def] \\  irule maxAbs \\ fs[]
+QED
 
-val validErrorboundCorrectAddition = store_thm (
-  "validErrorboundCorrectAddition",
-  ``!(E1 E2:env) (A:analysisResult) (e1:real expr) (e2:real expr)
+Theorem validErrorboundCorrectAddition:
+  ∀ (E1 E2:env) (A:analysisResult) (e1:real expr) (e2:real expr)
      (nR nR1 nR2 nF nF1 nF2:real) (e err1 err2:real)
      (alo ahi e1lo e1hi e2lo e2hi :real) dVars m m1 m2 Gamma.
-       eval_Real E1 Gamma e1 nR1 /\
-       eval_Real E1 Gamma e2 nR2 /\
-       eval_Real E1 Gamma (Binop Plus e1 e2) nR /\
-       eval_Fin E2 Gamma e1 nF1 m1 /\
-       eval_Fin E2 Gamma e2 nF2 m2 /\
-       eval_expr (updEnv 2 nF2 (updEnv 1 nF1 emptyEnv))
-          (updDefVars (Binop Plus (Var 1) (Var 2)) m
-            (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 (toRExpMap Gamma))))
-                (Binop Plus (Var 1) (Var 2)) nF m /\
-       validErrorbound (Binop Plus e1 e2) Gamma A dVars /\
-       (e1lo <= nR1 /\ nR1 <= e1hi) /\
-       (e2lo <= nR2 /\ nR2 <= e2hi) /\
-       FloverMapTree_find e1 A = SOME ((e1lo, e1hi), err1) /\
-       FloverMapTree_find e2 A = SOME ((e2lo, e2hi), err2) /\
-       FloverMapTree_find (Binop Plus e1 e2) A = SOME ((alo, ahi), e) /\
-       FloverMapTree_find (Binop Plus e1 e2) Gamma = SOME m /\
-       abs (nR1 - nF1) <= err1 /\
-       abs (nR2 - nF2) <= err2 ==>
-       abs (nR - nF) <= e``,
+    eval_Real E1 Gamma e1 nR1 /\
+    eval_Real E1 Gamma e2 nR2 /\
+    eval_Real E1 Gamma (Binop Plus e1 e2) nR /\
+    eval_Fin E2 Gamma e1 nF1 m1 /\
+    eval_Fin E2 Gamma e2 nF2 m2 /\
+    eval_expr (updEnv 2 nF2 (updEnv 1 nF1 emptyEnv))
+              (updDefVars (Binop Plus (Var 1) (Var 2)) m
+               (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 (toRExpMap Gamma))))
+              (Binop Plus (Var 1) (Var 2)) nF m /\
+    validErrorbound (Binop Plus e1 e2) Gamma A dVars /\
+    (e1lo <= nR1 /\ nR1 <= e1hi) /\
+    (e2lo <= nR2 /\ nR2 <= e2hi) /\
+    FloverMapTree_find e1 A = SOME ((e1lo, e1hi), err1) /\
+    FloverMapTree_find e2 A = SOME ((e2lo, e2hi), err2) /\
+    FloverMapTree_find (Binop Plus e1 e2) A = SOME ((alo, ahi), e) /\
+    FloverMapTree_find (Binop Plus e1 e2) Gamma = SOME m /\
+    abs (nR1 - nF1) <= err1 /\
+    abs (nR2 - nF2) <= err2 ==>
+    abs (nR - nF) <= e
+Proof
   rpt strip_tac \\ simp[Once toREval_def]
   \\ Flover_compute ``validErrorbound`` \\ rveq
   \\ fs[toRExpMap_def]
   \\ irule REAL_LE_TRANS
   \\ qexists_tac `err1 + err2 + (computeError (nF1 + nF2) m)`
   \\ conj_tac
-     >- (irule add_abs_err_bounded
-         \\ rpt (find_exists_tac \\ fs[]))
+  >- (irule add_abs_err_bounded
+      \\ rpt (find_exists_tac \\ fs[]))
   \\ irule REAL_LE_TRANS \\ find_exists_tac
   \\ conj_tac \\ fs[]
   \\ qmatch_abbrev_tac `_ <= computeError (maxAbs iv) _`
@@ -306,7 +316,8 @@ val validErrorboundCorrectAddition = store_thm (
            \\ qexists_tac `nR2` \\ conj_tac
            \\ simp[contained_def])
   \\ assume_tac (REWRITE_RULE [validIntervalAdd_def, contained_def] interval_addition_valid)
-  \\ fs[contained_def, widenInterval_def]);
+  \\ fs[contained_def, widenInterval_def]
+QED
 
 val validErrorboundCorrectSubtraction = store_thm ("validErrorboundCorrectSubtraction",
   ``!(E1 E2:env) (A:analysisResult) (e1:real expr) (e2:real expr)
@@ -2219,21 +2230,22 @@ val validErrorboundCorrectRounding = store_thm ("validErrorboundCorrectRounding"
    Proof is by induction on the expression e.
    Each case requires the application of one of the soundness lemmata proven before
  **)
-val validErrorbound_sound = store_thm ("validErrorbound_sound",
-  ``! (e:real expr) (E1 E2:env) (A:analysisResult) (nR err:real)
-     (elo ehi:real) (fVars:num_set) dVars Gamma.
-       validTypes e Gamma /\
-       approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
-       ((domain (usedVars e)) DIFF (domain dVars)) SUBSET (domain fVars) /\
-       eval_Real E1 Gamma e nR /\
-       validErrorbound e Gamma A dVars /\
-       validRanges e A E1 (toRTMap (toRExpMap Gamma)) /\
-       FloverMapTree_find e A = SOME ((elo,ehi),err) ==>
-       (?nF m.
-         eval_Fin E2 Gamma e nF m) /\
-       (! nF m.
-         eval_Fin E2 Gamma e nF m ==>
-         abs (nR - nF) <= err)``,
+Theorem validErrorbound_sound:
+  ∀ (e:real expr) (E1 E2:env) (A:analysisResult) (nR err:real)
+    (elo ehi:real) (fVars:num_set) dVars Gamma.
+    validTypes e Gamma /\
+    approxEnv E1 (toRExpMap Gamma) A fVars dVars E2 /\
+    ((domain (usedVars e)) DIFF (domain dVars)) SUBSET (domain fVars) /\
+    eval_Real E1 Gamma e nR /\
+    validErrorbound e Gamma A dVars /\
+    validRanges e A E1 (toRTMap (toRExpMap Gamma)) /\
+    FloverMapTree_find e A = SOME ((elo,ehi),err) ==>
+    (?nF m.
+      eval_Fin E2 Gamma e nF m) /\
+    (! nF m.
+      eval_Fin E2 Gamma e nF m ==>
+      abs (nR - nF) <= err)
+Proof
   Induct_on `e`
   \\ rpt gen_tac
   \\ rpt (disch_then assume_tac)
@@ -2355,7 +2367,7 @@ val validErrorbound_sound = store_thm ("validErrorbound_sound",
       \\ fs[]
       \\ qpat_assum `eval_expr E2 _ (Binop op e1 e2) _ _` assume_tac
       \\ inversion `eval_expr E2 _ (Binop op e1 e2) _ _` eval_expr_cases
-      \\ rename1 `abs delta2 <= mTypeToR mF`
+      \\ rename1 `abs delta2 <= mTypeToR mF _`
       \\ rveq
       \\ `m2' = m2`
         by (irule validTypes_exec
@@ -2470,7 +2482,7 @@ val validErrorbound_sound = store_thm ("validErrorbound_sound",
         by (irule validTypes_exec
             \\ qexistsl_tac [`E2`, `Gamma`, `e3`, `v3`] \\ fs[])
       \\ rveq
-      \\ rename1 `abs delta2 <= mTypeToR mF`
+      \\ rename1 `abs delta2 <= mTypeToR mF _`
       \\ `mF = x'`
         by (irule validTypes_exec
             \\ qexistsl_tac [`E2`, `Gamma`, `Fma e1 e2 e3`,
@@ -2543,7 +2555,8 @@ val validErrorbound_sound = store_thm ("validErrorbound_sound",
       \\ qexistsl_tac [`delta`, `v1`]
       \\ fs[]
       \\ irule Var_load
-      \\ fs[updDefVars_def, updEnv_def]));
+      \\ fs[updDefVars_def, updEnv_def])
+QED
 
 Theorem validErrorboundCmd_gives_eval:
   !(f:real cmd) (A:analysisResult) (E1 E2:env)
diff --git a/hol4/ExpressionSemanticsScript.sml b/hol4/ExpressionSemanticsScript.sml
index 88e20861317099fda428921e881dbb33cad00a7d..738f39ecc323c148f42caecc384063b81ab25281 100644
--- a/hol4/ExpressionSemanticsScript.sml
+++ b/hol4/ExpressionSemanticsScript.sml
@@ -12,13 +12,17 @@ val _ = temp_overload_on("abs",``real$abs``);
 (**
   Define a perturbation function to ease writing of basic definitions
 **)
-val perturb_def = Define `
+Definition perturb_def:
   (* The Real-type has no error *)
   perturb (rVal:real) (REAL) (delta:real) = rVal /\
   (* Fixed-point numbers have an absolute error *)
   perturb rVal (F w f) delta = rVal + delta /\
-  (* Floating-point numbers have a relative error *)
-  perturb rVal _ delta = rVal * (1 + delta)`;
+  (* Floating-point numbers have a relative error in the normal range, and an
+     absolute error in the subnormal range *)
+  perturb rVal m delta = if (denormal rVal m)
+                         then rVal + delta
+                         else rVal * (1 + delta)
+End
 
 (**
 Define exprression evaluation relation parametric by an "error" epsilon.
@@ -26,13 +30,13 @@ The result value exprresses float computations according to the IEEE standard,
 using a perturbation of the real valued computation by (1 + delta), where
 |delta| <= machine epsilon.
 **)
-val (eval_expr_rules, eval_expr_ind, eval_expr_cases) = Hol_reln `
+Inductive eval_expr:
   (!E Gamma m x v.
      Gamma (Var x) = SOME m /\
      E x = SOME v ==>
      eval_expr E Gamma (Var x) v m) /\
   (!E Gamma m n delta.
-      abs delta <= (mTypeToR m) ==>
+      abs delta <= (mTypeToR m) n ==>
      eval_expr E Gamma (Const m n) (perturb n m delta) m) /\
   (!E Gamma m mN f1 v1.
      Gamma (Unop Neg f1) = SOME mN /\
@@ -42,20 +46,20 @@ val (eval_expr_rules, eval_expr_ind, eval_expr_cases) = Hol_reln `
   (!E Gamma m mN f1 v1 delta.
      Gamma (Unop Inv f1) = SOME mN /\
      isCompat m mN /\
-     abs delta <= (mTypeToR m) /\
+     abs delta <= (mTypeToR m) (evalUnop Inv v1) /\
      eval_expr E Gamma f1 v1 m /\
      (v1 <> 0) ==>
      eval_expr E Gamma (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) mN) /\
   (!E Gamma m m1 f1 v1 delta.
      Gamma (Downcast m f1) = SOME m /\
      isMorePrecise m1 m /\
-     abs delta <= (mTypeToR m) /\
+     abs delta <= (mTypeToR m) v1 /\
      eval_expr E Gamma f1 v1 m1 ==>
      eval_expr E Gamma (Downcast m f1) (perturb v1 m delta) m) /\
   (!E Gamma m1 m2 m op f1 f2 v1 v2 delta.
      Gamma (Binop op f1 f2) = SOME m /\
      isJoin m1 m2 m /\
-     abs delta <= (mTypeToR m) /\
+     abs delta <= (mTypeToR m) (evalBinop op v1 v2) /\
      eval_expr E Gamma f1 v1 m1 /\
      eval_expr E Gamma f2 v2 m2 /\
      ((op = Div) ==> (v2 <> 0)) ==>
@@ -63,18 +67,19 @@ val (eval_expr_rules, eval_expr_ind, eval_expr_cases) = Hol_reln `
   (!E Gamma m1 m2 m3 m f1 f2 f3 v1 v2 v3 delta.
      Gamma (Fma f1 f2 f3) = SOME m /\
      isJoin3 m1 m2 m3 m /\
-     abs delta <= (mTypeToR m) /\
+     abs delta <= (mTypeToR m) (evalFma v1 v2 v3) /\
      eval_expr E Gamma f1 v1 m1 /\
      eval_expr E Gamma f2 v2 m2 /\
      eval_expr E Gamma f3 v3 m3 ==>
-     eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m)`;
+     eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m)
+End
 
-val eval_expr_cases_old = save_thm ("eval_expr_cases_old", eval_expr_cases);
+Theorem eval_expr_cases_old = eval_expr_cases
 
 (**
   Generate a better case lemma
 **)
-val eval_expr_cases =
+Theorem eval_expr_cases =
   map (GEN_ALL o SIMP_CONV (srw_ss()) [Once eval_expr_cases])
     [``eval_expr E Gamma (Var v) res m``,
      ``eval_expr E Gamma (Const m1 n) res m2``,
@@ -87,140 +92,153 @@ val eval_expr_cases =
 val [Var_load, Const_dist, Unop_neg, Unop_inv, Downcast_dist, Binop_dist, Fma_dist] =
   CONJ_LIST 7 eval_expr_rules;
 
-save_thm ("Var_load", Var_load);
-save_thm ("Const_dist", Const_dist);
-save_thm ("Unop_neg", Unop_neg);
-save_thm ("Unop_inv", Unop_inv);
-save_thm ("Binop_dist", Binop_dist);
-save_thm ("Fma_dist", Fma_dist);
-save_thm ("Downcast_dist", Downcast_dist);
+Theorem Var_load = Var_load
+Theorem Const_dist = Const_dist
+Theorem Unop_neg = Unop_neg
+Theorem Unop_inv = Unop_inv
+Theorem Binop_dist = Binop_dist
+Theorem Fma_dist = Fma_dist
+Theorem Downcast_dist = Downcast_dist
 
 (**
   Show some simpler (more general) rule lemmata
 **)
+Theorem Const_dist':
+  ∀ m n delta v m' E Gamma.
+    abs delta <= mTypeToR m n /\
+    v = perturb n m delta /\
+    m' = m ==>
+    eval_expr E Gamma (Const m n) v m'
+Proof
+  fs [Const_dist]
+QED
 
-val Const_dist' = store_thm (
-  "Const_dist'",
-  ``!m n delta v m' E Gamma.
-      (abs delta) <= (mTypeToR m) /\
-      v = perturb n m delta /\
-      m' = m ==>
-      eval_expr E Gamma (Const m n) v m'``,
-  fs [Const_dist]);
-
-val Unop_neg' = store_thm (
-  "Unop_neg'",
-  ``!m f1 v1 v m' mN E Gamma.
-      eval_expr E Gamma f1 v1 m /\
-      v = evalUnop Neg v1 /\
-      Gamma (Unop Neg f1) = SOME mN /\
-      isCompat m mN /\
-      m' = mN ==>
-      eval_expr E Gamma (Unop Neg f1) v m'``,
+Theorem Unop_neg':
+  ∀ m f1 v1 v m' mN E Gamma.
+    eval_expr E Gamma f1 v1 m /\
+    v = evalUnop Neg v1 /\
+    Gamma (Unop Neg f1) = SOME mN /\
+    isCompat m mN /\
+    m' = mN ==>
+    eval_expr E Gamma (Unop Neg f1) v m'
+Proof
   rpt strip_tac \\ rveq \\ irule Unop_neg \\ fs[]
-  \\ asm_exists_tac \\ fs[]);
+  \\ asm_exists_tac \\ fs[]
+QED
 
-val Unop_inv' = store_thm (
-  "Unop_inv'",
-  ``!m f1 v1 delta v m' mN E Gamma fBits.
-      (abs delta) <= (mTypeToR m) /\
-      eval_expr E Gamma f1 v1 m /\
-      (v1 <> 0) /\
-      v = perturb (evalUnop Inv v1) m delta /\
-      Gamma (Unop Inv f1) = SOME mN /\
-      isCompat m mN /\
-      m' = mN ==>
-      eval_expr E Gamma (Unop Inv f1) v m'``,
-  rpt strip_tac \\ rveq \\ irule Unop_inv \\ fs[]);
+Theorem Unop_inv':
+  ∀ m f1 v1 delta v m' mN E Gamma fBits.
+    (abs delta) <= mTypeToR m (evalUnop Inv v1) /\
+    eval_expr E Gamma f1 v1 m /\
+    (v1 <> 0) /\
+    v = perturb (evalUnop Inv v1) m delta /\
+    Gamma (Unop Inv f1) = SOME mN /\
+    isCompat m mN /\
+    m' = mN ==>
+    eval_expr E Gamma (Unop Inv f1) v m'
+Proof
+  rpt strip_tac \\ rveq \\ irule Unop_inv \\ fs[]
+QED
 
-val Downcast_dist' = store_thm ("Downcast_dist'",
-  ``!m m1 f1 v1 delta v m' E Gamma fBits.
-      isMorePrecise m1 m /\
-      (abs delta) <= (mTypeToR m) /\
-      eval_expr E Gamma f1 v1 m1 /\
-      v = perturb v1 m delta /\
-      Gamma (Downcast m f1) = SOME m /\
-      m' = m ==>
-      eval_expr E Gamma (Downcast m f1) v m'``,
+Theorem Downcast_dist':
+  ∀ m m1 f1 v1 delta v m' E Gamma fBits.
+    isMorePrecise m1 m /\
+    (abs delta) <= mTypeToR m v1 /\
+    eval_expr E Gamma f1 v1 m1 /\
+    v = perturb v1 m delta /\
+    Gamma (Downcast m f1) = SOME m /\
+    m' = m ==>
+    eval_expr E Gamma (Downcast m f1) v m'
+Proof
   rpt strip_tac
   \\ rw []
   \\ irule Downcast_dist
-  \\ fs[] \\ qexists_tac `m1` \\ fs[]);
+  \\ fs[] \\ qexists_tac `m1` \\ fs[]
+QED
 
-val Binop_dist' = store_thm ("Binop_dist'",
-  ``!m1 m2 op f1 f2 v1 v2 delta v m m' E Gamma fBit.
-      (abs delta) <= (mTypeToR m') /\
-      eval_expr E Gamma f1 v1 m1 /\
-      eval_expr E Gamma f2 v2 m2 /\
-      ((op = Div) ==> (v2 <> 0)) /\
-      v = perturb (evalBinop op v1 v2) m' delta /\
-      Gamma (Binop op f1 f2) = SOME m /\
-      isJoin m1 m2 m /\
-      m = m' ==>
-      eval_expr E Gamma (Binop op f1 f2) v m'``,
+Theorem Binop_dist':
+  ∀ m1 m2 op f1 f2 v1 v2 delta v m m' E Gamma fBit.
+    abs delta <= mTypeToR m' (evalBinop op v1 v2) /\
+    eval_expr E Gamma f1 v1 m1 /\
+    eval_expr E Gamma f2 v2 m2 /\
+    ((op = Div) ==> (v2 <> 0)) /\
+    v = perturb (evalBinop op v1 v2) m' delta /\
+    Gamma (Binop op f1 f2) = SOME m /\
+    isJoin m1 m2 m /\
+    m = m' ==>
+    eval_expr E Gamma (Binop op f1 f2) v m'
+Proof
   rpt strip_tac \\ rveq
   \\ irule Binop_dist \\ fs[]
-  \\ asm_exists_tac \\ fs[]);
+  \\ asm_exists_tac \\ fs[]
+QED
 
-val Fma_dist' = store_thm ("Fma_dist'",
-  ``!m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m m' E Gamma fBit.
-      (abs delta) <= (mTypeToR m') /\
-      eval_expr E Gamma f1 v1 m1 /\
-      eval_expr E Gamma f2 v2 m2 /\
-      eval_expr E Gamma f3 v3 m3 /\
-      v = perturb (evalFma v1 v2 v3) m' delta /\
-      Gamma (Fma f1 f2 f3) = SOME m /\
-      isJoin3 m1 m2 m3 m /\
-      m' = m ==>
-      eval_expr E Gamma (Fma f1 f2 f3) v m'``,
+Theorem Fma_dist':
+  ∀ m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m m' E Gamma fBit.
+    abs delta <= mTypeToR m' (evalFma v1 v2 v3) /\
+    eval_expr E Gamma f1 v1 m1 /\
+    eval_expr E Gamma f2 v2 m2 /\
+    eval_expr E Gamma f3 v3 m3 /\
+    v = perturb (evalFma v1 v2 v3) m' delta /\
+    Gamma (Fma f1 f2 f3) = SOME m /\
+    isJoin3 m1 m2 m3 m /\
+    m' = m ==>
+    eval_expr E Gamma (Fma f1 f2 f3) v m'
+Proof
   rpt strip_tac \\ rveq
   \\ irule Fma_dist \\ fs[]
-  \\ asm_exists_tac \\ fs[]);
+  \\ asm_exists_tac \\ fs[]
+QED
 
-val Gamma_det = store_thm (
-  "Gamma_det",
-  ``!e1 e2 E1 E2 Gamma v1 v2 m1 m2.
-      eval_expr E1 Gamma e1 v1 m1 /\
-      e1 = e2 /\
-      eval_expr E2 Gamma e2 v2 m2 ==>
-      m1 = m2``,
+Theorem Gamma_det:
+  ∀ e1 e2 E1 E2 Gamma v1 v2 m1 m2.
+    eval_expr E1 Gamma e1 v1 m1 /\
+    e1 = e2 /\
+    eval_expr E2 Gamma e2 v2 m2 ==>
+    m1 = m2
+Proof
   Induct_on `e1` \\ rpt strip_tac \\ fs[eval_expr_cases]
   \\ rveq \\ fs[eval_expr_cases]
-  \\ fs[]);
+  \\ fs[]
+QED
 
-val toRTMap_eval_REAL = store_thm (
-  "toRTMap_eval_REAL",
-  ``!e v E Gamma m.
-      eval_expr E (toRTMap Gamma) (toREval e) v m ==> m = REAL``,
+Theorem toRTMap_eval_REAL:
+  ∀ e v E Gamma m.
+    eval_expr E (toRTMap Gamma) (toREval e) v m ==> m = REAL
+Proof
   Induct_on `e` \\ rpt strip_tac \\ fs[toREval_def]
   \\ inversion `eval_expr _ _ _ _ _` eval_expr_cases
   \\ fs[toRTMap_def, option_case_eq]
-  \\ res_tac \\ fs[]);
+  \\ res_tac \\ fs[]
+QED
 
 (**
   If |delta| <= 0 then perturb v delta is exactly v.
 **)
-val delta_0_deterministic = store_thm(
-  "delta_0_deterministic",
-  ``!(v:real) (m:mType) (delta:real).
-      abs delta <= 0 ==> perturb v m delta = v``,
+Theorem delta_0_deterministic:
+  ∀ (v:real) (m:mType) (delta:real).
+    abs delta <= 0 ==> perturb v m delta = v
+Proof
   Cases_on `m` \\
-  fs [perturb_def,ABS_BOUNDS,REAL_LE_ANTISYM]);
+  fs [perturb_def,ABS_BOUNDS,REAL_LE_ANTISYM]
+QED
 
-val delta_REAL_deterministic = store_thm(
-  "delta_REAL_deterministic",
-  ``!(v:real) (m:mType) (delta:real).
-      abs delta <= mTypeToR REAL ==> perturb v m delta = v``,
-  Cases_on `m` \\ fs[mTypeToR_def, delta_0_deterministic]);
+Theorem delta_REAL_deterministic:
+  ∀ (v:real) (m:mType) (delta:real).
+    abs delta <= mTypeToR REAL v ==> perturb v m delta = v
+Proof
+  Cases_on `m` \\ fs[mTypeToR_def, delta_0_deterministic]
+QED
 
 (**
 Evaluation with 0 as machine epsilon is deterministic
 **)
-val meps_0_deterministic = store_thm("meps_0_deterministic",
-  ``!(f: real expr) v1:real v2:real E defVars fBits.
-      eval_expr E (toRTMap defVars) (toREval f) v1 REAL /\
-      eval_expr E (toRTMap defVars) (toREval f) v2 REAL ==>
-      v1 = v2``,
+Theorem meps_0_deterministic:
+  ∀ (f: real expr) v1:real v2:real E defVars fBits.
+    eval_expr E (toRTMap defVars) (toREval f) v1 REAL /\
+    eval_expr E (toRTMap defVars) (toREval f) v2 REAL ==>
+    v1 = v2
+Proof
   Induct_on `f` \\ fs[toREval_def] \\ rpt strip_tac \\ fs[eval_expr_cases]
   \\ rveq \\ fs[delta_REAL_deterministic]
   >- (IMP_RES_TAC toRTMap_eval_REAL \\ rveq \\ res_tac \\ fs[])
@@ -231,7 +249,8 @@ val meps_0_deterministic = store_thm("meps_0_deterministic",
   >- (IMP_RES_TAC toRTMap_eval_REAL \\ rveq \\ res_tac
       \\ fs[delta_REAL_deterministic])
   \\ IMP_RES_TAC toRTMap_eval_REAL \\ rveq \\ res_tac
-  \\ fs[delta_REAL_deterministic]);
+  \\ fs[delta_REAL_deterministic]
+QED
 
 (**
 Helping lemmas. Needed in soundness proof.
@@ -239,44 +258,50 @@ For each evaluation of using an arbitrary epsilon, we can replace it by
 evaluating the subExpressions and then binding the result values to different
 variables in the Environment.
 **)
-val binary_unfolding = store_thm("binary_unfolding",
-  ``!(b:binop) (f1:(real)expr) (f2:(real)expr) E Gamma v1 v2 m1 m2 delta m.
-      (b = Div ==> (v2 <> 0)) /\
-      (abs delta) <= (mTypeToR m) /\
-      eval_expr E Gamma f1 v1 m1 /\
-      eval_expr E Gamma f2 v2 m2 /\
-      eval_expr E Gamma (Binop b f1 f2) (perturb (evalBinop b v1 v2) m delta) m ==>
-      eval_expr (updEnv 2 v2 (updEnv 1 v1 emptyEnv))
-        (updDefVars (Binop b (Var 1) (Var 2)) m (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 Gamma)))
-        (Binop b (Var 1) (Var 2)) (perturb (evalBinop b v1 v2) m delta) m``,
+Theorem binary_unfolding:
+  ∀ (b:binop) (f1:(real)expr) (f2:(real)expr) E Gamma v1 v2 m1 m2 delta m.
+    (b = Div ==> (v2 <> 0)) /\
+    (abs delta) <= (mTypeToR m) (evalBinop b v1 v2) /\
+    eval_expr E Gamma f1 v1 m1 /\
+    eval_expr E Gamma f2 v2 m2 /\
+    eval_expr E Gamma (Binop b f1 f2) (perturb (evalBinop b v1 v2) m delta) m ==>
+    eval_expr (updEnv 2 v2 (updEnv 1 v1 emptyEnv))
+              (updDefVars (Binop b (Var 1) (Var 2)) m
+               (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 Gamma)))
+              (Binop b (Var 1) (Var 2)) (perturb (evalBinop b v1 v2) m delta) m
+Proof
   fs [updEnv_def,updDefVars_def,join_fl_def,eval_expr_cases,APPLY_UPDATE_THM,PULL_EXISTS]
   \\ rpt strip_tac
   \\ qexists_tac `delta` \\ fs[]
-  \\ IMP_RES_TAC Gamma_det \\ fs[]);
+  \\ IMP_RES_TAC Gamma_det \\ fs[]
+QED
 
-val fma_unfolding = store_thm("fma_unfolding",
-  ``!(f1:(real)expr) (f2:(real)expr) (f3:(real)expr) E Gamma (v:real) v1 v2 v3
-     m1 m2 m3 delta fBit m.
-      (abs delta) <= (mTypeToR m) /\
-      eval_expr E Gamma f1 v1 m1 /\
-      eval_expr E Gamma f2 v2 m2 /\
-      eval_expr E Gamma f3 v3 m3 /\
-      eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m ==>
-      eval_expr (updEnv 3 v3 (updEnv 2 v2 (updEnv 1 v1 emptyEnv)))
-        (updDefVars (Fma (Var 1) (Var 2) (Var 3)) m
-          (updDefVars (Var 3) m3 (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 Gamma))))
-        (Fma (Var 1) (Var 2) (Var 3)) (perturb (evalFma v1 v2 v3) m delta) m``,
+Theorem fma_unfolding:
+  ∀(f1:(real)expr) (f2:(real)expr) (f3:(real)expr) E Gamma (v:real) v1 v2 v3
+     m1 m2 m3 delta m.
+    (abs delta) <= (mTypeToR m) (evalFma v1 v2 v3) /\
+    eval_expr E Gamma f1 v1 m1 /\
+    eval_expr E Gamma f2 v2 m2 /\
+    eval_expr E Gamma f3 v3 m3 /\
+    eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m ==>
+    eval_expr (updEnv 3 v3 (updEnv 2 v2 (updEnv 1 v1 emptyEnv)))
+              (updDefVars (Fma (Var 1) (Var 2) (Var 3)) m
+               (updDefVars (Var 3) m3
+                (updDefVars (Var 2) m2 (updDefVars (Var 1) m1 Gamma))))
+              (Fma (Var 1) (Var 2) (Var 3)) (perturb (evalFma v1 v2 v3) m delta) m
+Proof
   fs [updEnv_def,updDefVars_def,join_fl3_def,join_fl_def,eval_expr_cases,APPLY_UPDATE_THM,PULL_EXISTS]
   \\ rpt strip_tac
   \\ qexists_tac `delta` \\ fs[]
-  \\ IMP_RES_TAC Gamma_det \\ fs[]);
+  \\ IMP_RES_TAC Gamma_det \\ fs[]
+QED
 
-val eval_eq_env = store_thm (
-  "eval_eq_env",
-  ``!e E1 E2 Gamma v m.
-      (!x. E1 x = E2 x) /\
-      eval_expr E1 Gamma e v m ==>
-      eval_expr E2 Gamma e v m``,
+Theorem eval_eq_env:
+  ∀ e E1 E2 Gamma v m.
+    (!x. E1 x = E2 x) /\
+    eval_expr E1 Gamma e v m ==>
+    eval_expr E2 Gamma e v m
+Proof
   Induct \\ rpt strip_tac \\ fs[eval_expr_cases]
   >- (`E1 n = E2 n` by (first_x_assum irule)
       \\ fs[])
@@ -291,15 +316,17 @@ val eval_eq_env = store_thm (
   >- (rveq \\ qexistsl_tac [`m1`, `m2`, `m3`, `v1`, `v2`, `v3`, `delta`]
       \\ fs[] \\ prove_tac [])
   >- (rveq \\ qexistsl_tac [`m1'`, `v1`, `delta`] \\ fs[]
-      \\ first_x_assum drule \\ disch_then irule \\ fs[]));
+      \\ first_x_assum drule \\ disch_then irule \\ fs[])
+QED
 
-val swap_Gamma_eval_expr = store_thm (
-  "swap_Gamma_eval_expr",
-  ``!e E vR m Gamma1 Gamma2.
+Theorem swap_Gamma_eval_expr:
+  ∀ e E vR m Gamma1 Gamma2.
       (!e. Gamma1 e = Gamma2 e) /\
       eval_expr E Gamma1 e vR m ==>
-      eval_expr E Gamma2 e vR m``,
+      eval_expr E Gamma2 e vR m
+Proof
   Induct_on `e` \\ fs[eval_expr_cases] \\ rpt strip_tac
-  \\ metis_tac[]);
+  \\ metis_tac[]
+QED
 
 val _ = export_theory();
diff --git a/hol4/Infra/MachineTypeScript.sml b/hol4/Infra/MachineTypeScript.sml
index 03eb5d928d5392bda7c7373c4d89b927df1c5c3e..7945c5e75c126e52247e5e7d7164f306ae516d22 100644
--- a/hol4/Infra/MachineTypeScript.sml
+++ b/hol4/Infra/MachineTypeScript.sml
@@ -4,7 +4,7 @@
   @author: Raphael Monat
   @maintainer: Heiko Becker
  **)
-open realTheory realLib sptreeTheory;
+open realTheory realLib;
 open RealSimpsTheory;
 open preambleFloVer;
 
@@ -14,6 +14,7 @@ val _ = temp_overload_on("abs",``real$abs``);
 val _ = monadsyntax.enable_monadsyntax();
 val _ = List.app monadsyntax.enable_monad ["option"];
 
+
 Datatype:
   mType = REAL | M16 | M32 | M64 | F num num (* first num is word length, second is fractional bits *)
 End
@@ -23,16 +24,6 @@ Definition isFixedPoint_def :
   isFixedPoint _ = F
 End
 
-Definition mTypeToR_def:
-  mTypeToR (m:mType) : real =
-    case m of
-      | REAL => 0
-      | M16 => 1 / (2 pow 11)
-      | M32 => 1 / (2 pow 24)
-      | M64 => 1 / (2 pow 53)
-      | F w f => 1 / (2 pow f)
-End
-
 Definition maxExponent_def:
   (maxExponent (REAL) = 0n) /\
   (maxExponent (M16) = 15) /\
@@ -54,7 +45,6 @@ Goldberg - Handbook of Floating-Point Arithmetic: (p.183)
 (𝛃 - 𝛃^(1 - p)) * 𝛃^(e_max)
 which simplifies to 2^(e_max) for base 2
 **)
-
 Definition maxValue_def:
   maxValue (m:mType) =
     case m of
@@ -71,33 +61,70 @@ Definition minValue_pos_def:
     | _ =>  inv (&(2n ** (minExponentPos m)))
 End
 
+(** Goldberg - Handbook of Floating-Point Arithmetic: (p.183)
+  𝛃^(e_min -p + 1) = 𝛃^(e_min -1) for base 2
+**)
+Definition normal_def:
+  normal (v:real) (m:mType) =
+  (minValue_pos m <= abs v /\ abs v <= maxValue m)
+End
+
+Definition denormal_def:
+  denormal (v:real) (m:mType) =
+    case m of
+      | REAL => F
+      | _ => ((abs v) < (minValue_pos m) /\ v <> 0)
+End
+
+(** Return an upperbound on the error for value v committed by rounding
+  For real numbers this is 0,
+  For floating-point numbers, we have to first check whether the the value is
+    denormal or normal. For denormal numbers, we use the upper bound
+    1/2 * 2^(-e_min - p + 1) = 2^(-1) * 2^(-e_min - p + 1)) =
+    2 ^(-1 - e_min - p + 1) = 2^(-(e_min + p)) = inv (2 ^(e_min + p))
+    where p is the number of precision bits for the floating-point type
+    For normal numbers the roundoff error is at most
+      1/2 * 2^(1 - p) = 2^(-1) * 2^(1-p) = 2^(-1 + 1 -p) = 2^(-p) = inv (2^p)
+    where p is the number of precision bits.
+  For fixed-point numbers with f fractional bits, the truncation error is at
+    maximum inv (2 pow f)
+  see p. 24 of Handbook of Floating-Point Arithmetic for details *)
+Definition mTypeToR_def:
+  mTypeToR (m:mType) (v:real):real =
+    case m of
+      | REAL => 0
+      | M16 => if (denormal v M16)
+               then inv (2 pow (minExponentPos M16 + 11))
+               else inv (2 pow 11)
+      | M32 => if (denormal v M32)
+               then inv(2 pow (minExponentPos M32 + 24))
+               else inv (2 pow 24)
+      | M64 => if (denormal v M64)
+               then inv (2 pow (minExponentPos M64 + 53))
+               else inv (2 pow 53)
+      | F w f => inv (2 pow f)
+End
+
+Theorem mTypeToR_simp =
+  SIMP_RULE std_ss
+    (List.map EVAL [“minExponentPos M16”, “minExponentPos M32”, “minExponentPos M64”])
+    mTypeToR_def
 
+(** Compute the maximum error contributed by rounding a value v to type m,
+    does not account for propagation errors, these must be accounted elsewhere *)
 Definition computeError_def:
   computeError v m =
     case m of
     | REAL => 0
-    | F w f => mTypeToR m
-    | _ => abs v * mTypeToR m
+    | F w f => mTypeToR m v
+    | _ => if (denormal v m) then mTypeToR m v else abs v * mTypeToR m v
 End
 
 Theorem mTypeToR_pos:
-  ! e. 0 <= mTypeToR e
+  ∀ m v. 0 ≤ mTypeToR m v
 Proof
-  Cases_on `e` \\ EVAL_TAC \\ fs[]
-  (* \\ qspec_then `n0` (fn thm => once_rewrite_tac [GSYM thm]) POW_ONE *)
-  (* \\ rewrite_tac [GSYM REAL_POW_DIV] *)
-  (* \\ irule POW_POS \\ fs[] *)
-QED
-
-Theorem computeError_up:
-  ! v a b m.
-    abs v <= maxAbs (a,b) ==>
-    computeError v m <= computeError (maxAbs (a,b)) m
-Proof
-  rpt strip_tac \\ Cases_on `m`
-  \\ fs[mTypeToR_def, computeError_def, maxAbs_def]
-  \\ irule REAL_LE_TRANS \\ asm_exists_tac
-  \\ fs[ABS_LE]
+  Cases_on `m` \\ fs[mTypeToR_def]
+  \\ rpt strip_tac \\ COND_CASES_TAC \\ fs[]
 QED
 
 Theorem computeError_pos:
@@ -105,8 +132,48 @@ Theorem computeError_pos:
 Proof
   fs[computeError_def]
   \\ Cases_on ‘m’ \\ fs[mTypeToR_def]
-  \\ irule REAL_LE_ADD \\ fs[maxAbs_def, max_def]
-  \\ TOP_CASE_TAC \\ fs[minValue_pos_def]
+  \\ COND_CASES_TAC \\ fs[]
+QED
+
+Theorem normal_not_denormal:
+  normal x m ⇒ ~ denormal x m
+Proof
+  fs[normal_def, denormal_def] \\ Cases_on ‘m’ \\ fs[REAL_NOT_LT]
+QED
+
+(** Show that the roundoff error contributed by normal numbers is worse than the
+    roundoff error contributed of any possible denormal number.
+    This is the key theorem for integrating denormals with the 1+delta model *)
+Theorem computeError_up:
+  ∀ v a b m.
+    abs v <= maxAbs (a,b) ⇒
+    (* normal (maxAbs(a,b)) m ⇒ *)
+    computeError v m <= computeError (maxAbs (a,b)) m
+Proof
+  rpt strip_tac \\ Cases_on ‘m’
+  \\ fs[Excl "RMUL_LEQNORM", mTypeToR_def, computeError_def, maxAbs_def]
+  \\ imp_res_tac normal_not_denormal
+  \\ fs[Excl "RMUL_LEQNORM"]
+  \\ ntac 2 COND_CASES_TAC
+  \\ TRY (irule REAL_LE_LMUL_IMP \\ fs[] \\ irule REAL_LE_TRANS
+          \\ qexists_tac ‘max (abs a) (abs b)’ \\ fs[ABS_LE])
+  \\ fs[Excl "RMUL_LEQNORM", normal_def]
+  \\ TRY (fs[Excl "RMUL_LEQNORM", denormal_def, REAL_NOT_LT]
+          \\ ‘abs (max (abs a) (abs b)) = max (abs a) (abs b)’
+            by (fs[REAL_LE_MAX])
+          \\ ‘max (abs a) (abs b) < max (abs a) (abs b)’ suffices_by (fs[])
+          \\ RealArith.REAL_ASM_ARITH_TAC)
+  \\ fs[Excl "RMUL_LEQNORM", denormal_def]
+  \\ TRY (‘abs v ≤ 0’ suffices_by (fs[])
+          \\ pop_assum (rewrite_tac o single o GSYM)
+          \\ fs[] \\ NO_TAC)
+  \\ irule REAL_LE_TRANS
+  THENL (map qexists_tac
+         [‘1 / 2048 * minValue_pos M16’,
+          ‘1 / 16777216 * minValue_pos M32’,
+          ‘1 / 9007199254740992 * minValue_pos M64’])
+  \\ conj_tac \\ TRY (irule REAL_LE_LMUL_IMP \\ fs[REAL_NOT_LT])
+  \\ fs[minExponentPos_def, minValue_pos_def]
 QED
 
 (**
@@ -126,7 +193,12 @@ Definition isMorePrecise_def:
     | F w1 f1, F w2 f2 => w2 <= w1
     | F w f, _ => F
     | _, F w f => F
-    | _, _ => (mTypeToR (m1) <= mTypeToR (m2))
+    | M16, M16 => T
+    | M32, M32 => T
+    | M32, M16 => T
+    | M64, REAL => F
+    | M64, _ => T
+    | _, _ => F
 End
 
 (**
@@ -235,21 +307,6 @@ Definition isJoin3_def:
        |SOME mNew => mNew = mj)
 End
 
-(** Goldberg - Handbook of Floating-Point Arithmetic: (p.183)
-  𝛃^(e_min -p + 1) = 𝛃^(e_min -1) for base 2
-**)
-Definition normal_def:
-  normal (v:real) (m:mType) =
-  (minValue_pos m <= abs v /\ abs v <= maxValue m)
-End
-
-Definition denormal_def:
-  denormal (v:real) (m:mType) =
-    case m of
-      | REAL => F
-      | _ => ((abs v) < (minValue_pos m) /\ v <> 0)
-End
-
 (**
   Predicate that is true if and only if the given value v is a valid
   floating-point value according to the the type m.