Commit 8883bb13 authored by Heiko Becker's avatar Heiko Becker

Prove subtraction error bound

parent e7378eba
......@@ -147,3 +147,100 @@ let add_abs_err_bounded =
THEN split
THENL [
REWRITE_TAC[REAL_ABS_POS]; auto]]]]]);;
let sub_abs_err_bounded =
prove (
`!(e1:(real)exp) (e1R:real) (e1F:real) (e2:(real)exp) (e2R:real) (e2F:real)
(vR:real) (vF:real) (cenv:num->real) (err1:real) (err2:real).
eval_exp (&0) cenv e1 e1R /\
eval_exp machineEpsilon cenv e1 e1F /\
eval_exp (&0) cenv e2 e2R /\
eval_exp machineEpsilon cenv e2 e2F /\
eval_exp (&0) cenv (Binop Sub e1 e2) vR /\
eval_exp machineEpsilon (updEnv 2 e2F (updEnv 1 e1F cenv)) (Binop Sub (Var 1) (Var 2)) vF /\
abs (e1R - e1F) <= err1 /\
abs (e2R - e2F) <= err2 ==>
abs (vR - vF) <= err1 + err2 + (((abs e1F) + (abs e2F)) * machineEpsilon)`,
(* Prove that e1R and e2R are the correct values and that vR is e1R + e2R *)
intros "!e1 e1R e1F e2 e2R e2F vR vF cenv err1 err2; e1_real e1_float e2_real e2_float plus_real plus_float abs_e1 abs_e2"
THEN USE_THEN "plus_real"
(fun th -> LABEL_TAC "plus_real_inv"
(MP (SPECL [`&0:real`;`cenv:num->real`; `Sub:binop`;`e1:(real)exp`; `e2:(real)exp`; `vR:real`] binop_inv) th))
THEN destruct "plus_real_inv" "@delta. @v1. @v2. vR_eq eval_e1_v1 eval_e2_v2 abs_0"
THEN ASM_SIMP_TAC[]
THEN (USE_THEN "abs_0"
(fun th -> REWRITE_TAC [MP (SPECL [`vR:real`; `delta:real`] perturb_0_val) th]))
(* FIXME: UGLY REWRITES *)
THEN LABEL_TAC "eval_e1_0_det"
(SPECL [`e1:(real)exp`; `v1:real`; `e1R:real`; `cenv:num->real`] eval_0_det)
THEN SUBGOAL_TAC "eval_e1_conj" `eval_exp (&0) cenv e1 v1 /\ eval_exp (&0) cenv e1 e1R` [split THEN auto]
THEN mp_spec "eval_e1_0_det" "eval_e1_conj"
THEN LABEL_TAC "eval_e2_0_det" (SPECL [`e2:(real)exp`; `v2:real`; `e2R:real`; `cenv:num->real`] eval_0_det)
THEN SUBGOAL_TAC "eval_e2_conj" `eval_exp (&0) cenv e2 v2 /\ eval_exp (&0) cenv e2 e2R` [split THEN auto]
THEN mp_spec "eval_e2_0_det" "eval_e2_conj"
THEN ASM_SIMP_TAC[eval_binop]
(* Now unfold the float valued evaluation to get the deltas we need for the inequality *)
THEN USE_THEN "plus_float"
(fun th -> LABEL_TAC "plus_float_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`Sub:binop`;`(Var 1):(real)exp`; `(Var 2):(real)exp`; `vF:real`] binop_inv) th))
THEN destruct "plus_float_inv" "@delta2. @v1F. @v2F. vF_eq eval_e1_v1F eval_e2_v2F abs_mEps"
THEN USE_THEN "eval_e1_v1F"
(fun th -> LABEL_TAC "v1F_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`1:num`; `v1F:real`] var_inv) th))
THEN USE_THEN "eval_e2_v2F"
(fun th -> LABEL_TAC "v2F_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`2:num`; `v2F:real`] var_inv) th))
THEN ASM_REWRITE_TAC[updEnv; eval_binop; perturb]
(* TODO: Find out how to evaluate the conditional here! *)
THEN SUBGOAL_TAC "1_neq_2" `1:num = 2:num <=> F` [ARITH_TAC]
THEN ASM_REWRITE_TAC[]
THEN REWRITE_TAC [REAL_ADD_LDISTRIB; REAL_MUL_RID; REAL_NEG_ADD; real_sub]
THEN SUBGOAL_TAC "bounds_simplify"
`((e1R:real) + --(e2R:real)) + (--(e1F:real) + -- --(e2F:real)) + --(((e1F:real) + --(e2F:real)) * (delta2:real)) =
((e1R:real) + -- e1F) + ((-- e2R) + e2F) + -- ((e1F + -- e2F) * delta2)` [REAL_ARITH_TAC]
THEN REMOVE_THEN "bounds_simplify" (fun th -> REWRITE_TAC[th])
THEN mp REAL_LE_TRANS
THEN EXISTS_TAC `abs (((e1R:real) + --(e1F:real))) + abs ((--(e2R:real) + (e2F:real)) + --(((e1F:real) + --(e2F:real)) * (delta2:real)))`
THEN split
THENL[
REWRITE_TAC[REAL_ABS_TRIANGLE];
mp REAL_LE_ADD2
THEN split
THENL [
REWRITE_TAC[GSYM real_sub] THEN auto;
mp REAL_LE_TRANS
THEN EXISTS_TAC `abs ((-- (e2R:real) + (e2F:real))) + abs (--(((e1F:real) + -- (e2F:real)) * (delta2:real)))`
THEN split
THENL [
REWRITE_TAC [REAL_ABS_TRIANGLE];
mp REAL_LE_ADD2
THEN split
THENL [
ONCE_REWRITE_TAC[REAL_ADD_SYM] THEN ONCE_REWRITE_TAC[GSYM real_sub]
THEN ONCE_REWRITE_TAC[SPECL [`e2F:real`; `e2R:real`]REAL_ABS_SUB]
THEN auto;
ASM_REWRITE_TAC[REAL_ABS_NEG; REAL_ABS_MUL]
THEN mp REAL_LE_TRANS
THEN EXISTS_TAC `(abs (e1F:real) + abs (--(e2F:real))) * abs (delta2:real)`
THEN split
THENL [
mp REAL_LE_RMUL
THEN split
THENL [REWRITE_TAC[REAL_ABS_TRIANGLE]; REWRITE_TAC[REAL_ABS_POS]];
REWRITE_TAC [REAL_ABS_NEG]
THEN mp REAL_LE_LMUL
THEN split
THENL [
REAL_ARITH_TAC; auto]]]]]]);;
e (intros "!e1 e1R e1F e2 e2R e2F vR vF cenv err1 err2; e1_real e1_float e2_real e2_float plus_real plus_float abs_e1 abs_e2");;
e (USE_THEN "plus_real"
(fun th -> LABEL_TAC "plus_real_inv"
(MP (SPECL [`&0:real`;`cenv:num->real`; `Sub:binop`;`e1:(real)exp`; `e2:(real)exp`; `vR:real`] binop_inv) th)));;
e (destruct "plus_real_inv" "@delta. @v1. @v2. vR_eq eval_e1_v1 eval_e2_v2 abs_0");;
e (ASM_SIMP_TAC[]
THEN (USE_THEN "abs_0"
(fun th -> REWRITE_TAC [MP (SPECL [`vR:real`; `delta:real`] perturb_0_val) th])));;
(* FIXME: UGLY REWRITES *)
e (LABEL_TAC "eval_e1_0_det"
(SPECL [`e1:(real)exp`; `v1:real`; `e1R:real`; `cenv:num->real`] eval_0_det));;
e (SUBGOAL_TAC "eval_e1_conj" `eval_exp (&0) cenv e1 v1 /\ eval_exp (&0) cenv e1 e1R` [split THEN auto]);;
e (mp_spec "eval_e1_0_det" "eval_e1_conj");;
e (LABEL_TAC "eval_e2_0_det" (SPECL [`e2:(real)exp`; `v2:real`; `e2R:real`; `cenv:num->real`] eval_0_det));;
e (SUBGOAL_TAC "eval_e2_conj" `eval_exp (&0) cenv e2 v2 /\ eval_exp (&0) cenv e2 e2R` [split THEN auto]);;
e (mp_spec "eval_e2_0_det" "eval_e2_conj");;
e (ASM_SIMP_TAC[eval_binop]);;
e (clear ["eval_e2_0_det"; "eval_e2_conj"; "eval_e1_0_det"; "eval_e1_conj"; "eval_e1_v1"; "eval_e2_v2"]);;
(* Now unfold the float valued evaluation to get the deltas we need for the inequality *)
e (USE_THEN "plus_float"
(fun th -> LABEL_TAC "plus_float_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`Sub:binop`;`(Var 1):(real)exp`; `(Var 2):(real)exp`; `vF:real`] binop_inv) th)));;
e (destruct "plus_float_inv" "@delta2. @v1F. @v2F. vF_eq eval_e1_v1F eval_e2_v2F abs_mEps");;
e (USE_THEN "eval_e1_v1F"
(fun th -> LABEL_TAC "v1F_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`1:num`; `v1F:real`] var_inv) th)));;
e (USE_THEN "eval_e2_v2F"
(fun th -> LABEL_TAC "v2F_inv"
(MP (SPECL [`machineEpsilon:real`;
`(updEnv:num->real->(num->real)->num->real) 2 (e2F:real)
((updEnv:num->real->(num->real)->num->real) 1 (e1F:real)
(cenv:num->real))`;
`2:num`; `v2F:real`] var_inv) th)));;
e (ASM_REWRITE_TAC[updEnv; eval_binop; perturb]);;
(* TODO: Find out how to evaluate the conditional here! *)
e (SUBGOAL_TAC "1_neq_2" `1:num = 2:num <=> F` [ARITH_TAC]);;
e (ASM_REWRITE_TAC[]);;
(* We have now obtained all necessary values from the evaluations --> remove them for readability *)
e (clear ["1_neq_2"; "v2F_inv"; "v1F_inv"; "eval_e2_v2F"; "eval_e1_v1F"; "vF_eq"; "vR_eq"; "e1_real"; "e2_real"; "plus_real"; "plus_float"]);;
e (REWRITE_TAC [REAL_ADD_LDISTRIB; REAL_MUL_RID; real_sub]);;
e (SUBGOAL_TAC "bounds_simplify" `((e1R:real) + e2R) + -- ((e1F + e2F) + (e1F + e2F) * delta2) =
((e1R:real) + -- e1F) + (e2R + -- e2F) + -- ((e1F + e2F) * delta2)` [REAL_ARITH_TAC]);;
e (ASM_SIMP_TAC[]);;
e (mp REAL_LE_TRANS);;
e (EXISTS_TAC `abs (((e1R:real) + --(e1F:real))) + abs (((e2R:real) + --(e2F:real)) + --(((e1F:real) + (e2F:real)) * (delta2:real)))`);;
e (split);;
e (REWRITE_TAC[REAL_ABS_TRIANGLE]);;
e (mp REAL_LE_ADD2);;
e (split THENL [REWRITE_TAC[GSYM real_sub] THEN auto; ALL_TAC]);;
e (mp REAL_LE_TRANS);;
e (EXISTS_TAC `abs (((e2R:real) + --(e2F:real))) + abs (--(((e1F:real) + (e2F:real)) * (delta2:real)))`);;
e (split THENL [REWRITE_TAC [REAL_ABS_TRIANGLE]; ALL_TAC]);;
e (mp REAL_LE_ADD2);;
e (split THENL [REWRITE_TAC[GSYM real_sub] THEN auto; REWRITE_TAC[REAL_ABS_NEG; REAL_ABS_MUL]]);;
e (mp REAL_LE_LMUL);;
e (split THENL [REWRITE_TAC[REAL_ABS_POS]; auto]);;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment