diff --git a/hol4/ErrorBoundsScript.sml b/hol4/ErrorBoundsScript.sml
index f9b2a3ac18cf192bc65c660597a555cc165986fa..de5066aff8284f56497f1d83fabaabf983a454d3 100644
--- a/hol4/ErrorBoundsScript.sml
+++ b/hol4/ErrorBoundsScript.sml
@@ -440,6 +440,42 @@ Proof
\\ rewrite_tac [ABS_MUL] \\ irule REAL_LE_LMUL_IMP \\ fs[]
QED
+Theorem sqrt_abs_err_bounded:
+ !(e:real expr) (nR:real) nR1 (nF1:real) (nF:real) (E1:env) (E2:env) (err:real)
+ (m1 m:mType) (Gamma: real expr -> mType option).
+ eval_expr E1 (toRTMap Gamma) (toREval e) nR1 REAL /\
+ eval_expr E2 Gamma e nF1 m1 /\
+ eval_expr E1 (toRTMap Gamma) (toREval (Unop Sqrt e)) nR REAL /\
+ eval_expr (updEnv 1 nF1 emptyEnv)
+ (updDefVars (Unop Sqrt (Var 1)) m
+ (updDefVars (Var 1) m1 Gamma))
+ (Unop Sqrt (Var 1)) nF m /\
+ abs (nR1 - nF1) <= err ⇒
+ abs (nR - nF) <= abs (sqrt nR1 - sqrt nF1) + computeError (sqrt nF1) m
+Proof
+ rpt strip_tac
+ \\ inversion `eval_expr (updEnv _ _ _) _ _ _ _` eval_expr_cases \\ gs[]
+ \\ inversion `eval_expr (updEnv _ _ _) _ _ _ _` eval_expr_cases \\ gs[]
+ \\ rveq \\ gs[toREval_def, updDefVars_def] \\ rveq
+ \\ inversion `eval_expr _ _ (Unop Sqrt _) _ _` eval_expr_cases \\ gs[]
+ \\ rveq \\ gs[]
+ \\ ‘m1' = REAL’ by (Cases_on ‘m1'’ \\ gs[isCompat_def])
+ \\ rveq \\ gs[perturb_def]
+ \\ ‘v1 = nR1’ by (imp_res_tac meps_0_deterministic \\ gs[])
+ \\ rveq
+ \\ Cases_on ‘m’ \\ gs[perturb_def, evalUnop_def, computeError_pos, real_sub]
+ \\ TRY (COND_CASES_TAC \\ gs[])
+ \\ rewrite_tac [REAL_NEG_ADD, REAL_ADD_ASSOC, REAL_ADD_LDISTRIB]
+ \\ triangle_tac \\ gs[computeError_def, REAL_ABS_MUL]
+ \\ irule REAL_LE_TRANS THENL [
+ qexists_tac ‘mTypeToR M16 (sqrt nF1) * abs (sqrt nF1)’,
+ qexists_tac ‘mTypeToR M32 (sqrt nF1) * abs (sqrt nF1)’,
+ qexists_tac ‘mTypeToR M64 (sqrt nF1) * abs (sqrt nF1)’]
+ \\ conj_tac
+ \\ TRY (irule REAL_LE_RMUL_IMP \\ gs[])
+ \\ gs[REAL_MUL_COMM, mTypeToR_def]
+QED
+
Theorem nonzerop_EQ1_I'[simp]:
0 < r ⇒ (nonzerop r = 1)
Proof
diff --git a/hol4/ErrorIntervalInferenceScript.sml b/hol4/ErrorIntervalInferenceScript.sml
index 32574f07357013014aac692274efc589ae62245d..b140f267be9c32ef7230d96e598f4255c3aacd04 100644
--- a/hol4/ErrorIntervalInferenceScript.sml
+++ b/hol4/ErrorIntervalInferenceScript.sml
@@ -37,6 +37,22 @@ Definition inferErrorbound_def:
SOME (FloverMapTree_insert e (iv_f, err1) recRes);
od
| Unop Inv e1 => NONE
+ | Unop Sqrt e1 =>
+ do
+ recRes <- inferErrorbound e1 typeMap I akk;
+ (iv1, err1) <- FloverMapTree_find e1 recRes;
+ let errIve1 = widenInterval iv1 err1;
+ sqrtRealLo = newton ITERCOUNT (IVlo iv1 * 0.99) (IVlo iv1 * 0.99);
+ sqrtRealHi = newton ITERCOUNT (IVhi iv1 * 1.01) (IVhi iv1 * 1.01);
+ sqrtRealIv = (sqrtRealLo, sqrtRealHi);
+ sqrtFinLo= newton ITERCOUNT (IVlo errIve1 * 0.99) (IVlo errIve1 * 0.99);
+ sqrtFinHi = newton ITERCOUNT (IVhi errIve1 * 1.01) (IVhi errIve1 * 1.01);
+ sqrtFinIv = (sqrtFinLo, sqrtFinHi);
+ propError = maxAbs (subtractInterval sqrtRealIv sqrtFinIv);
+ mAbs = maxAbs sqrtFinIv
+ in
+ SOME (FloverMapTree_insert e (iv_f, propError + computeError mAbs m) recRes);
+ od
| Binop op e1 e2 =>
do
recRes1 <- inferErrorbound e1 typeMap I akk;
diff --git a/hol4/ErrorValidationScript.sml b/hol4/ErrorValidationScript.sml
index 5ad5d16b06967e07fc462ac387dbbba9773135fe..fbb521f9c23ef560c7859f5481cff9e5e5cd1967 100644
--- a/hol4/ErrorValidationScript.sml
+++ b/hol4/ErrorValidationScript.sml
@@ -11,7 +11,7 @@ open AbbrevsTheory ExpressionsTheory ExpressionSemanticsTheory RealSimpsTheory
RealRangeArithTheory FloverTactics MachineTypeTheory
ExpressionAbbrevsTheory ErrorBoundsTheory IntervalArithTheory
TypeValidatorTheory IntervalValidationTheory EnvironmentsTheory
- CommandsTheory ssaPrgsTheory FloverMapTheory;
+ CommandsTheory ssaPrgsTheory FloverMapTheory sqrtApproxTheory;
open preambleFloVer;
val _ = new_theory "ErrorValidation";
@@ -29,26 +29,51 @@ val triangle_tac =
fun real_rewrite t =
rewrite_tac [REAL_ARITH (Parse.Term t)];
+val real_prove =
+ rpt (qpat_x_assum `!x. _` kall_tac)
+ \\ rpt (qpat_x_assum `_ ==> ! x. _` kall_tac)
+ \\ REAL_ASM_ARITH_TAC;
+
Definition validErrorbound_def:
validErrorbound e (typeMap: typeMap) (A:analysisResult) (dVars:num_set)=
case FloverMapTree_find e A, FloverMapTree_find e typeMap of
| NONE, _ => F
| _, NONE => F
| SOME (intv, err), SOME m =>
- (if (0 <= err)
- then
+ (if (0 <= err) then
case e of
| Var v =>
- if (lookup v dVars = SOME ())
- then T
+ if (lookup v dVars = SOME ()) then T
else computeError (maxAbs intv) m <= err
| Const _ n => (computeError (maxAbs intv) m <= err)
| Unop Neg e1 =>
- if (validErrorbound e1 typeMap A dVars)
- then
+ if (validErrorbound e1 typeMap A dVars) then
case (FloverMapTree_find e1 A) of
- | SOME (_, err1) => err = err1
- | _ => F
+ | SOME (_, err1) => err = err1
+ | _ => F
+ else F
+ | Unop Sqrt e1 =>
+ if (validErrorbound e1 typeMap A dVars) then
+ case (FloverMapTree_find e1 A) of
+ | SOME (iv1, err1) =>
+ let errIve1 = widenInterval iv1 err1;
+ sqrtRealLo = newton ITERCOUNT (IVlo iv1 * 0.99) (IVlo iv1 * 0.99);
+ sqrtRealHi = newton ITERCOUNT (IVhi iv1 * 1.01) (IVhi iv1 * 1.01);
+ sqrtRealIv = (sqrtRealLo, sqrtRealHi);
+ sqrtFinLo= newton ITERCOUNT (IVlo errIve1 * 0.99) (IVlo errIve1 * 0.99);
+ sqrtFinHi = newton ITERCOUNT (IVhi errIve1 * 1.01) (IVhi errIve1 * 1.01);
+ sqrtFinIv = (sqrtFinLo, sqrtFinHi);
+ propError = maxAbs (subtractInterval sqrtRealIv sqrtFinIv);
+ mAbs = maxAbs sqrtFinIv
+ in
+ if (validate_newton_down sqrtRealLo (IVlo iv1) ∧
+ validate_newton_up sqrtRealHi (IVhi iv1) ∧
+ validate_newton_down sqrtFinLo (IVlo errIve1) ∧
+ validate_newton_up sqrtFinHi (IVhi errIve1) ∧
+ 0 < IVlo errIve1) then
+ propError + computeError mAbs m <= err
+ else F
+ | _ => F
else F
| Unop Inv e1 => F
| Binop op e1 e2 =>
@@ -2232,6 +2257,80 @@ Proof
\\ rpt (first_x_assum (destruct) \\ fs[])
QED
+Theorem validErrorboundCorrectSqrt:
+ ∀(E1 E2:env) (A:analysisResult) (e:real expr)
+ (nR nR1 nF nF1:real) (err err1:real) (alo ahi elo ehi:real) dVars m m1 Gamma.
+ eval_Real E1 Gamma e nR1 ∧
+ eval_Fin E2 Gamma e nF1 m1 ∧
+ eval_Real E1 Gamma (Unop Sqrt e) nR ∧
+ eval_expr (updEnv 1 nF1 emptyEnv)
+ (updDefVars (Unop Sqrt (Var 1)) m
+ (updDefVars (Var 1) m1 (toRExpMap Gamma)))
+ (Unop Sqrt (Var 1)) nF m /\
+ validErrorbound (Unop Sqrt e) Gamma A dVars ∧
+ elo <= nR1 ∧ nR1 <= ehi ∧
+ 0 < elo ∧
+ FloverMapTree_find e A = SOME ((elo, ehi), err1) ∧
+ FloverMapTree_find (Unop Sqrt e) A = SOME ((alo, ahi), err) ∧
+ FloverMapTree_find (Unop Sqrt e) Gamma = SOME m ∧
+ abs (nR1 - nF1) <= err1 ⇒
+ abs (nR - nF) <= err
+Proof
+ fs [toREval_def] \\ rpt strip_tac
+ \\ fs[Once validErrorbound_eq] \\ rveq
+ \\ irule REAL_LE_TRANS
+ \\ qexists_tac ‘abs (sqrt nR1 - sqrt nF1) + computeError (sqrt nF1) m’
+ \\ conj_tac
+ >- (drule sqrt_abs_err_bounded \\ rpt $ disch_then drule
+ \\ disch_then $ qspecl_then [‘nR’, ‘nF’, ‘err1’, ‘m’] mp_tac
+ \\ gs[toREval_def])
+ \\ ‘contained nF1 (widenInterval (elo, ehi) err1)’
+ by (irule distance_gives_iv \\ find_exists_tac \\ gs[contained_def])
+ \\ rpt (
+ qpat_x_assum `validate_newton_down _ _` $ mp_then Any mp_tac validate_newton_lb_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
+ \\ real_prove)
+ \\ strip_tac)
+ \\ ‘0 < ehi’
+ by (
+ irule REAL_LTE_TRANS \\ qexists_tac ‘elo’ \\ gs[]
+ \\ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[])
+ \\ ‘0 < SND (widenInterval (elo, ehi) err1)’
+ by (
+ irule REAL_LTE_TRANS \\ qexists_tac ‘FST (widenInterval (elo, ehi) err1)’ \\ gs[widenInterval_def]
+ \\ REAL_ASM_ARITH_TAC)
+ \\ rpt (
+ qpat_x_assum `validate_newton_up _ _` $ mp_then Any mp_tac validate_newton_ub_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
+ \\ real_prove)
+ \\ strip_tac)
+ \\ irule REAL_LE_TRANS \\ once_rewrite_tac[CONJ_COMM] \\ asm_exists_tac \\ gs[]
+ \\ irule REAL_LE_ADD2 \\ conj_tac
+ (* Prove propagation error *)
+ >- (
+ irule contained_leq_maxAbs
+ \\ assume_tac $ SIMP_RULE std_ss [validIntervalSub_def] interval_subtraction_valid
+ \\ pop_assum irule \\ conj_tac \\ gs[contained_def]
+ \\ conj_tac \\ irule REAL_LE_TRANS
+ >- (asm_exists_tac \\ gs[] \\ irule transcTheory.SQRT_MONO_LE \\ real_prove)
+ >- (once_rewrite_tac[CONJ_COMM] \\ asm_exists_tac \\ gs[]
+ \\ irule transcTheory.SQRT_MONO_LE \\ gs[] \\ REAL_ASM_ARITH_TAC)
+ >- (asm_exists_tac \\ gs[] \\ irule transcTheory.SQRT_MONO_LE \\ real_prove)
+ \\ once_rewrite_tac[CONJ_COMM] \\ asm_exists_tac \\ gs[]
+ \\ irule transcTheory.SQRT_MONO_LE \\ gs[] \\ REAL_ASM_ARITH_TAC)
+ (* New error *)
+ \\ irule computeError_up
+ \\ irule contained_leq_maxAbs
+ \\ gs[contained_def] \\ conj_tac \\ irule REAL_LE_TRANS
+ >- (asm_exists_tac \\ gs[] \\ irule transcTheory.SQRT_MONO_LE \\ real_prove)
+ \\ once_rewrite_tac[CONJ_COMM] \\ asm_exists_tac \\ gs[]
+ \\ irule transcTheory.SQRT_MONO_LE \\ gs[] \\ REAL_ASM_ARITH_TAC
+QED
+
(**
Soundness theorem for the error bound validator.
Proof is by induction on the expression e.
@@ -2277,31 +2376,79 @@ Proof
\\ qpat_x_assum `validErrorbound _ _ _ _` mp_tac
\\ simp[Once validErrorbound_eq] \\ strip_tac
\\ Cases_on `u` \\ fs[toREval_def] \\ rveq
+ >- (
+ rw_thm_asm `(domain _) DIFF _ SUBSET _` usedVars_def
+ \\ inversion `eval_expr E1 _ _ _ _` eval_expr_cases \\ fs[]
+ \\ `m' = REAL` by (irule toRTMap_eval_REAL \\ find_exists_tac \\ fs[])
+ \\ PairCases_on `v4` \\ fs[] \\ rveq
+ \\ rename1 `FloverMapTree_find e A = SOME ((e_lo, e_hi) , err)`
+ \\ once_rewrite_tac [eval_expr_cases] \\ fs[]
+ \\ first_x_assum
+ (qspecl_then
+ [`E1`, `E2`, `A`, `v1`, `err`, `e_lo`, `e_hi`, `fVars`,
+ `dVars`, `Gamma`]
+ destruct)
+ >- (fs[Once validTypes_def, Once validRanges_def])
+ \\ fs[Once validTypes_def, toRExpMap_def] \\ conj_tac \\ rpt strip_tac
+ >- (qexistsl_tac [`me`, `nF`] \\ fs[]
+ \\ `m' = me`
+ by (irule validTypes_exec
+ \\ rpt (find_exists_tac \\ fs[toRExpMap_def]))
+ \\ rveq \\ fs[])
+ \\ fs[evalUnop_def]
+ \\ once_rewrite_tac [real_sub]
+ \\ rewrite_tac [GSYM REAL_NEG_ADD, ABS_NEG, GSYM real_sub]
+ \\ first_x_assum drule
+ \\ rpt (disch_then drule)
+ \\ disch_then assume_tac \\ fs[]
+ \\ first_x_assum irule \\ asm_exists_tac \\ fs[])
\\ rw_thm_asm `(domain _) DIFF _ SUBSET _` usedVars_def
\\ inversion `eval_expr E1 _ _ _ _` eval_expr_cases \\ fs[]
- \\ `m' = REAL` by (irule toRTMap_eval_REAL \\ find_exists_tac \\ fs[])
- \\ PairCases_on `v4` \\ fs[] \\ rveq
- \\ rename1 `FloverMapTree_find e A = SOME ((e_lo, e_hi) , err)`
+ \\ `m1 = REAL` by (irule toRTMap_eval_REAL \\ find_exists_tac \\ fs[])
+ \\ PairCases_on `iv1` \\ fs[] \\ rveq
+ \\ rename1 `FloverMapTree_find e A = SOME ((e_lo, e_hi) , err1)`
\\ once_rewrite_tac [eval_expr_cases] \\ fs[]
\\ first_x_assum
- (qspecl_then
- [`E1`, `E2`, `A`, `v1`, `err`, `e_lo`, `e_hi`, `fVars`,
- `dVars`, `Gamma`]
- destruct)
- >- (fs[Once validTypes_def, Once validRanges_def])
- \\ fs[Once validTypes_def, toRExpMap_def] \\ conj_tac \\ rpt strip_tac
- >- (qexistsl_tac [`me`, `nF`] \\ fs[]
+ (qspecl_then
+ [`E1`, `E2`, `A`, `v1`, `err1`, `e_lo`, `e_hi`, `fVars`,
+ `dVars`, `Gamma`]
+ destruct)
+ >- (fs[Once validTypes_def, Once validRanges_def])
+ \\ qpat_x_assum `validRanges _ _ _ _` mp_tac
+ \\ qpat_x_assum `validTypes _ _` mp_tac
+ \\ simp[Once validRanges_def, Once validTypes_def]
+ \\ rpt $ disch_then strip_assume_tac
+ \\ imp_res_tac validRanges_single
+ \\ imp_res_tac meps_0_deterministic \\ rveq
+ \\ gs[] \\ rveq \\ gs[]
+ \\ fs[toRExpMap_def] \\ conj_tac \\ rpt strip_tac
+ >- (qexistsl_tac [`me`, `nF`, ‘0’] \\ fs[]
\\ `m' = me`
by (irule validTypes_exec
\\ rpt (find_exists_tac \\ fs[toRExpMap_def]))
- \\ rveq \\ fs[])
- \\ fs[evalUnop_def]
- \\ once_rewrite_tac [real_sub]
- \\ rewrite_tac [GSYM REAL_NEG_ADD, ABS_NEG, GSYM real_sub]
- \\ first_x_assum drule
- \\ rpt (disch_then drule)
- \\ disch_then assume_tac \\ fs[]
- \\ first_x_assum irule \\ asm_exists_tac \\ fs[])
+ \\ rveq \\ fs[] \\ conj_tac
+ >- (irule mTypeToR_pos)
+ \\ res_tac \\ first_x_assum $ mp_then Any mp_tac distance_gives_iv
+ \\ gs[contained_def]
+ \\ disch_then $ qspec_then ‘(e_lo, e_hi)’ mp_tac \\ gs[]
+ \\ strip_tac \\ irule REAL_LTE_TRANS
+ \\ qexists_tac ‘FST (widenInterval (e_lo, e_hi) err')’ \\ gs[])
+ \\ rveq \\ res_tac
+ \\ irule validErrorboundCorrectSqrt
+ \\ qexistsl_tac [‘A’, ‘E1’, ‘E2’, ‘Gamma’, ‘ehi’, ‘elo’, ‘dVars’, ‘e’,
+ ‘e_hi’, ‘e_lo’, ‘err'’, ‘m’, ‘m1’, ‘v1'’, ‘v1’]
+ \\ gs[toRExpMap_def] \\ rpt conj_tac
+ >- (gs[perturb_def, toRExpMap_def, toREval_def]
+ \\ ‘vR = evalUnop Sqrt v1’
+ by (
+ inversion `eval_expr E1 _ (Unop Sqrt _) _ _` eval_expr_cases \\ fs[]
+ \\ Cases_on ‘m1'’ \\ gs[isCompat_def, perturb_def]
+ \\ imp_res_tac meps_0_deterministic \\ gs[])
+ \\ rveq \\ gs[])
+ >- simp[Once validErrorbound_def]
+ \\ irule Unop_sqrt' \\ fsrw_tac [SATISFY_ss] [updDefVars_def]
+ \\ qexistsl_tac [‘delta'’, ‘m1’, ‘v1'’] \\ gs[]
+ \\ simp[Once eval_expr_cases] \\ gs[updDefVars_def])
>- (rename1 `Binop op e1 e2` \\ fs[]
\\ qpat_x_assum `validErrorbound _ _ _ _` mp_tac
\\ simp[Once validErrorbound_eq] \\ strip_tac
diff --git a/hol4/IEEE_connectionScript.sml b/hol4/IEEE_connectionScript.sml
index 376f9791b94728827851e1141d10ec78f6f8af94..8428a2cf8fbbb0856ab26c829be52e7fa6c1926e 100644
--- a/hol4/IEEE_connectionScript.sml
+++ b/hol4/IEEE_connectionScript.sml
@@ -33,8 +33,12 @@ Definition eval_expr_float_def:
(eval_expr_float (Unop Neg e) E =
case eval_expr_float e E of
| SOME v => SOME (fp64_negate v)
- | _ => NONE) /\
- (eval_expr_float (Unop Inv e) E = NONE) /\
+ | _ => NONE) ∧
+ (eval_expr_float (Unop Inv e) E = NONE) ∧
+ (eval_expr_float (Unop Sqrt e) E =
+ case eval_expr_float e E of
+ | SOME v => SOME (fp64_sqrt dmode v)
+ | _ => NONE) ∧
(eval_expr_float (Binop b e1 e2) E =
(case (eval_expr_float e1 E), (eval_expr_float e2 E) of
| SOME v1, SOME v2 =>
@@ -115,6 +119,26 @@ Proof
\\ fs[REAL_MUL_RINV]
QED
+Theorem zero_lt_sign_zero:
+ 0 < float_to_real fp ⇒
+ fp.Sign = 0w
+Proof
+ rpt strip_tac \\ CCONTR_TAC
+ \\ ‘fp.Sign = 1w’ by (Cases_on ‘fp.Sign’ \\ gs[])
+ \\ gs[float_to_real_def]
+ \\ ‘0 < 0:real’ suffices_by gs[]
+ \\ irule REAL_LTE_TRANS
+ \\ asm_exists_tac \\ gs[REAL_MUL_LNEG]
+ \\ Cases_on ‘fp.Exponent = 0w’ \\ gs[real_div]
+ >- (
+ irule realTheory.REAL_LE_MUL \\ conj_tac
+ \\ irule realTheory.REAL_LE_MUL \\ gs[] )
+ \\ gs[real_div] \\ irule REAL_LE_MUL \\ conj_tac
+ \\ TRY (irule REAL_LE_MUL \\ gs[])
+ \\ irule REAL_LE_ADD \\ conj_tac \\ gs[]
+ \\ irule REAL_LE_MUL \\ conj_tac \\ gs[]
+QED
+
Theorem pow_simp1[local] = Q.prove (`2 pow 2047 / 2 pow 1023 = 2 pow 1024`,
qspecl_then [`2`, `2047`,`1023`] destruct real_div_pow \\ fs[]);
@@ -529,7 +553,6 @@ Theorem eval_expr_gives_IEEE:
is64BitEval (toRExp e) /\
is64BitEnv Gamma /\
noDowncast (toRExp e) /\
- (* eval_expr_valid e E2 /\ *)
(∀v.
v ∈ domain dVars ⇒
∃vF m.
@@ -604,6 +627,189 @@ Proof
qpat_x_assum ‘validErrorbound _ _ _ _’
(fn thm => mp_tac (ONCE_REWRITE_RULE [validErrorbound_def] thm))
\\ fs[option_case_eq] \\ rpt (TOP_CASE_TAC \\ fs[]))
+ >- (
+ fs[Once validTypes_def] \\ rveq
+ \\ imp_res_tac validTypes_single
+ \\ ‘M64 = m1’
+ by (qpat_x_assum ‘isCompat m1 _’ mp_tac \\ Cases_on ‘m1’
+ \\ simp[isCompat_def])
+ \\ rveq
+ \\ ‘M64 = mG’
+ by (qpat_x_assum ‘toRExpMap _ _ = SOME _’ mp_tac
+ \\ qpat_x_assum ‘FloVerMapTree_find _ _ = SOME mG’ mp_tac
+ \\ simp[toRExpMap_def])
+ \\ ‘me = M64’
+ by (qpat_x_assum ‘isCompat me _’ mp_tac \\ rveq \\ Cases_on ‘me’
+ \\ simp[isCompat_def])
+ \\ rveq \\ fs[isCompat_def] \\ rpt (qpat_x_assum ‘T’ kall_tac)
+ \\ rveq
+ \\ fs[eval_expr_float_def]
+ \\ first_x_assum
+ (qspecl_then
+ [`E1`, `E2`, `E2_real`, `Gamma`, `v1`, `A`, `fVars`, `dVars`]
+ destruct)
+ >- (
+ fs[] \\ rpt conj_tac
+ >- (fs[Once validTypes_def])
+ >- (rveq \\ rpt strip_tac
+ \\ drule validTypes_single \\ strip_tac \\ rfs[] \\ rveq
+ \\ first_x_assum irule
+ \\ find_exists_tac \\ fs[]
+ \\ asm_exists_tac \\ fs[])
+ >- (fs[Once validRanges_def])
+ >- (
+ qpat_x_assum ‘validErrorbound _ _ _ _’
+ (fn thm => assume_tac (ONCE_REWRITE_RULE [validErrorbound_def] thm))
+ \\ fs[option_case_eq]
+ \\ pop_assum mp_tac \\ rpt (TOP_CASE_TAC \\ fs[]))
+ >- (
+ rw_thm_asm `FPRangeValidator _ _ _ _` FPRangeValidator_def
+ \\ fs[]
+ \\ Cases_on `FloverMapTree_find (Unop Sqrt (toRExp e)) A` \\ fs[]
+ \\ PairCases_on `x` \\ fs[]
+ \\ rw_asm_star `FloverMapTree_find (Unop _ _) Gamma = _`)
+ >- (rw_thm_asm `domain (usedVars _) DIFF _ SUBSET _` usedVars_def \\ fs[])
+ \\ rw_thm_asm `is64BitEval _` is64BitEval_def \\ fs[])
+ \\ fs[fp64_sqrt_def, fp64_to_float_float_to_fp64]
+ \\ qpat_x_assum `validRanges _ _ _ _` mp_tac
+ \\ simp[Once validRanges_def] \\ rpt strip_tac
+ \\ imp_res_tac validRanges_single
+ \\ rename1 ‘FloverMapTree_find (toRExp e) A = SOME (iv1, err1)’
+ \\ ‘0 < IVlo iv1’ by (res_tac \\ gs[IVlo_def])
+ \\ rename1 ‘IVlo iv1 ≤ vR1’
+ (* Obtain evaluation for E2_real*)
+ \\ ‘!vF1 m1. eval_expr E2_real (toRExpMap Gamma) (toRExp e) vF1 m1 ==>
+ abs (vR1 - vF1) <= err1’
+ by (qspecl_then [`toRExp e`, `E1`, `E2_real`, `A`,`vR1`,
+ `err1`, `IVlo iv1`, `IVhi iv1`, `fVars`,
+ `dVars`,`Gamma`] destruct validErrorbound_sound
+ \\ rpt conj_tac \\ fs[]
+ >- (fs [DIFF_DEF, SUBSET_DEF]
+ \\ rpt strip_tac \\ first_x_assum irule
+ \\ once_rewrite_tac [usedVars_def] \\ fs[domain_union])
+ \\ qpat_x_assum ‘validErrorbound _ _ _ _’
+ (fn thm => mp_tac (ONCE_REWRITE_RULE [validErrorbound_def] thm))
+ \\ fs[option_case_eq] \\ rpt (TOP_CASE_TAC \\ fs[]))
+ \\ ‘validFloatValue (float_to_real (fp64_to_float v)) M64’
+ by (drule FPRangeValidator_sound
+ \\ disch_then $ qspecl_then [‘toRExp e’, ‘fp64_to_real v’, ‘M64’] mp_tac
+ \\ gs[] \\ impl_tac
+ >- (rpt conj_tac
+ >- (drule eval_eq_env
+ \\ rpt $ disch_then drule \\ gs[fp64_to_real_def])
+ >- (
+ qpat_x_assum ‘validErrorbound _ _ _ _’
+ (fn thm => mp_tac (ONCE_REWRITE_RULE [validErrorbound_def] thm))
+ \\ fs[option_case_eq] \\ rpt (TOP_CASE_TAC \\ fs[]))
+ >- (
+ qpat_x_assum ‘FPRangeValidator _ _ _ _’
+ (fn thm => mp_tac (ONCE_REWRITE_RULE [FPRangeValidator_def] thm))
+ \\ fs[option_case_eq] \\ rpt (TOP_CASE_TAC \\ fs[]))
+ \\ fs [DIFF_DEF, SUBSET_DEF]
+ \\ rpt strip_tac \\ first_x_assum irule
+ \\ once_rewrite_tac [usedVars_def] \\ fs[domain_union])
+ \\ simp[fp64_to_real_def])
+ \\ ‘contained (float_to_real (fp64_to_float v))
+ (widenInterval (FST iv1, SND iv1) err1)’
+ by (
+ irule distance_gives_iv
+ \\ qexists_tac `vR1` \\ fs [contained_def]
+ \\ first_x_assum irule
+ \\ qexists_tac `M64`
+ \\ drule eval_eq_env
+ \\ rpt (disch_then drule) \\ fs[])
+ \\ ‘0 < FST (widenInterval (FST iv1, SND iv1) err1)’
+ by (
+ qpat_x_assum ‘validErrorbound _ _ _ _’
+ (fn thm => mp_tac (ONCE_REWRITE_RULE [validErrorbound_def] thm))
+ \\ fs[option_case_eq] \\ rpt (TOP_CASE_TAC \\ fs[]))
+ \\ ‘0 < float_to_real (fp64_to_float v)’
+ by (gs[contained_def, widenInterval_def] \\ irule REAL_LTE_TRANS
+ \\ asm_exists_tac \\ gs[])
+ \\ ‘(fp64_to_float v).Sign = 0w’
+ by imp_res_tac zero_lt_sign_zero
+ \\ ‘validFloatValue (evalUnop Sqrt (float_to_real (fp64_to_float v))) M64’
+ by (
+ drule FPRangeValidator_sound
+ \\ disch_then
+ (qspecl_then
+ [`Unop Sqrt (toRExp e)`,
+ `evalUnop Sqrt (float_to_real (fp64_to_float v))`, `M64`]
+ irule)
+ \\ fs[]
+ \\ qexists_tac ‘e’ \\ fs[]
+ \\ rpt conj_tac
+ >- (simp[Once validTypes_def, isCompat_def] \\ first_x_assum MATCH_ACCEPT_TAC)
+ >- (simp[Once validRanges_def] \\ asm_exists_tac \\ fs[]
+ \\ fs[] \\ rveq \\ fs[])
+ \\ irule eval_eq_env
+ \\ qexists_tac ‘toREnv E2’ \\ rpt conj_tac >- fs[]
+ \\ irule Unop_sqrt'
+ \\ qexistsl_tac [‘0’, `M64`, ‘M64’, `float_to_real (fp64_to_float v)`]
+ \\ fs[perturb_def, evalUnop_def, mTypeToR_pos, isCompat_def])
+ \\ qpat_x_assum `validFloatValue (evalUnop _ _) M64` (assume_tac o SIMP_RULE std_ss [validFloatValue_def])
+ \\ gs[]
+ (* normal sqrt *)
+ >- (
+ Q.ISPEC_THEN `(fp64_to_float v):(52,11) float`
+ impl_subgoal_tac
+ float_sqrt_relative
+ >- (rpt conj_tac
+ \\ fs[validFloatValue_def,
+ normalTranslatedValue_implies_finiteness,
+ denormalTranslatedValue_implies_finiteness,
+ normalValue_implies_normalization,
+ GSYM float_is_zero_to_real, float_is_finite, evalUnop_def,
+ normalizes_def])
+ \\ fs[dmode_def] \\ simp[Once eval_expr_cases]
+ \\ qexistsl_tac [`M64`, `float_to_real (fp64_to_float v)`, ‘e'’]
+ \\ fs[perturb_def, evalUnop_def]
+ \\ imp_res_tac normal_not_denormal
+ \\ fs[REAL_INV_1OVER, mTypeToR_def, isCompat_def])
+ (* denormal sqrt *)
+ >- (
+ Q.ISPEC_THEN `(fp64_to_float v):(52,11) float`
+ impl_subgoal_tac
+ float_sqrt_relative_denorm
+ >- (rpt conj_tac
+ >- fs[validFloatValue_def,
+ normalTranslatedValue_implies_finiteness,
+ denormalTranslatedValue_implies_finiteness,
+ normalValue_implies_normalization, float_is_finite,
+ GSYM float_is_zero_to_real, evalUnop_def]
+ >- fs[validFloatValue_def,
+ normalTranslatedValue_implies_finiteness,
+ denormalTranslatedValue_implies_finiteness,
+ normalValue_implies_normalization, float_is_finite,
+ GSYM float_is_zero_to_real, evalUnop_def]
+ >- (
+ fs[normalTranslatedValue_implies_finiteness,
+ denormalTranslatedValue_implies_finiteness,
+ normalValue_implies_normalization, float_is_finite,
+ GSYM float_is_zero_to_real, evalUnop_def, denormal_def, minValue_pos_def]
+ \\ rewrite_tac [INT_MAX_def, INT_MIN_def, dimindex_11, EVAL “2 ** (11 - 1) - 1 - 1â€]
+ \\ irule REAL_LT_TRANS
+ \\ asm_exists_tac \\ fs[GSYM REAL_OF_NUM_POW, minExponentPos_def]
+ \\ irule REAL_LET_TRANS \\ qexists_tac ‘1 * inv (2 pow 1022)’
+ \\ conj_tac
+ >- (rewrite_tac[real_div] \\ irule REAL_LT_RMUL_IMP \\ fs[])
+ \\ fs[])
+ >- (irule REAL_LT_TRANS \\ qexists_tac ‘maxValue M64’
+ \\ fs[threshold_64_bit_lt_maxValue, denormal_def]
+ \\ irule REAL_LTE_TRANS \\ qexists_tac ‘minValue_pos M64’
+ \\ fs[minValue_pos_def, maxValue_def, GSYM REAL_OF_NUM_POW, evalUnop_def]
+ \\ irule REAL_LE_TRANS \\ qexists_tac ‘1’ \\ conj_tac
+ >- (once_rewrite_tac[GSYM REAL_INV1] \\ irule REAL_INV_LE_ANTIMONO_IMPR
+ \\ fs[POW_2_LE1])
+ \\ fs[POW_2_LE1])
+ \\ fs[INT_MAX_def, INT_MIN_def, dimindex_11])
+ \\ gs[dmode_def] \\ simp[Once eval_expr_cases]
+ \\ qexistsl_tac [`M64`, `float_to_real (fp64_to_float v)`, ‘e'’]
+ \\ fs[perturb_def, evalUnop_def]
+ \\ fs[REAL_INV_1OVER, mTypeToR_def, isCompat_def, minExponentPos_def])
+ (* sqrt 0 *)
+ \\ ‘0 < sqrt (float_to_real (fp64_to_float v))’ by (irule transcTheory.SQRT_POS_LT \\ gs[])
+ \\ gs[evalUnop_def])
>- (
rename1 `Binop b (toRExp e1) (toRExp e2)` \\ rveq
\\ IMP_RES_TAC validRanges_single
@@ -1786,9 +1992,8 @@ Proof
\\ rpt conj_tac \\ fs[is64BitBstep_def]
QED
-val IEEE_connection_expr = store_thm (
- "IEEE_connection_expr",
- ``!(e:word64 expr) (A:analysisResult) (P:precond) E1 E2 defVars fVars Gamma.
+Theorem IEEE_connection_expr:
+ !(e:word64 expr) (A:analysisResult) (P:precond) E1 E2 defVars fVars Gamma.
is64BitEval (toRExp e) /\
is64BitEnv defVars /\
noDowncast (toRExp e) /\
@@ -1802,7 +2007,8 @@ val IEEE_connection_expr = store_thm (
eval_expr_float e E2 = SOME vF /\
eval_expr (toREnv E2) (toRExpMap Gamma) (toRExp e)
(float_to_real (fp64_to_float vF)) M64 /\
- abs (vR - (float_to_real (fp64_to_float vF))) <= err``,
+ abs (vR - (float_to_real (fp64_to_float vF))) <= err
+Proof
simp [CertificateChecker_def]
\\ rpt strip_tac
\\ Cases_on `getValidMap defVars (toRExp e) FloverMapTree_empty` \\ fs[]
@@ -1837,11 +2043,11 @@ val IEEE_connection_expr = store_thm (
\\ rveq \\ fs[])
\\ fs[is64BitEnv_def] \\ first_x_assum MATCH_ACCEPT_TAC)
\\ rpt (find_exists_tac \\ fs[])
- \\ first_x_assum irule \\ find_exists_tac \\ fs[]);
+ \\ first_x_assum irule \\ find_exists_tac \\ fs[]
+QED
-val IEEE_connection_cmds = store_thm (
- "IEEE_connection_cmds",
- ``! (f:word64 cmd) (A:analysisResult) (P:precond) E1 E2 defVars (fVars:num_set)
+Theorem IEEE_connection_cmds:
+ ! (f:word64 cmd) (A:analysisResult) (P:precond) E1 E2 defVars (fVars:num_set)
Gamma.
is64BitBstep (toRCmd f) /\
is64BitEnv defVars /\
@@ -1856,7 +2062,8 @@ val IEEE_connection_cmds = store_thm (
bstep_float f E2 = SOME vF /\
bstep (toRCmd f) (toREnv E2) (toRExpMap Gamma)
(float_to_real (fp64_to_float vF)) M64 /\
- abs (vR - (float_to_real (fp64_to_float vF))) <= err``,
+ abs (vR - (float_to_real (fp64_to_float vF))) <= err
+Proof
simp [CertificateCheckerCmd_def]
\\ rpt strip_tac
\\ Cases_on `getValidMapCmd defVars (toRCmd f) FloverMapTree_empty` \\ fs[]
@@ -1908,63 +2115,7 @@ val IEEE_connection_cmds = store_thm (
\\ disch_then
(qspecl_then
[`outVars`, `vR`, `vF2`, `FST iv_e`, `SND iv_e`, `err_e`, `m2`] destruct)
- \\ fs[]);
-
-(** UNUSED:
-val normal_or_zero_def = Define `
- normal_or_zero (v:real) =
- (minValue_pos M64 <= abs v \/ v = 0)`;
-
-val isValid_def = Define `
- isValid e =
- let trans_e = optionLift e (\ v. SOME (float_to_real (fp64_to_float v))) NONE in
- optionLift trans_e normal_or_zero F`;
-
-val eval_expr_valid_def = Define `
- (eval_expr_valid (Var n) E = T) /\
- (eval_expr_valid (Const m v) E = T) /\
- (eval_expr_valid (Unop u e) E = eval_expr_valid e E) /\
- (eval_expr_valid (Binop b e1 e2) E =
- (eval_expr_valid e1 E /\ eval_expr_valid e2 E /\
- let e1_res = eval_expr_float e1 E in
- let e2_res = eval_expr_float e2 E in
- optionLift (e1_res)
- (\ v1. let v1_real = float_to_real (fp64_to_float v1)
- in
- optionLift e2_res
- (\ v2.
- let v2_real = float_to_real (fp64_to_float v2)
- in
- normal_or_zero (evalBinop b v1_real v2_real))
- T)
- T)) /\
- (eval_expr_valid (Fma e1 e2 e3) E =
- (eval_expr_valid e1 E /\ eval_expr_valid e2 E /\ eval_expr_valid e3 E /\
- let e1_res = eval_expr_float e1 E in
- let e2_res = eval_expr_float e2 E in
- let e3_res = eval_expr_float e3 E in
- optionLift (e1_res)
- (\ v1. let v1_real = float_to_real (fp64_to_float v1)
- in
- optionLift e2_res
- (\ v2.
- let v2_real = float_to_real (fp64_to_float v2)
- in
- optionLift e3_res
- (\ v3. let v3_real = float_to_real (fp64_to_float v3) in
- normal_or_zero (evalFma v1_real v2_real v3_real))
- T)
- T)
- T)) /\
- (eval_expr_valid (Downcast m e) E = eval_expr_valid e E)`;
-
-val bstep_valid_def = Define `
- (bstep_valid (Let m x e g) E =
- (eval_expr_valid e E /\
- optionLift (eval_expr_float e E)
- (\v_e. bstep_valid g (updFlEnv x v_e E))
- T)) /\
- (bstep_valid (Ret e) E = eval_expr_valid e E)`;
-*)
+ \\ fs[]
+QED
val _ = export_theory ();
diff --git a/hol4/IntervalArithScript.sml b/hol4/IntervalArithScript.sml
index 9757d1ad4acc4a5818666b38ff99465f2ba051f8..3419d12a6a7c2a008da83c8ef66e6422c34915eb 100644
--- a/hol4/IntervalArithScript.sml
+++ b/hol4/IntervalArithScript.sml
@@ -92,6 +92,10 @@ Definition invertInterval_def:
invertInterval (iv:interval) = (1 /(IVhi iv), 1 /(IVlo iv))
End
+Definition sqrtInterval_def:
+ sqrtInterval (iv:interval) = (sqrt (IVlo iv), sqrt (IVhi iv))
+End
+
Definition addInterval_def:
addInterval (iv1:interval) (iv2:interval) = absIntvUpd (+) iv1 iv2
End
@@ -271,6 +275,14 @@ Proof
\\ asm_exists_tac \\ fs []
QED
+Theorem iv_sqrt_preserves_valid:
+ ∀ iv.
+ 0 ≤ IVlo iv ∧ valid iv ⇒ valid (sqrtInterval iv)
+Proof
+ gs[valid_def, sqrtInterval_def] \\ rpt strip_tac
+ \\ irule transcTheory.SQRT_MONO_LE \\ gs[]
+QED
+
val interval_addition_valid = store_thm ("interval_addition_valid",
``!iv1 iv2. validIntervalAdd iv1 iv2 (addInterval iv1 iv2)``,
fs iv_ss \\ rpt strip_tac
diff --git a/hol4/IntervalValidationScript.sml b/hol4/IntervalValidationScript.sml
index 8b42fa98d5771d1690958c0c137448149cf050c8..62c342676294460b3886a42adfa202c8520af23d 100644
--- a/hol4/IntervalValidationScript.sml
+++ b/hol4/IntervalValidationScript.sml
@@ -10,7 +10,7 @@ open simpLib realTheory realLib RealArith pred_setTheory sptreeTheory
open AbbrevsTheory ExpressionsTheory RealSimpsTheory FloverTactics
ExpressionAbbrevsTheory IntervalArithTheory CommandsTheory ssaPrgsTheory
MachineTypeTheory FloverMapTheory TypeValidatorTheory RealRangeArithTheory
- ExpressionSemanticsTheory;
+ ExpressionSemanticsTheory sqrtApproxTheory;
open preambleFloVer;
val _ = new_theory "IntervalValidation";
@@ -22,78 +22,90 @@ val _ = temp_overload_on("abs",``real$abs``);
val _ = temp_overload_on("max",``real$max``);
val _ = temp_overload_on("min",``real$min``);
+(** Define a global iteration count for square roots **)
+Definition ITERCOUNT_def:
+ (ITERCOUNT:num) = 4
+End
+
Definition validIntervalbounds_def:
validIntervalbounds f (absenv:analysisResult) (P:precond) (validVars:num_set) =
case (FloverMapTree_find f absenv) of
- | SOME (intv, _) =>
- (case f of
- | Var v =>
- if (lookup v validVars = SOME ())
- then T
- else (isSupersetInterval (P v) intv /\ IVlo (P v) <= IVhi (P v))
- | Const m n => isSupersetInterval (n,n) intv
- | Unop op f1 =>
- (case (FloverMapTree_find f1 absenv) of
- | SOME (iv, _) =>
- (if validIntervalbounds f1 absenv P validVars
- then
- case op of
- | Neg =>
- let new_iv = negateInterval iv in
- isSupersetInterval new_iv intv
- | Inv =>
- if IVhi iv < 0 \/ 0 < IVlo iv
- then
- let new_iv = invertInterval iv in
- isSupersetInterval new_iv intv
- else
- F
- else F)
- | _ => F)
- | Downcast m f1 =>
- (case (FloverMapTree_find f1 absenv) of
- | SOME (iv1, _) =>
- (if (validIntervalbounds f1 absenv P validVars) then
- ((isSupersetInterval intv iv1) /\ (isSupersetInterval iv1 intv))
- else F)
- | _ => F)
- | Binop op f1 f2 =>
- (if (validIntervalbounds f1 absenv P validVars /\
- validIntervalbounds f2 absenv P validVars)
- then
- (case (FloverMapTree_find f1 absenv, FloverMapTree_find f2 absenv) of
- | (SOME (iv1, _), SOME (iv2, _)) =>
- (case op of
- | Plus =>
- let new_iv = addInterval iv1 iv2 in
- isSupersetInterval new_iv intv
- | Sub =>
- let new_iv = subtractInterval iv1 iv2 in
- isSupersetInterval new_iv intv
- | Mult =>
- let new_iv = multInterval iv1 iv2 in
- isSupersetInterval new_iv intv
- | Div =>
- if (IVhi iv2 < 0 \/ 0 < IVlo iv2)
- then
- let new_iv = divideInterval iv1 iv2 in
- isSupersetInterval new_iv intv
- else F)
- | _ => F)
- else F)
- | Fma f1 f2 f3 =>
- (if (validIntervalbounds f1 absenv P validVars /\
- validIntervalbounds f2 absenv P validVars /\
- validIntervalbounds f3 absenv P validVars)
- then
- case FloverMapTree_find f1 absenv, FloverMapTree_find f2 absenv,
- FloverMapTree_find f3 absenv of
- | SOME (iv1, _), SOME (iv2, _), SOME (iv3, _) =>
- let new_iv = addInterval (multInterval iv1 iv2) iv3 in
+ | SOME (intv, _) =>
+ (case f of
+ | Var v =>
+ if (lookup v validVars = SOME ())
+ then T
+ else (isSupersetInterval (P v) intv /\ IVlo (P v) <= IVhi (P v))
+ | Const m n => isSupersetInterval (n,n) intv
+ | Unop op f1 =>
+ (case (FloverMapTree_find f1 absenv) of
+ | SOME (iv, _) =>
+ (if validIntervalbounds f1 absenv P validVars then
+ case op of
+ | Neg =>
+ let new_iv = negateInterval iv in
+ isSupersetInterval new_iv intv
+ | Inv =>
+ if IVhi iv < 0 \/ 0 < IVlo iv then
+ let new_iv = invertInterval iv in
+ isSupersetInterval new_iv intv
+ else F
+ (* Sqrt is transcendental -> We cannot directly compute it
+ thus we do a newton approximation of the lower and upper bounds *)
+ | Sqrt =>
+ if 0 < IVlo iv then
+ let sqrtLo = newton ITERCOUNT (IVlo iv * 0.99) (IVlo iv * 0.99);
+ sqrtHi = newton ITERCOUNT (IVhi iv * 1.01) (IVhi iv * 1.01);
+ new_iv = (sqrtLo, sqrtHi)
+ in
+ if (validate_newton_down sqrtLo (IVlo iv) ∧
+ validate_newton_up sqrtHi (IVhi iv)) then
isSupersetInterval new_iv intv
- | _, _, _ => F
- else F))
- | _ => F
+ else F
+ else F
+ else F)
+ | _ => F)
+ | Downcast m f1 =>
+ (case (FloverMapTree_find f1 absenv) of
+ | SOME (iv1, _) =>
+ (if (validIntervalbounds f1 absenv P validVars) then
+ ((isSupersetInterval intv iv1) /\ (isSupersetInterval iv1 intv))
+ else F)
+ | _ => F)
+ | Binop op f1 f2 =>
+ (if (validIntervalbounds f1 absenv P validVars /\
+ validIntervalbounds f2 absenv P validVars) then
+ (case (FloverMapTree_find f1 absenv, FloverMapTree_find f2 absenv) of
+ | (SOME (iv1, _), SOME (iv2, _)) =>
+ (case op of
+ | Plus =>
+ let new_iv = addInterval iv1 iv2 in
+ isSupersetInterval new_iv intv
+ | Sub =>
+ let new_iv = subtractInterval iv1 iv2 in
+ isSupersetInterval new_iv intv
+ | Mult =>
+ let new_iv = multInterval iv1 iv2 in
+ isSupersetInterval new_iv intv
+ | Div =>
+ if (IVhi iv2 < 0 \/ 0 < IVlo iv2) then
+ let new_iv = divideInterval iv1 iv2 in
+ isSupersetInterval new_iv intv
+ else F)
+ | _ => F)
+ else F)
+ | Fma f1 f2 f3 =>
+ (if (validIntervalbounds f1 absenv P validVars /\
+ validIntervalbounds f2 absenv P validVars /\
+ validIntervalbounds f3 absenv P validVars) then
+ case FloverMapTree_find f1 absenv, FloverMapTree_find f2 absenv,
+ FloverMapTree_find f3 absenv of
+ | SOME (iv1, _), SOME (iv2, _), SOME (iv3, _) =>
+ let new_iv = addInterval (multInterval iv1 iv2) iv3 in
+ isSupersetInterval new_iv intv
+ | _, _, _ => F
+ else F))
+ | _ => F
End
Theorem validIntervalbounds_eq =
@@ -127,8 +139,34 @@ val real_prove =
\\ rpt (qpat_x_assum `_ ==> ! x. _` kall_tac)
\\ REAL_ASM_ARITH_TAC;
+Theorem newton_method_pos:
+ ∀ x m.
+ 0 ≤ m ∧
+ 0 ≤ x ⇒ 0 ≤ newton n x m
+Proof
+ Induct_on ‘n’ >> gs[newton_def]
+ >> rpt strip_tac
+ >> first_x_assum irule >> rewrite_tac[real_div] >> gs[]
+ >> irule REAL_LE_MUL >> gs[]
+ >> irule REAL_LE_ADD >> gs[]
+ >> irule REAL_LE_MUL >> gs[]
+QED
+
+Theorem newton_method_lt_pos:
+ ∀ x m.
+ 0 < m ∧
+ 0 < x ⇒ 0 < newton n m x
+Proof
+ Induct_on ‘n’ >> gs[newton_def]
+ >> rpt strip_tac
+ >> first_x_assum irule >> rewrite_tac[real_div] >> gs[]
+ >> irule REAL_LT_ADD >> gs[]
+ >> irule REAL_LT_MUL >> gs[nonzerop_def]
+ >> COND_CASES_TAC >> gs[]
+QED
+
Theorem validIntervalbounds_sound:
- !(f:real expr) (A:analysisResult) (P:precond) (fVars:num_set)
+ ∀ (f:real expr) (A:analysisResult) (P:precond) (fVars:num_set)
(dVars:num_set) E Gamma.
validIntervalbounds f A P dVars /\ (* The checker succeeded *)
(* All defined vars have already been analyzed *)
@@ -176,7 +214,7 @@ Proof
(qspecl_then [`A`, `P`, `fVars`, `dVars`, `E`, `Gamma`] destruct)
\\ fs[] \\ rveq
\\ IMP_RES_TAC validRanges_single
- \\ simp[Once validRanges_def] \\ Cases_on `u`
+ \\ simp[Once validRanges_def] \\ Cases_on `u` \\ simp[]
(* Negation case *)
>- (qexists_tac `- vR` \\ fs[negateInterval_def, isSupersetInterval_def]
\\ rveq \\ Cases_on `iv` \\ fs[]
@@ -184,30 +222,72 @@ Proof
\\ irule Unop_neg' \\ qexistsl_tac [`REAL`, `REAL`, `vR`]
\\ fs[evalUnop_def, isCompat_def, toRTMap_def])
(* Inversion case *)
- \\ qexists_tac `1 / vR`
- \\ conj_tac
- >- (simp[toREval_def]
- \\ irule Unop_inv' \\ qexistsl_tac [`0`, `REAL`, `REAL`, `vR`]
- \\ simp[evalUnop_def, mTypeToR_def, perturb_def, REAL_INV_1OVER,
- mTypeToR_pos, isCompat_def, toRTMap_def]
- \\ CCONTR_TAC \\ fs[] \\ rveq
- \\ `0 < 0:real` by (real_prove)
- \\ fs[])
+ >- (
+ qexists_tac `1 / vR`
+ \\ conj_tac
+ >- (
+ simp[toREval_def]
+ \\ irule Unop_inv' \\ gs[] \\ qexistsl_tac [`0`, `REAL`, `vR`]
+ \\ simp[evalUnop_def, mTypeToR_def, perturb_def, REAL_INV_1OVER,
+ mTypeToR_pos, isCompat_def, toRTMap_def]
+ \\ CCONTR_TAC \\ fs[] \\ rveq
+ \\ `0 < 0:real` by (real_prove)
+ \\ fs[])
+ \\ conj_tac
+ \\ fs[invertInterval_def, isSupersetInterval_def] \\ rveq
+ \\ Cases_on `iv` \\ fs[]
+ >- (irule REAL_LE_TRANS \\ qexists_tac `1 / r` \\ conj_tac \\ fs[]
+ \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPL \\ fs[]
+ \\ real_prove)
+ >- (irule REAL_LE_TRANS \\ qexists_tac `1 / r` \\ conj_tac \\ fs[]
+ \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPR \\ fs[]
+ \\ real_prove)
+ >- (irule REAL_LE_TRANS \\ qexists_tac `1 / q` \\ conj_tac \\ fs[]
+ \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPL \\ fs[]
+ \\ real_prove)
+ \\ irule REAL_LE_TRANS \\ qexists_tac `1/q` \\ conj_tac \\ fs[]
+ \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPR \\ fs[]
+ \\ real_prove)
+ (* Sqrt case *)
+ \\ conj_tac >- gs[]
+ \\ qexists_tac `sqrt vR`
\\ conj_tac
- \\ fs[invertInterval_def, isSupersetInterval_def] \\ rveq
- \\ Cases_on `iv` \\ fs[]
- >- (irule REAL_LE_TRANS \\ qexists_tac `1 / r` \\ conj_tac \\ fs[]
- \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPL \\ fs[]
- \\ real_prove)
- >- (irule REAL_LE_TRANS \\ qexists_tac `1 / r` \\ conj_tac \\ fs[]
- \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPR \\ fs[]
+ >- (
+ simp[toREval_def]
+ \\ irule Unop_sqrt' \\ gs[] \\ qexistsl_tac [`0`, `REAL`, `vR`]
+ \\ simp[evalUnop_def, mTypeToR_def, perturb_def, REAL_INV_1OVER,
+ mTypeToR_pos, isCompat_def, toRTMap_def]
+ \\ CCONTR_TAC \\ fs[] \\ rveq
+ \\ `0 < 0:real` by (real_prove)
+ \\ fs[])
+ \\ gs[] \\ rveq
+ \\ qpat_x_assum `validate_newton_down _ _` $ mp_then Any mp_tac validate_newton_lb_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
\\ real_prove)
- >- (irule REAL_LE_TRANS \\ qexists_tac `1 / q` \\ conj_tac \\ fs[]
- \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPL \\ fs[]
+ \\ strip_tac
+ \\ ‘0 < SND iv ’
+ by (
+ irule REAL_LTE_TRANS \\ qexists_tac ‘FST iv’ \\ gs[]
+ \\ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[])
+ \\ qpat_x_assum `validate_newton_up _ _` $ mp_then Any mp_tac validate_newton_ub_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
\\ real_prove)
- \\ irule REAL_LE_TRANS \\ qexists_tac `1/q` \\ conj_tac \\ fs[]
- \\ fs[GSYM REAL_INV_1OVER] \\ irule REAL_INV_LE_ANTIMONO_IMPR \\ fs[]
- \\ real_prove)
+ \\ strip_tac
+ \\ conj_tac
+ \\ fs[isSupersetInterval_def] \\ rveq
+ >- (
+ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[]
+ \\ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[]
+ \\ irule transcTheory.SQRT_MONO_LE \\ real_prove)
+ \\ irule REAL_LE_TRANS \\ qexists_tac ‘sqrt (SND iv)’ \\ conj_tac
+ >- (
+ irule transcTheory.SQRT_MONO_LE \\ gs[]
+ \\ irule REAL_LE_TRANS \\ qexists_tac ‘FST iv’ \\ real_prove)
+ \\ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[])
(* Binary operator case *)
>- (rename1 `Binop b f1 f2` \\ simp[Once validRanges_def, toREval_def]
\\ rw_thm_asm `validTypes _ _` validTypes_def \\ fs[]
@@ -457,7 +537,29 @@ Proof
\\ irule REAL_LE_TRANS
\\ asm_exists_tac \\ fs[]
\\ irule REAL_LE_TRANS
- \\ find_exists_tac \\ fs[invertInterval_def]))
+ \\ find_exists_tac \\ fs[invertInterval_def])
+ \\ gs[valid_def, isSupersetInterval_def]
+ \\ irule REAL_LE_TRANS
+ \\ asm_exists_tac \\ gs[]
+ \\ irule REAL_LE_TRANS \\ once_rewrite_tac [CONJ_COMM]
+ \\ asm_exists_tac \\ gs[]
+ \\ qpat_x_assum `validate_newton_down _ _` $ mp_then Any mp_tac validate_newton_lb_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
+ \\ real_prove)
+ \\ strip_tac
+ \\ ‘0 < SND iv ’
+ by (irule REAL_LTE_TRANS \\ qexists_tac ‘FST iv’ \\ gs[])
+ \\ qpat_x_assum `validate_newton_up _ _` $ mp_then Any mp_tac validate_newton_ub_valid
+ \\ impl_tac
+ >- (conj_tac
+ >- (irule newton_method_pos \\ conj_tac \\ irule REAL_LE_MUL \\ gs[] \\ real_prove)
+ \\ real_prove)
+ \\ strip_tac
+ \\ irule REAL_LE_TRANS \\ asm_exists_tac \\ gs[]
+ \\ irule REAL_LE_TRANS \\ qexists_tac ‘sqrt(SND iv)’\\ gs[]
+ \\ irule transcTheory.SQRT_MONO_LE \\ real_prove)
>- (rename1 `Binop b f1 f2`
\\ rpt (first_x_assum (qspecl_then [`A`, `P`, `dVars`] destruct) \\ fs[])
\\ rveq \\ fs[]
diff --git a/hol4/RealIntervalInferenceScript.sml b/hol4/RealIntervalInferenceScript.sml
index aea0e78cb8e534fe8f687789559e89b69fe2c7fd..ae49afea261295f055dfd8ef6c3409fbb7814d30 100644
--- a/hol4/RealIntervalInferenceScript.sml
+++ b/hol4/RealIntervalInferenceScript.sml
@@ -6,7 +6,7 @@ open simpLib realTheory realLib RealArith;
open AbbrevsTheory ExpressionsTheory RealSimpsTheory FloverTactics
ExpressionAbbrevsTheory IntervalArithTheory CommandsTheory ssaPrgsTheory
MachineTypeTheory FloverMapTheory TypeValidatorTheory RealRangeArithTheory
- ExpressionSemanticsTheory;
+ ExpressionSemanticsTheory sqrtApproxTheory;
open preambleFloVer;
val _ = new_theory "RealIntervalInference";
@@ -15,6 +15,10 @@ val _ = monadsyntax.enable_monadsyntax();
val _ = List.app monadsyntax.enable_monad ["option"];
val _ = type_abbrev ("ivMap", ``:(real # real) fMap``);
+Definition ITERCOUNT_def:
+ ITERCOUNT = 4:num
+End
+
(** Unverified, trusted inferencer for interval bounds on expressions **)
Definition inferIntervalbounds_def:
inferIntervalbounds e (P:precond) (akk:ivMap) :ivMap option =
@@ -33,6 +37,12 @@ Definition inferIntervalbounds_def:
case op of
| Neg =>
SOME (FloverMapTree_insert e (negateInterval iv_e1) recRes)
+ | Sqrt =>
+ let sqrtLo = newton ITERCOUNT (IVlo iv_e1 * 0.99) (IVlo iv_e1 * 0.99);
+ sqrtHi = newton ITERCOUNT (IVhi iv_e1 * 1.01) (IVhi iv_e1 * 1.01);
+ new_iv = (sqrtLo, sqrtHi);
+ in
+ SOME (FloverMapTree_insert e new_iv recRes)
| _ => NONE
od
| Downcast m e1 =>
diff --git a/hol4/RealRangeArithScript.sml b/hol4/RealRangeArithScript.sml
index e365b97bf2ddaab6432fc5caee29291a43c34cd9..0c68021cd384ad965b86151630acb9f068f68915 100644
--- a/hol4/RealRangeArithScript.sml
+++ b/hol4/RealRangeArithScript.sml
@@ -24,7 +24,12 @@ val validRanges_def = Define `
((case e of
| Var x => T
| Const m v => T
- | Unop u e1 => validRanges e1 A E Gamma
+ | Unop u e1 =>
+ (u = Sqrt ==>
+ (! iv1 err.
+ FloverMapTree_find e1 A = SOME (iv1, err) ==>
+ (0 < IVlo iv1))) /\
+ validRanges e1 A E Gamma
| Binop b e1 e2 =>
(b = Div ==>
(! iv2 err.
diff --git a/hol4/semantics/ExpressionSemanticsScript.sml b/hol4/semantics/ExpressionSemanticsScript.sml
index c88abbdc8b4be4e2bc509e9519855fb9768321c9..d2db65dcfdaddd0d771e426ed0a0e26fee87eda5 100644
--- a/hol4/semantics/ExpressionSemanticsScript.sml
+++ b/hol4/semantics/ExpressionSemanticsScript.sml
@@ -43,13 +43,21 @@ Inductive eval_expr:
isCompat m mN /\
eval_expr E Gamma f1 v1 m ==>
eval_expr E Gamma (Unop Neg f1) (evalUnop Neg v1) mN) /\
- (!E Gamma m mN f1 v1 delta.
- Gamma (Unop Inv f1) = SOME mN /\
- isCompat m mN /\
+ (!E Gamma m m1 f1 v1 delta.
+ Gamma (Unop Inv f1) = SOME m /\
+ isCompat m1 m /\
abs delta <= (mTypeToR m) (evalUnop Inv v1) /\
- eval_expr E Gamma f1 v1 m /\
+ eval_expr E Gamma f1 v1 m1 /\
(v1 <> 0) ==>
- eval_expr E Gamma (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) mN) /\
+ eval_expr E Gamma (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) m) /\
+ (!E Gamma m m1 f1 v1 delta.
+ Gamma (Unop Sqrt f1) = SOME m /\
+ isCompat m1 m /\
+ abs delta <= (mTypeToR m) (evalUnop Sqrt v1) /\
+ eval_expr E Gamma f1 v1 m1 /\
+ (* leq would be fine here, but IEEE semantics are currently broken, thus we need lt *)
+ (0 < v1) ==>
+ eval_expr E Gamma (Unop Sqrt f1) (perturb (evalUnop Sqrt v1) m delta) m) /\
(!E Gamma m m1 f1 v1 delta.
Gamma (Downcast m f1) = SOME m /\
isMorePrecise m1 m /\
@@ -89,8 +97,8 @@ Theorem eval_expr_cases =
``eval_expr E Gamma (Downcast m1 e1) res m2``]
|> LIST_CONJ |> curry save_thm "eval_expr_cases";
-val [Var_load, Const_dist, Unop_neg, Unop_inv, Downcast_dist, Binop_dist, Fma_dist] =
- CONJ_LIST 7 eval_expr_rules;
+val [Var_load, Const_dist, Unop_neg, Unop_inv, Unop_sqrt, Downcast_dist, Binop_dist, Fma_dist] =
+ CONJ_LIST 8 eval_expr_rules;
Theorem Var_load = Var_load
Theorem Const_dist = Const_dist
@@ -127,17 +135,33 @@ Proof
QED
Theorem Unop_inv':
- ∀ m f1 v1 delta v m' mN E Gamma fBits.
+ ∀ m f1 v1 delta v mN m1 E Gamma fBits.
(abs delta) <= mTypeToR m (evalUnop Inv v1) /\
- eval_expr E Gamma f1 v1 m /\
+ eval_expr E Gamma f1 v1 m1 /\
(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'
+ isCompat m1 mN /\
+ m = mN ==>
+ eval_expr E Gamma (Unop Inv f1) v m
Proof
rpt strip_tac \\ rveq \\ irule Unop_inv \\ fs[]
+ \\ asm_exists_tac \\ gs[]
+QED
+
+Theorem Unop_sqrt':
+ ∀ m f1 v1 delta v m1 mN E Gamma fBits.
+ (abs delta) <= mTypeToR m (evalUnop Sqrt v1) /\
+ eval_expr E Gamma f1 v1 m1 /\
+ (0 < v1) /\
+ v = perturb (evalUnop Sqrt v1) mN delta /\
+ Gamma (Unop Sqrt f1) = SOME mN /\
+ isCompat m1 mN /\
+ m = mN ==>
+ eval_expr E Gamma (Unop Sqrt f1) v m
+Proof
+ rpt strip_tac \\ rveq \\ irule Unop_sqrt \\ fs[]
+ \\ asm_exists_tac \\ gs[]
QED
Theorem Downcast_dist':
@@ -307,7 +331,10 @@ Proof
>- (qexists_tac `delta` \\ fs[])
>- (rveq \\ qexistsl_tac [`m'`, `v1`] \\ fs[]
\\ first_x_assum drule \\ disch_then irule \\ fs[])
- >- (rveq \\ qexistsl_tac [`m'`, `v1`] \\ fs[]
+ >- (rveq \\ qexistsl_tac [`m1`, `v1`] \\ fs[]
+ \\ qexists_tac `delta` \\ fs[]
+ \\ first_x_assum drule \\ disch_then irule \\ fs[])
+ >- (rveq \\ qexistsl_tac [`m1`, `v1`] \\ fs[]
\\ qexists_tac `delta` \\ fs[]
\\ first_x_assum drule \\ disch_then irule \\ fs[])
>- (rveq \\ qexistsl_tac [`m1`, `m2`, `v1`, `v2`, `delta`]
diff --git a/hol4/semantics/ExpressionsScript.sml b/hol4/semantics/ExpressionsScript.sml
index 6bda980183f9395af350ab430b9a165af436bc4c..0283cdbedfe9e217bc23ed2c60140098a8a92a53 100644
--- a/hol4/semantics/ExpressionsScript.sml
+++ b/hol4/semantics/ExpressionsScript.sml
@@ -30,54 +30,60 @@ Expressions will use unary operators.
Define them first
**)
val _ = Datatype `
- unop = Neg | Inv`
+ unop = Neg | Inv | Sqrt`
(**
Define evaluation for unary operators on reals.
Errors are added in the exprression evaluation level later
**)
-val evalUnop_def = Define `
- evalUnop Neg v = - v /\
- evalUnop Inv v = inv(v:real)`
+Definition evalUnop_def:
+ evalUnop Neg v = - v ∧
+ evalUnop Inv v = inv(v:real) ∧
+ evalUnop Sqrt v = sqrt v
+End
(**
Define Expressions parametric over some value type 'v.
Will ease reasoning about different instantiations later.
**)
-val _ = Datatype `
+Datatype:
expr = Var num
| Const mType 'v
| Unop unop expr
| Binop binop expr expr
| Fma expr expr expr
- | Downcast mType expr`
+ | Downcast mType expr
+End
(**
Define evaluation of FMA (fused multiply-add) on reals
Errors are added on the exprression evaluation level later
**)
-val evalFma_def = Define `
- evalFma v1 v2 v3 = evalBinop Plus (evalBinop Mult v1 v2) v3`;
+Definition evalFma_def:
+ evalFma v1 v2 v3 = evalBinop Plus (evalBinop Mult v1 v2) v3
+End
-val toREval_def = Define `
+Definition toREval_def:
(toREval (Var n) = Var n) /\
(toREval (Const m n) = Const REAL n) /\
(toREval (Unop u e1) = Unop u (toREval e1)) /\
(toREval (Binop b e1 e2) = Binop b (toREval e1) (toREval e2)) /\
(toREval (Fma e1 e2 e3) = Fma (toREval e1) (toREval e2) (toREval e3)) /\
- (toREval (Downcast m e1) = toREval e1)`;
+ (toREval (Downcast m e1) = toREval e1)
+End
-val toRMap_def = Define `
+Definition toRMap_def:
toRMap (d:num -> mType option) (n:num) : mType option =
case d n of
| SOME m => SOME REAL
- | NONE => NONE`;
+ | NONE => NONE
+End
(**
Define the set of "used" variables of an exprression to be the set of variables
occuring in it
**)
-val usedVars_def = Define `
+Definition usedVars_def:
usedVars (e: 'a expr) :num_set =
case e of
| Var x => insert x () (LN)
@@ -85,7 +91,8 @@ val usedVars_def = Define `
| Binop b e1 e2 => union (usedVars e1) (usedVars e2)
| Fma e1 e2 e3 => union (usedVars e1) (union (usedVars e2) (usedVars e3))
| Downcast m e1 => usedVars e1
- | _ => LN`;
+ | _ => LN
+End
(* (** *)
(* Analogous lemma for unary exprressions *)
diff --git a/hol4/sqrtApproxScript.sml b/hol4/sqrtApproxScript.sml
new file mode 100644
index 0000000000000000000000000000000000000000..847d6360b6712ad38e55cb1deb0e25213a3d7c22
--- /dev/null
+++ b/hol4/sqrtApproxScript.sml
@@ -0,0 +1,61 @@
+open transcTheory realTheory realLib RealArith bossLib derivativeTheory;
+open preambleFloVer;
+
+val _ = new_theory "sqrtApprox";
+
+(*
+Definition newton_helper_def:
+ newton_helper 0 n (x:real) = x ∧
+ newton_helper (SUC n) (init:real) x = newton_helper n init ((x + (init / x)) / 2)
+End
+
+Definition newton_def:
+ newton (n:num) (x:real) = newton_helper n x x
+End
+
+Definition newton_def:
+ newton 0 n (x:real) = x ∧
+ newton (SUC n) m x = newton n m ((x + (&m / x)) / 2)
+End
+*)
+
+Definition newton_def:
+ newton 0 n (x:real) = x ∧
+ newton (SUC n) (m:real) x = newton n m ((x + (m / x)) / 2)
+End
+
+Definition validate_newton_down_def:
+ validate_newton_down estimate lb ⇔ (estimate pow 2) ≤ lb
+End
+
+Theorem validate_newton_lb_valid:
+ 0 ≤ estimate ∧
+ 0 ≤ lb ∧
+ validate_newton_down estimate lb ⇒
+ estimate ≤ sqrt lb
+Proof
+ rw[validate_newton_down_def]
+ >> qspecl_then [‘estimate pow 2’, ‘estimate’] mp_tac SQRT_POS_UNIQ
+ >> impl_tac >- gs[POW_POS]
+ >> disch_then $ once_rewrite_tac o single o GSYM
+ >> irule SQRT_MONO_LE >> gs[POW_POS]
+QED
+
+Definition validate_newton_up_def:
+ validate_newton_up estimate ub ⇔ ub ≤ (estimate pow 2)
+End
+
+Theorem validate_newton_ub_valid:
+ 0 ≤ estimate ∧
+ 0 ≤ ub ∧
+ validate_newton_up estimate ub ⇒
+ sqrt ub ≤ estimate
+Proof
+ rw[validate_newton_up_def]
+ >> qspecl_then [‘estimate pow 2’, ‘estimate’] mp_tac SQRT_POS_UNIQ
+ >> impl_tac >- gs[POW_POS]
+ >> disch_then $ once_rewrite_tac o single o GSYM
+ >> irule SQRT_MONO_LE >> gs[POW_POS]
+QED
+
+val _ = export_theory();