Expressions.hl 13 KB
Newer Older
1
2
3
4
(*
  Formalization of the base expression language for the daisy framework
*)
needs "Infra/tactics.hl";;
5
needs "Infra/RealConstruction.hl";;
6
7
8
9
(*
  Expressions will use binary operators.
  Define them first
*)
10
let binop_IND, binop_REC = define_type
11
12
13
14
15
16
17
18
19
20
  "binop = Plus | Sub | Mult | Div ";;
(*
  Define an evaluation function for binary operators.
  Errors are added on the expression evaluation level later
*)
let eval_binop = new_recursive_definition binop_REC
  `(eval_binop Plus v1 v2 = v1 + v2) /\
  (eval_binop Sub v1 v2 = v1 - v2) /\
  (eval_binop Mult v1 v2 = v1 * v2) /\
  (eval_binop Div v1 v2 = v1 / v2)`;;
21
22
23
24
(*
  Define expressions parametric over some value type V.
  Will ease reasoning about different instantiations later.
*)
25
let exp_IND, exp_REC = define_type
26
  "exp = Var num
27
  |Param num
28
  | Const V
29
  | Binop binop exp exp";;
30
31
(*
  Define the machine epsilon for floating point operations.
32
  Current value: 2^(-53)
33
*)
34
let machineEpsilon =  define `machineEpsilon:real = realFromNum 1 1 53`;;
35
36
37
(*
  Define a perturbation function to ease writing of basic definitions
*)
38
let perturb = define `(perturb:real->real->real) (r:real) (e:real) = r * ((&1) + e)`;;
39
(*
40
41
  Define expression evaluation relation parametric by an "error" delta.
  This value will be used later to express float computations using a perturbation
42
43
  of the real valued computation by (1 + d)
*)
44
let eval_exp_RULES, eval_exp_IND, eval_exp_CASES = new_inductive_definition
45
46
  `(!eps env v.
      eval_exp eps env (Var v) (env v)) /\
47
48
49
  (!eps env v delta.
     abs delta <= eps ==>
     eval_exp eps env (Param v) (perturb (env v) delta)) /\
50
51
52
53
54
55
56
57
  (!eps env n delta.
     abs delta <= eps ==>
     eval_exp eps env (Const n) (perturb n delta)) /\
  (!eps env b e1 e2 v1 v2 delta.
     eval_exp eps env e1 v1 /\
     eval_exp eps env e2 v2 /\
     abs delta <= eps ==>
     eval_exp eps env (Binop b e1 e2) (perturb (eval_binop b v1 v2) delta))`;;
58

59
(* Inversion Lemmata, to do proofs similar to Coq proofs *)
60
let var_inv =
61
62
63
64
65
66
  prove(
    `!eps env v val.
       eval_exp eps env (Var v) val ==> val = env v`,
     ONCE_REWRITE_TAC
       [eval_exp_CASES] THEN
       INTRO_TAC "!eps env v val; cases_exp" THEN
67
68
69
70
71
72
73
       REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc")
       THENL [
         REMOVE_THEN "varc" (DESTRUCT_TAC "@v2. var_eq val_eq") THEN
           RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
         REMOVE_THEN "paramc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
         REMOVE_THEN "constc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN
           lcontra "var_eq" "exp";
74
75
        REMOVE_THEN "binopc"
          (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. var_eq _ _ _ _") THEN
76
          lcontra "var_eq" "exp"] THEN
77
78
79
80
81
82
83
84
85
86
87
88
89
90
       FIRST_ASSUM CONTR_TAC);;

let var_equiv =
  prove (
    `!eps env v val. eval_exp eps env (Var v) val <=> val = env v`,
    intros "!eps env v val"
      THEN EQ_TAC
      THENL[
        intros "eval_exp_var"
          THEN MATCH_MP_TAC var_inv
          THEN EXISTS_TAC `eps:real`
          THEN ASM_REWRITE_TAC[];
        intros "val_eq"
          THEN ASM_REWRITE_TAC[eval_exp_RULES]]);;
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
let param_inv =
  prove
 (`!eps env n v.
  eval_exp eps env (Param n) v ==> ?d. v = perturb (env n) d /\ abs d <= eps`,
  ONCE_REWRITE_TAC
    [eval_exp_CASES] THEN
  INTRO_TAC "!eps env n v; cases_exp" THEN
  REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc") THENL
   [REMOVE_THEN "varc" (DESTRUCT_TAC "@v. const_eq _") THEN
      lcontra "const_eq" "exp";
    REMOVE_THEN "paramc" (DESTRUCT_TAC "@n2. @d. c_eq v_eq abs_less") THEN
    EXISTS_TAC `d:real` THEN
    RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
    REMOVE_THEN "constc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
    REMOVE_THEN "binopc"
      (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. const_eq _ _ _ _") THEN
      lcontra "const_eq" "exp"]);;

110
111
112
113
114
115
116
let const_inv =
  prove
 (`!eps env n v.
  eval_exp eps env (Const n) v ==> ?d. v = perturb n d /\ abs d <= eps`,
  ONCE_REWRITE_TAC
    [eval_exp_CASES] THEN
  INTRO_TAC "!eps env n v; cases_exp" THEN
117
  REMOVE_THEN "cases_exp" (DESTRUCT_TAC "varc | paramc | constc | binopc") THENL
118
   [REMOVE_THEN "varc" (DESTRUCT_TAC "@v. const_eq _") THEN
119
120
   lcontra "const_eq" "exp";
    REMOVE_THEN "paramc" (DESTRUCT_TAC "@n. @d. var_eq _ _") THEN lcontra "var_eq" "exp";
121
122
123
124
125
    REMOVE_THEN "constc" (DESTRUCT_TAC "@n2. @d. c_eq v_eq abs_less") THEN
    EXISTS_TAC `d:real` THEN
    RULE_ASSUM_TAC (REWRITE_RULE [injectivity "exp"]) THEN ASM_REWRITE_TAC [];
    REMOVE_THEN "binopc"
      (DESTRUCT_TAC "@b. @e1. @e2. @v1. @v2. @d. const_eq _ _ _ _") THEN
126
   lcontra "const_eq" "exp"]);;
127

128
129
130
131
132
133
134
135
136
137
138
let binop_inv =
  prove
    ( `!eps env b e1 e2 v.
        eval_exp eps env (Binop b e1 e2) v ==>
        ?d v1 v2.
        v = perturb (eval_binop b v1 v2) d /\
            eval_exp eps env e1 v1 /\
            eval_exp eps env e2 v2 /\
            abs d <= eps`,
      INTRO_TAC "!eps env b e1 e2 v; cases_exp"
        THEN RULE_ASSUM_TAC (ONCE_REWRITE_RULE[eval_exp_CASES])
139
        THEN destruct "cases_exp" "varc | paramc | constc | binopc "
140
        THENL [
141
142
143
          destruct "varc" "@v. b_eq _" THEN lcontra "b_eq" "exp";
          destruct "paramc" "@n. @d. b_eq _ _" THEN lcontra "b_eq" "exp";
          destruct "constc" "@n. @d. b_eq _ _" THEN lcontra "b_eq" "exp";
144
145
146
          destruct "binopc" "@b2. @e1'. @e2'. @v1. @v2. @d. b_eq  v_eq eval_e1 eval_e2 abs_delta"
            THEN EXISTS_TAC `d:real` THEN EXISTS_TAC `v1:real` THEN EXISTS_TAC `v2:real`
            THEN RULE_ASSUM_TAC (REWRITE_RULE[injectivity "exp"]) THEN ASM_REWRITE_TAC[]]);;
147
148
149
150
151
152
153
154
155
156
157
158

let abs_leq_0_impl_zero =
  prove (
    `!d:real. abs d <= &0 ==> d = &0`,
    INTRO_TAC "!d; abs_leq_0"
      THEN SUBGOAL_TAC "abs_geq_0" `&0 <= abs d` [(REWRITE_TAC[REAL_ABS_POS])]
      THEN SUBGOAL_TAC "abs_geq_leq_0" `abs d <= &0 /\ &0 <= abs d` [CONJ_TAC THEN (ASM_REWRITE_TAC[])]
      THEN SUBGOAL_TAC "abs_eq_0" `abs d = &0` [ALL_TAC]
          THEN MP_TAC (ASSUME `abs d <= &0 /\ &0 <= abs d`)
          THEN ASM_REWRITE_TAC [SPECL [`abs d`; `&0:real`] REAL_LE_ANTISYM]
          THEN MP_TAC (ASSUME `abs d = &0`)
          THEN ASM_REWRITE_TAC [REAL_ABS_ZERO]);;
159

160

161
162
163
164
165
166
167
168
169
let perturb_0_val =
  prove (
    `!(v:real) (delta:real). abs delta <= &0 ==> perturb v delta = v`,
    intros "!v delta; abs_leq_0"
      THEN SIMP_TAC[perturb]
      THEN RULE_ASSUM_TAC (MATCH_MP abs_leq_0_impl_zero)
      THEN ASM_REWRITE_TAC[]
      THEN REAL_ARITH_TAC);;

170
171
g  (`!e:(real)exp v1:real v2:real env.
      eval_exp (&0) env e v1 /\ eval_exp (&0) env e v2 ==> v1 = v2`);;
172
e (MATCH_MP_TAC exp_IND THEN STRIP_TAC);;
173
(* Var Case *)
174
175
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
176
177
178
179
e (SUBGOAL_TAC "v1_eq_env_a" `v1:real = (env (a:num)):real`
    [MATCH_MP_TAC var_inv THEN EXISTS_TAC `&0` THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_env_a" `v2:real = (env (a:num)):real`
    [MATCH_MP_TAC var_inv THEN EXISTS_TAC `&0` THEN ASM_SIMP_TAC[]]);;
180
e (ASM_REWRITE_TAC[]);;
181
182
183
184
185
186
187
188
189
190
191
192
(* Param Case *)
e (STRIP_TAC);;
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
e (SUBGOAL_TAC "v1_eq_a" `?d. v1:real = perturb ((env (a:num)):real) d /\ abs d <= &0`
    [MATCH_MP_TAC param_inv THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_a" `?d. v2:real = perturb ((env (a:num)):real) d /\ abs d <= &0`
    [MATCH_MP_TAC param_inv THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_a" "@d1. v1_eq abs");;
e (destruct "v2_eq_a" "@d2. v2_eq abs2");;
e (USE_THEN "abs" (fun th -> ASM_REWRITE_TAC[MP (SPECL [`((env:num->real) a):real`; `d1:real`] perturb_0_val) th]));;
e (USE_THEN "abs2" (fun th -> ASM_REWRITE_TAC[MP (SPECL [`((env:num->real) a):real`; `d2:real`] perturb_0_val) th]));;
193
(* Const Case *)
194
e (STRIP_TAC);;
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
e (INTRO_TAC "!a v1 v2 env; eval_exps");;
e (destruct "eval_exps" "eval_v1 eval_v2");;
e (SUBGOAL_TAC "v1_eq_a" `?d. v1:real = perturb (a:real) d /\ abs d <= &0`
    [MATCH_MP_TAC const_inv THEN EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_a" `?d. v2:real = perturb (a:real) d /\ abs d <= &0`
    [MATCH_MP_TAC const_inv THEN EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_a" "@d1. v1_eq abs");;
e (destruct "v2_eq_a" "@d2. v2_eq abs2");;
e (SUBGOAL_TAC "abs_d1_0" `d1 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "abs_d2_0" `d2 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (ASM_REWRITE_TAC[]);;
(* Binop Case *)
e (REPEAT STRIP_TAC);;
e (FIND_ASSUM (LABEL_TAC "IH1") `!(v1:real) (v2:real) (env:num->real).
          eval_exp (&0) (env:num->real) (a1:(real)exp) (v1:real) /\
          eval_exp (&0) (env:num->real) (a1:(real)exp) (v2:real)
   ==> (v1:real) = (v2:real)`);;
e (FIND_ASSUM (LABEL_TAC "IH2") `!(v1:real) (v2:real) (env:num->real).
          eval_exp (&0) (env:num->real) (a2:(real)exp) (v1:real) /\
          eval_exp (&0) (env:num->real) (a2:(real)exp) (v2:real)
   ==> (v1:real) = (v2:real)`);;
e (SUBGOAL_TAC "v1_eq_op_a" `?d v1' v2'. v1:real = perturb (eval_binop (a0:binop) (v1':real) (v2':real)) d /\ eval_exp (&0) env a1 v1' /\ eval_exp (&0) env a2 v2' /\ abs d <= &0`
    [MATCH_MP_TAC binop_inv THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v2_eq_op_a" `?d v1' v2'. v2:real = perturb (eval_binop (a0:binop) (v1':real) (v2':real)) d /\ eval_exp (&0) env a1 v1' /\ eval_exp (&0) env a2 v2' /\ abs d <= &0`
    [MATCH_MP_TAC binop_inv THEN ASM_SIMP_TAC[]]);;
e (destruct "v1_eq_op_a" "@d1. @v11. @v21. v1_eq eval_a11 eval_a21 abs_0");;
e (destruct "v2_eq_op_a" "@d2. @v12. @v22. v2_eq eval_a12 eval_a22 abs_0");;
e (SUBGOAL_TAC "abs_d1_0" `d1 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "abs_d2_0" `d2 = &0` [MATCH_MP_TAC abs_leq_0_impl_zero THEN ASM_SIMP_TAC[]]);;
e (SUBGOAL_TAC "v11_eq_v12" `v11:real = v12:real` [REMOVE_THEN "IH1" (MATCH_MP_TAC)]);;
e (EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]);;
e (SUBGOAL_TAC "v21_eq_v22" `v21:real = v22:real` [REMOVE_THEN "IH2" (MATCH_MP_TAC)]);;
e (EXISTS_TAC `env:num->real` THEN ASM_SIMP_TAC[]);;
228
e (ASM_REWRITE_TAC[]);;
229
let eval_0_det = top_thm();;
230

231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
(* Alias for proof similarity *)
let Rsub_eq_Ropp_Rplus = real_sub;;

(** weird: cannot name lemma Rabs_err_simpl **)
let abs_err_simpl  =
  prove (
    `!(a:real) (b:real).
     abs (a - (a * (&1 + b))) = abs (a * b)`,
    intros "!a b"
      THEN REWRITE_TAC [REAL_ADD_LDISTRIB; REAL_MUL_RID]
      THEN SUBGOAL_TAC "arith_simp" `(a:real) - ((a:real) + (a:real) * (b:real)) = -- (a * b)` [REAL_ARITH_TAC]
      THEN REMOVE_THEN "arith_simp" (fun th -> REWRITE_TAC [th])
      THEN REWRITE_TAC [REAL_ABS_NEG]);;

(**
  TODO: Check wether we need Rabs (n * machineEpsilon) instead
**)

let const_abs_err_bounded =
  prove (
    `!(n:real) (nR:real) (nF:real) (cenv:num ->real).
      eval_exp (&0) cenv (Const n) nR /\
      eval_exp machineEpsilon cenv (Const n) nF ==>
      abs (nR - nF) <= abs n * machineEpsilon`,
    intros "!n nR nF cenv; eval_real eval_float"
      THEN RULE_ASSUM_TAC (MATCH_MP const_inv)
      THEN destruct "eval_real" "@d1. nR_eq abs_d1_0"
      THEN destruct "eval_float" "@d2. nF_eq abs_d2_leq"
      THEN ASM_SIMP_TAC []
      THEN USE_THEN "abs_d1_0"
      (fun th ->
         ASM_REWRITE_TAC [MP (SPECL [`(n:real)`; `d1:real`] perturb_0_val) th])
      THEN REWRITE_TAC [perturb;abs_err_simpl; REAL_ABS_MUL]
      THEN MATCH_MP_TAC REAL_LE_LMUL
      THEN ASM_REWRITE_TAC [REAL_ABS_POS]);;

(**
  TODO: Maybe improve bound by adding interval constraint and proving that it is leq than maxAbs of bounds
  (nlo <= cenv n <= nhi)%R ->   (Rabs (nR - nF) <= Rabs ((Rmax (Rabs nlo) (Rabs nhi)) * machineEpsilon))%R.
**)
let param_abs_err_bounded =
  prove (
    `!(n:num) (nR:real) (nF:real) (cenv:num->real).
      eval_exp (&0) cenv (Param n) nR /\
      eval_exp machineEpsilon cenv (Param n) nF ==>
      (abs (nR - nF) <= abs (cenv n) * machineEpsilon)`,
    intros "!n nR nF cenv; eval_real eval_float"
      THEN RULE_ASSUM_TAC (MATCH_MP param_inv)
      THEN destruct "eval_real" "@d1. nR_eq abs_d1_0"
      THEN destruct "eval_float" "@d2. nF_eq abs_d1_leq"
      THEN ASM_SIMP_TAC[]
      THEN USE_THEN "abs_d1_0"
      (fun th ->
         ASM_REWRITE_TAC [MP (SPECL [`(cenv n:real)`; `d1:real`] perturb_0_val) th])
      THEN REWRITE_TAC [perturb; abs_err_simpl; REAL_ABS_MUL]
      THEN MATCH_MP_TAC REAL_LE_LMUL THEN ASM_REWRITE_TAC[REAL_ABS_POS]);;
287
(*
Heiko Becker's avatar
Heiko Becker committed
288
  Using the parametric expressions, define boolean expressions for conditionals
289
*)
Heiko Becker's avatar
Heiko Becker committed
290
291
292
let bexp_INDUCT, bexp_REC = define_type
  "bexp = leq (V)exp (V)exp
        | less (V)exp (V)exp";;
293
294

(*
295
  Define evaluation of boolean expressions
296
*)
297
298
299
300
301
302
303
304
305
let bval_exp_RULES, bval_exp_IND, bval_exp_CASES = new_inductive_definition
  `(!eps env e1 e2 v1 v2.
      eval_exp eps env e1 v1 /\
      eval_exp eps env e2 v2 ==>
      bval eps env (leq e1 e2) (v1 <= v2))/\
  (!eps env e1 e2 v1 v2.
      eval_exp eps env e1 v1 /\
      eval_exp eps env e2 v2 ==>
      bval eps env (less e1 e2) (v1 < v2))`;;
306

Heiko Becker's avatar
Heiko Becker committed
307
308
309
310
311
(* Simplify arithmetic later by making > >= only abbreviations *)
let gr_simps  = define
  `gr = \(e1:(V)exp). \(e2:(V)exp). less e2 e1`;;
let greq_simps = define
  `greq = \(e1:(V)exp). \(e2:(V)exp). leq e2 e1`;;