Expressions.v 14.3 KB
Newer Older
1
(**
2
  Formalization of the base exprression language for the flover framework
3
 **)
4 5
From Flover
     Require Import Infra.RealRationalProps Infra.Ltacs.
6

7 8
From Flover
     Require Export Infra.Abbrevs Infra.NatSet Infra.MachineType.
9

10 11 12 13 14
(**
  Expressions will use binary operators.
  Define them first
**)
Inductive binop : Type := Plus | Sub | Mult | Div.
15

16
Definition binopEq (b1:binop) (b2:binop) :=
='s avatar
= committed
17 18 19 20 21 22
  match b1, b2 with
  | Plus, Plus => true
  | Sub,  Sub  => true
  | Mult, Mult => true
  | Div,  Div  => true
  | _,_ => false
23 24
  end.

25 26
(**
  Next define an evaluation function for binary operators on reals.
27
  Errors are added on the exprression evaluation level later.
28
 **)
29
Definition evalBinop (o:binop) (v1:R) (v2:R) :=
30 31 32 33 34 35
  match o with
  | Plus => Rplus v1 v2
  | Sub => Rminus v1 v2
  | Mult => Rmult v1 v2
  | Div => Rdiv v1 v2
  end.
36

37 38
Lemma binopEq_refl b:
  binopEq b b = true.
39 40 41 42
Proof.
  case b; auto.
Qed.

43 44
Lemma binopEq_compat_eq b1 b2:
  binopEq b1 b2 = true <-> b1 = b2.
45
Proof.
46 47 48 49 50 51 52 53 54 55
  split; case b1; case b2; intros; simpl in *; congruence.
Qed.

Lemma binopEq_compat_eq_false b1 b2:
  binopEq b1 b2 = false <-> ~ (b1 = b2).
Proof.
  split; intros neq.
  - hnf; intros; subst. rewrite binopEq_refl in neq.
    congruence.
  - destruct b1; destruct b2; cbv; congruence.
56 57
Qed.

58 59 60 61 62 63
(**
   Expressions will use unary operators.
   Define them first
 **)
Inductive unop: Type := Neg | Inv.

64
Definition unopEq (o1:unop) (o2:unop) :=
='s avatar
= committed
65 66 67 68
  match o1, o2 with
  | Neg, Neg => true
  | Inv, Inv => true
  | _ , _ => false
69 70
  end.

71 72
Lemma unopEq_refl b:
  unopEq b b = true.
73 74 75 76
Proof.
  case b; auto.
Qed.

77 78 79 80 81 82
Lemma unopEq_sym u1 u2:
  unopEq u1 u2 = unopEq u2 u1.
Proof.
  destruct u1,u2; compute; auto.
Qed.

83 84
Lemma unopEq_compat_eq b1 b2:
  unopEq b1 b2 = true <-> b1 = b2.
85
Proof.
86
  split; case b1; case b2; intros; simpl in *; congruence.
87 88
Qed.

89 90
(**
   Define evaluation for unary operators on reals.
91
   Errors are added in the exprression evaluation level later.
92
 **)
93
Definition evalUnop (o:unop) (v:R):=
94 95 96 97 98
  match o with
  |Neg => (- v)%R
  |Inv => (/ v)%R
  end .

99
Definition evalFma (v1:R) (v2:R) (v3:R):=
Nikita Zyuzin's avatar
Nikita Zyuzin committed
100
  evalBinop Plus v1 (evalBinop Mult v2 v3).
Nikita Zyuzin's avatar
Nikita Zyuzin committed
101

102
(**
103
  Define exprressions parametric over some value type V.
104
  Will ease reasoning about different instantiations later.
105
**)
106 107 108 109 110 111 112
Inductive expr (V:Type): Type :=
  Var: nat -> expr V
| Const: mType -> V -> expr V
| Unop: unop -> expr V -> expr V
| Binop: binop -> expr V -> expr V -> expr V
| Fma: expr V -> expr V -> expr V -> expr V
| Downcast: mType -> expr V -> expr V.
113

114
Fixpoint toRExp (e:expr Q) :=
115
  match e with
Nikita Zyuzin's avatar
Nikita Zyuzin committed
116 117 118 119
  | Var _ v => Var R v
  | Const m n => Const m (Q2R n)
  | Unop o e1 => Unop o (toRExp e1)
  | Binop o e1 e2 => Binop o (toRExp e1) (toRExp e2)
120
  | Fma e1 e2 e3 => Fma (toRExp e1) (toRExp e2) (toRExp e3)
Nikita Zyuzin's avatar
Nikita Zyuzin committed
121
  | Downcast m e1 => Downcast m (toRExp e1)
122
  end.
123

124
Fixpoint toREval (e:expr R) :=
125
  match e with
126
  | Var _ v => Var R v
127
  | Const _ n => Const REAL n
128 129
  | Unop o e1 => Unop o (toREval e1)
  | Binop o e1 e2 => Binop o (toREval e1) (toREval e2)
130
  | Fma e1 e2 e3 => Fma (toREval e1) (toREval e2) (toREval e3)
131
  | Downcast _ e1 =>   Downcast REAL (toREval e1)
132
  end.
133

134
Definition toRMap (d:nat -> option mType) (n:nat) :=
='s avatar
= committed
135
  match d n with
136
  | Some m => Some REAL
='s avatar
= committed
137
  | None => None
138
  end.
139

140 141
Arguments toRMap _ _/.

142
(**
143 144 145 146
  Finally, define an error function that computes an errorneous value for a
  given type. For a fixed-point datatype, truncation is used and any
  floating-point type is perturbed. As we need not compute on this function we
  define it in Prop.
147
**)
148 149 150
Definition perturb (rVal:R) (m:mType) (delta:R) :R :=
  match m with
  (* The Real-type has no error *)
151
  |REAL => rVal
152 153 154 155 156
  (* Fixed-point numbers have an absolute error *)
  |F w f => rVal + delta
  (* Floating-point numbers have a relative error *)
  | _ => rVal * (1 + delta)
  end.
Heiko Becker's avatar
Heiko Becker committed
157

158 159
Hint Unfold perturb.

160
(**
161
Define expression evaluation relation parametric by an "error" epsilon.
162
The result value exprresses float computations according to the IEEE standard,
163 164
using a perturbation of the real valued computation by (1 + delta), where
|delta| <= machine epsilon.
165
**)
166 167 168
Open Scope R_scope.
Inductive eval_expr (E:env) (Gamma: nat -> option mType)
  :(expr R) -> R -> mType -> Prop :=
169
| Var_load m x v:
170
    Gamma x = Some m ->
171
    E x = Some v ->
172
    eval_expr E Gamma (Var R x) v m
173
| Const_dist m n delta:
174 175
    Rabs delta <= mTypeToR m ->
    eval_expr E Gamma (Const m n) (perturb n m delta) m
176
| Unop_neg m f1 v1:
177 178
    eval_expr E Gamma f1 v1 m ->
    eval_expr E Gamma (Unop Neg f1) (evalUnop Neg v1) m
179
| Unop_inv m f1 v1 delta:
180
    Rabs delta <= mTypeToR m ->
181
    eval_expr E Gamma  f1 v1 m ->
182
    (~ v1 = 0)%R  ->
183
    eval_expr E Gamma (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) m
184 185
| Downcast_dist m m1 f1 v1 delta:
    isMorePrecise m1 m = true ->
186
    Rabs delta <= mTypeToR m ->
187
    eval_expr E Gamma f1 v1 m1 ->
188
    eval_expr E Gamma (Downcast m f1) (perturb v1 m delta) m
189
| Binop_dist m1 m2 op f1 f2 v1 v2 delta:
190
    Rabs delta <= mTypeToR (join m1 m2) ->
191 192
    eval_expr E Gamma f1 v1 m1 ->
    eval_expr E Gamma f2 v2 m2 ->
193
    ((op = Div) -> (~ v2 = 0)%R) ->
194
    eval_expr E Gamma (Binop op f1 f2) (perturb (evalBinop op v1 v2) (join m1 m2) delta) (join m1 m2)
195
| Fma_dist m1 m2 m3 f1 f2 f3 v1 v2 v3 delta:
196
    Rabs delta <= mTypeToR (join3 m1 m2 m3) ->
197 198 199
    eval_expr E Gamma f1 v1 m1 ->
    eval_expr E Gamma f2 v2 m2 ->
    eval_expr E Gamma f3 v3 m3 ->
200
    eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3).
201

202
Close Scope R_scope.
203

204
Hint Constructors eval_expr.
205

206 207 208
(** *)
(*   Show some simpler (more general) rule lemmata *)
(* **)
209
Lemma Const_dist' m n delta v m' E Gamma:
210 211
  Rle (Rabs delta) (mTypeToR m') ->
  v = perturb n m delta ->
212
  m' = m ->
213
  eval_expr E Gamma (Const m n) v m'.
214 215 216 217 218 219 220
Proof.
  intros; subst; auto.
Qed.

Hint Resolve Const_dist'.

Lemma Unop_neg' m f1 v1 v m' E Gamma:
221
  eval_expr E Gamma f1 v1 m ->
222 223
  v = evalUnop Neg v1 ->
  m' = m ->
224
  eval_expr E Gamma (Unop Neg f1) v m'.
225 226 227 228 229 230 231
Proof.
  intros; subst; auto.
Qed.

Hint Resolve Unop_neg'.

Lemma Unop_inv' m f1 v1 delta v m' E Gamma:
232
  Rle (Rabs delta) (mTypeToR m') ->
233
  eval_expr E Gamma  f1 v1 m ->
234
  (~ v1 = 0)%R  ->
235
  v = perturb (evalUnop Inv v1) m delta ->
236
  m' = m ->
237
  eval_expr E Gamma (Unop Inv f1) v m'.
238 239 240 241 242 243 244 245
Proof.
  intros; subst; auto.
Qed.

Hint Resolve Unop_inv'.

Lemma Downcast_dist' m m1 f1 v1 delta v m' E Gamma:
  isMorePrecise m1 m = true ->
246
  Rle (Rabs delta) (mTypeToR m') ->
247
  eval_expr E Gamma f1 v1 m1 ->
248
  v = (perturb v1 m delta) ->
249
  m' = m ->
250
  eval_expr E Gamma (Downcast m f1) v m'.
251 252 253 254 255 256 257
Proof.
  intros; subst; eauto.
Qed.

Hint Resolve Downcast_dist'.

Lemma Binop_dist' m1 m2 op f1 f2 v1 v2 delta v m' E Gamma:
258
  Rle (Rabs delta) (mTypeToR m') ->
259 260
  eval_expr E Gamma f1 v1 m1 ->
  eval_expr E Gamma f2 v2 m2 ->
261
  ((op = Div) -> (~ v2 = 0)%R) ->
262
  v = perturb (evalBinop op v1 v2) m' delta ->
263
  m' = join m1 m2 ->
264
  eval_expr E Gamma (Binop op f1 f2) v m'.
265 266 267 268 269
Proof.
  intros; subst; auto.
Qed.

Hint Resolve Binop_dist'.
270

271
Lemma Fma_dist' m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m' E Gamma:
272
  Rle (Rabs delta) (mTypeToR m') ->
273 274 275
  eval_expr E Gamma f1 v1 m1 ->
  eval_expr E Gamma f2 v2 m2 ->
  eval_expr E Gamma f3 v3 m3 ->
276
  v = perturb (evalFma v1 v2 v3) m' delta ->
Nikita Zyuzin's avatar
Nikita Zyuzin committed
277
  m' = join3 m1 m2 m3 ->
278
  eval_expr E Gamma (Fma f1 f2 f3) v m'.
Nikita Zyuzin's avatar
Nikita Zyuzin committed
279 280 281 282
Proof.
  intros; subst; auto.
Qed.

283
Hint Resolve Fma_dist'.
Nikita Zyuzin's avatar
Nikita Zyuzin committed
284

285
(**
286
  Define the set of "used" variables of an exprression to be the set of variables
287 288
  occuring in it
**)
289
Fixpoint usedVars (V:Type) (e:expr V) :NatSet.t :=
290 291
  match e with
  | Var _ x => NatSet.singleton x
292 293
  | Unop u e1 => usedVars e1
  | Binop b e1 e2 => NatSet.union (usedVars e1) (usedVars e2)
294
  | Fma e1 e2 e3 => NatSet.union (usedVars e1) (NatSet.union (usedVars e2) (usedVars e3))
295
  | Downcast _ e1 => usedVars e1
296 297
  | _ => NatSet.empty
  end.
298

299 300
Lemma toRMap_eval_REAL f v E Gamma m:
  eval_expr E (toRMap Gamma) (toREval f) v m -> m = REAL.
301 302 303 304 305 306 307 308 309 310 311 312
Proof.
  revert v E Gamma m.
  induction f; intros * eval_f; inversion eval_f; subst;
  repeat
    match goal with
    | H: context[toRMap _ _] |- _ => unfold toRMap in H
    | H: context[match ?Gamma ?v with | _ => _ end ] |- _ => destruct (Gamma v) eqn:?
    | H: Some ?m1 = Some ?m2 |- _ => inversion H; try auto
    | H: None = Some ?m |- _ => inversion H
    end; try auto.
  - eapply IHf; eauto.
  - eapply IHf; eauto.
313
  - assert (m1 = REAL)
314
      by (eapply IHf1; eauto).
315
    assert (m2 = REAL)
316 317
      by (eapply IHf2; eauto);
      subst; auto.
318
  - assert (m1 = REAL)
Nikita Zyuzin's avatar
Nikita Zyuzin committed
319
      by (eapply IHf1; eauto).
320
    assert (m2 = REAL)
Nikita Zyuzin's avatar
Nikita Zyuzin committed
321
      by (eapply IHf2; eauto).
322
    assert (m3 = REAL)
Nikita Zyuzin's avatar
Nikita Zyuzin committed
323 324
      by (eapply IHf3; eauto);
      subst; auto.
325 326
Qed.

327
(**
328
  If |delta| <= 0 then perturb v delta is exactly v.
329
**)
330
Lemma delta_0_deterministic (v:R) m (delta:R):
Heiko Becker's avatar
Heiko Becker committed
331
  (Rabs delta <= 0)%R ->
332
  perturb v m delta = v.
Heiko Becker's avatar
Heiko Becker committed
333 334
Proof.
  intros abs_0; apply Rabs_0_impl_eq in abs_0; subst.
335
  unfold perturb. destruct m; lra.
Heiko Becker's avatar
Heiko Becker committed
336 337
Qed.

338
(**
339
Evaluation with 0 as machine epsilon is deterministic
340
**)
341
Lemma meps_0_deterministic (f:expr R) (E:env) Gamma:
342
  forall v1 v2,
343 344
  eval_expr E (toRMap Gamma) (toREval f) v1 REAL ->
  eval_expr E (toRMap Gamma) (toREval f) v2 REAL ->
345 346
  v1 = v2.
Proof.
347 348 349 350 351 352
  induction f;
    intros v1 v2 ev1 ev2.
  - inversion ev1; inversion ev2; subst.
    rewrite H1 in H6.
    inversion H6; auto.
  - inversion ev1; inversion ev2; subst.
353
    simpl in *; subst; auto.
354 355
  - inversion ev1; inversion ev2; subst; try congruence.
    + rewrite (IHf v0 v3); eauto.
356
    + cbn in H1, H8. subst. rewrite (IHf v0 v3); eauto.
357
  - inversion ev1; inversion ev2; subst.
358 359 360 361
    assert (m0 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m3 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m1 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m2 = REAL) by (eapply toRMap_eval_REAL; eauto).
362
    subst.
363 364
    rewrite (IHf1 v0 v4); try auto.
    rewrite (IHf2 v3 v5); try auto.
Nikita Zyuzin's avatar
Nikita Zyuzin committed
365
  - inversion ev1; inversion ev2; subst.
366 367 368 369 370 371
    assert (m0 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m1 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m2 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m3 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m4 = REAL) by (eapply toRMap_eval_REAL; eauto).
    assert (m5 = REAL) by (eapply toRMap_eval_REAL; eauto).
Nikita Zyuzin's avatar
Nikita Zyuzin committed
372 373 374 375
    subst.
    rewrite (IHf1 v0 v5); try auto.
    rewrite (IHf2 v3 v6); try auto.
    rewrite (IHf3 v4 v7); try auto.
376
  - inversion ev1; inversion ev2; subst.
377 378
    apply REAL_least_precision in H1;
      apply REAL_least_precision in H7; subst.
379
    rewrite (IHf v0 v3); try auto.
380 381
Qed.

382
(**
Nikita Zyuzin's avatar
Nikita Zyuzin committed
383
Helping lemmas. Needed in soundness proof.
384
For each evaluation of using an arbitrary epsilon, we can replace it by
385
evaluating the subexprressions and then binding the result values to different
386
variables in the Environment.
387
 **)
388
Lemma binary_unfolding b f1 f2 E v1 v2 m1 m2 Gamma delta:
389
  (b = Div -> ~(v2 = 0 )%R) ->
390
  (Rabs delta <= (mTypeToR (join m1 m2)))%R ->
391 392
  eval_expr E Gamma f1 v1 m1 ->
  eval_expr E Gamma f2 v2 m2 ->
393
  eval_expr E Gamma (Binop b f1 f2) (perturb (evalBinop b v1 v2) (join m1 m2) delta) (join m1 m2) ->
394
  eval_expr (updEnv 2 v2 (updEnv 1 v1 emptyEnv))
395
           (updDefVars 2 m2 (updDefVars 1 m1 Gamma))
396
             (Binop b (Var R 1) (Var R 2)) (perturb (evalBinop b v1 v2) (join m1 m2) delta) (join m1 m2).
397
Proof.
398 399
  intros no_div_zero eval_f1 eval_f2 eval_float.
  econstructor; try auto.
400 401
Qed.

Nikita Zyuzin's avatar
Nikita Zyuzin committed
402
Lemma fma_unfolding f1 f2 f3 E v1 v2 v3 m1 m2 m3 Gamma delta:
403
  (Rabs delta <= (mTypeToR (join3 m1 m2 m3)))%R ->
404 405 406
  eval_expr E Gamma f1 v1 m1 ->
  eval_expr E Gamma f2 v2 m2 ->
  eval_expr E Gamma f3 v3 m3 ->
407
  eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3) ->
408
  eval_expr (updEnv 3 v3 (updEnv 2 v2 (updEnv 1 v1 emptyEnv)))
Nikita Zyuzin's avatar
Nikita Zyuzin committed
409
           (updDefVars 3 m3 (updDefVars 2 m2 (updDefVars 1 m1 Gamma)))
410
             (Fma (Var R 1) (Var R 2) (Var R 3)) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3).
Nikita Zyuzin's avatar
Nikita Zyuzin committed
411 412 413 414
Proof.
  econstructor; try auto.
Qed.

Heiko Becker's avatar
Heiko Becker committed
415 416 417
Lemma eval_eq_env e:
  forall E1 E2 Gamma v m,
    (forall x, E1 x = E2 x) ->
418 419
    eval_expr E1 Gamma e v m ->
    eval_expr E2 Gamma e v m.
Heiko Becker's avatar
Heiko Becker committed
420 421
Proof.
  induction e; intros;
422
    (match_pat (eval_expr _ _ _ _ _) (fun H => inversion H; subst; simpl in *));
Heiko Becker's avatar
Heiko Becker committed
423 424 425 426 427
    try eauto.
  eapply Var_load; auto.
  rewrite <- (H n); auto.
Qed.

428
Lemma eval_expr_ignore_bind e:
429
  forall x v m Gamma E,
430
    eval_expr E Gamma e v m ->
431 432
    ~ NatSet.In x (usedVars e) ->
    forall m_new v_new,
433
    eval_expr (updEnv x v_new E) (updDefVars x m_new Gamma) e v m.
434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455
Proof.
  induction e; intros * eval_e no_usedVar *; cbn in *;
    inversion eval_e; subst; try eauto.
  - assert (n <> x).
    { hnf. intros. subst. apply no_usedVar; set_tac. }
    rewrite <- Nat.eqb_neq in H.
    eapply Var_load.
    + unfold updDefVars.
      rewrite H; auto.
    + unfold updEnv.
      rewrite H; auto.
  - eapply Binop_dist'; eauto;
      [ eapply IHe1 | eapply IHe2];
      eauto;
      hnf; intros; eapply no_usedVar;
      set_tac.
  - eapply Fma_dist'; eauto;
      [eapply IHe1 | eapply IHe2 | eapply IHe3];
      eauto;
      hnf; intros; eapply no_usedVar;
        set_tac.
Qed.
456
(*
457 458 459
(**
Analogous lemma for unary expressions.
**)
460 461
Lemma unary_unfolding (e:expr R) (eps:R) (E:env) (v:R):
  (eval_expr eps E (Unop Inv e) v <->
462
   exists v1,
463 464
     eval_expr eps E e v1 /\
     eval_expr eps (updEnv 1 v1 E) (Unop Inv (Var R 1)) v).
465 466 467 468 469 470 471 472 473
Proof.
  split.
  - intros eval_un.
    inversion eval_un; subst.
    exists v1.
    repeat split; try auto.
    constructor; try auto.
    constructor; auto.
  - intros exists_val.
474 475
    destruct exists_val as [v1 [eval_f1 eval_e_E]].
    inversion eval_e_E; subst.
476 477 478
    inversion H1; subst.
    unfold updEnv in *; simpl in *.
    constructor; auto.
479
    inversion H3; subst; auto.
480
Qed. *)
481

482
(*   Using the parametric exprressions, define boolean exprressions for conditionals *)
483
(* **)
484 485 486
(* Inductive bexpr (V:Type) : Type := *)
(*   leq: expr V -> expr V -> bexpr V *)
(* | less: expr V -> expr V -> bexpr V. *)
487

488
(**
489
  Define evaluation of boolean exprressions
490
 **)
491 492 493 494
(* Inductive bval (E:env): (bexpr R) -> Prop -> Prop := *)
(*   leq_eval (f1:expr R) (f2:expr R) (v1:R) (v2:R): *)
(*     eval_expr E f1 v1 -> *)
(*     eval_expr E f2 v2 -> *)
495
(*     bval E (leq f1 f2) (Rle v1 v2) *)
496 497 498
(* |less_eval (f1:expr R) (f2:expr R) (v1:R) (v2:R): *)
(*     eval_expr E f1 v1 -> *)
(*     eval_expr E f2 v2 -> *)
499 500 501 502
(*     bval E (less f1 f2) (Rlt v1 v2). *)
(* (** *)
(*  Simplify arithmetic later by making > >= only abbreviations *)
(* **) *)
503 504
(* Definition gr := fun (V:Type) (f1: expr V) (f2: expr V) => less f2 f1. *)
(* Definition greq := fun (V:Type) (f1:expr V) (f2: expr V) => leq f2 f1. *)