IEEE_connection.v 18.3 KB
Newer Older
1
2
3
4
5
6
Require Import Coq.Reals.Reals Coq.QArith.QArith Coq.micromega.Psatz
        Coq.QArith.Qreals.
Require Import Daisy.Expressions Daisy.Infra.RationalSimps
        Daisy.Infra.RealRationalProps.
Require Import Flocq.Appli.Fappli_IEEE_bits Flocq.Appli.Fappli_IEEE
        Flocq.Core.Fcore_Raux Flocq.Prop.Fprop_relative.
7

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
Fixpoint eval_exp_float (e:exp (binary_float 53 1024)) (E:nat -> option (binary_float 53 1024)):=
  match e with
  | Var _ x => E x
  | Const m v => Some v
  | Unop Neg e =>
    match eval_exp_float e E with
    |Some v1 => Some (b64_opp v1)
    |_ => None
    end
  | Unop Inv e => None
  | Binop b e1 e2 =>
    match eval_exp_float e1 E, eval_exp_float e2 E with
    | Some f1, Some f2 =>
      match b with
      | Plus => Some (b64_plus mode_NE f1 f2)
      | Sub => Some (b64_minus mode_NE f1 f2)
      | Mult => Some (b64_mult mode_NE f1 f2)
      | Div => Some (b64_div mode_NE f1 f2)
      end
    |_ , _ => None
    end
  | _ => None
30
31
  end.

32
33
34
35
36
37
38
39
Definition optionLift (A B:Type) (e:option A) (some_cont:A -> B) (none_cont:B) :=
  match e with
  | Some v => some_cont v
  | None => none_cont
  end.

Definition normal_or_zero v :=
  Qeq_bool v 0 || Qle_bool (minValue M64) (Qabs v).
40

41
42
43
44
45
46
47
(* Lemma eval_exp_float_finite e E: *)
(*   forall v, eval_exp_float e E = Some v -> *)
(*        (forall n, exists v, E n = Some v -> Is_true (is_finite 53 1024 v)) -> *)
(*        Is_true (is_finite 53 1024 v). *)
(* Proof. *)
(*   induction e. *)
(*   - intros v eval_v; simpl in *)
48

49
50
51
52
53
54
55
(* Fixpoint toRExp e := *)
(*   match e with *)
(*   |Const c => Const (B2R 53 1024 c) *)
(*   |Var _ x => Var R x *)
(*   |Unop u e => Unop u (toRExp e) *)
(*   |Binop b e1 e2 => Binop b (toRExp e1) (toRExp e2) *)
(*   end. *)
56

57
58
59
60
(* Definition toREnv E :env:= fun x:nat => match E x with *)
(*                                |None => None *)
(*                                |Some c => Some (B2R 53 1024 c) *)
(*                                    end. *)
61

62
63
64
65
66
67
68
(* Definition Binop_to_Rop (b:binop) := *)
(*   match b with *)
(*   |Plus => Rplus *)
(*   |Sub => Rminus *)
(*   |Mult => Rmult *)
(*   |Div => Rdiv *)
(*   end. *)
69

70
71
72
73
74
75
76
77
78
(* Definition no_overflow b fl1 fl2 := *)
(*   let the_val := (Binop_to_Rop b) (B2R 53 1024 fl1) (B2R 53 1024 fl2) in *)
(*   (Rlt_bool *)
(*      (Rabs *)
(*         (Fcore_generic_fmt.round radix2 *)
(*                                  (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53) *)
(*                                  (round_mode mode_NE) *)
(*                                  the_val)) *)
(*      (bpow radix2 1024) = true). *)
79

80
81
82
83
84
(* Definition no_overflow_in_eval := *)
(*   forall b e1 e2 E fl1 fl2, *)
(*     eval_exp_float e1 E = Some fl1 -> *)
(*     eval_exp_float e2 E = Some fl2 -> *)
(*     no_overflow b fl1 fl2. *)
85

86
87
88
89
(* (* TODO: Maybe move to Prop *) *)
(* Definition no_underflow b fl1 fl2 := *)
(*   let the_val := (Binop_to_Rop b) (B2R 53 1024 fl1) (B2R 53 1024 fl2) in *)
(*   ((bpow radix2 (-1022)) <= (Rabs the_val))%R. *)
90

91
92
93
94
95
(* Definition no_underflow_in_eval := *)
(*   forall b e1 e2 E fl1 fl2, *)
(*     eval_exp_float e1 E = Some fl1 -> *)
(*     eval_exp_float e2 E = Some fl2 -> *)
(*     no_underflow b fl1 fl2. *)
96

97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
(* Lemma b64_plus_preserves_finite f1 f2: *)
(*   is_finite 53 1024 (b64_plus mode_NE f1 f2) = true -> *)
(*   is_finite 53 1024 f1 = true /\ is_finite 53 1024 f2 = true. *)
(* Proof. *)
(*   intros finite_add. *)
(*   case_eq (b64_plus mode_NE f1 f2); *)
(*     intros res add_res; *)
(*     unfold b64_plus, Bplus in *; *)
(*     case_eq f1; intros * f1_eq; *)
(*       try rewrite f1_eq in *; simpl in *; *)
(*         case_eq f2; intros * f2_eq; *)
(*           try rewrite f2_eq in *; simpl in *; *)
(*             try destruct (eqb b b0); *)
(*             try congruence; try auto. *)
(* Qed. *)
112
113


114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
(* Lemma b64_minus_preserves_finite f1 f2: *)
(*   is_finite 53 1024 (b64_minus mode_NE f1 f2) = true -> *)
(*   is_finite 53 1024 f1 = true /\ is_finite 53 1024 f2 = true. *)
(* Proof. *)
(*   intros finite_minus. *)
(*   case_eq (b64_minus mode_NE f1 f2); *)
(*     intros res add_res; *)
(*     unfold b64_minus, Bminus in *; *)
(*     case_eq f1; intros * f1_eq; *)
(*       try rewrite f1_eq in *; simpl in *; *)
(*         case_eq f2; intros * f2_eq; *)
(*           try rewrite f2_eq in *; simpl in *; *)
(*             try destruct (eqb b (negb b0)); *)
(*             try congruence; try auto. *)
(* Qed. *)
129

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
(* Lemma b64_mult_preserves_finite f1 f2: *)
(*   is_finite 53 1024 (b64_mult mode_NE f1 f2) = true -> *)
(*   is_finite 53 1024 f1 = true /\ is_finite 53 1024 f2 = true. *)
(* Proof. *)
(*   intros finite_mult. *)
(*   case_eq (b64_mult mode_NE f1 f2); *)
(*     intros res mult_res; *)
(*     unfold b64_mult, Bmult in *; *)
(*     case_eq f1; intros * f1_eq; *)
(*       try rewrite f1_eq in *; simpl in *; *)
(*         case_eq f2; intros * f2_eq; *)
(*           try rewrite f2_eq in *; simpl in *; *)
(*             try destruct (eqb b b0); *)
(*             try congruence; try auto. *)
(* Qed. *)
145

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
(* Lemma b64_div_preserves_finite f1 f2: *)
(*   is_finite 53 1024 (b64_div mode_NE f1 f2) = true -> *)
(*   (forall sgn, f2 = B754_infinity 53 1024 sgn -> False) -> *)
(*   is_finite 53 1024 f1 = true /\ is_finite 53 1024 f2 = true. *)
(* Proof. *)
(*   intros finite_div no_infty. *)
(*   case_eq (b64_div mode_NE f1 f2); *)
(*     intros res div_res; *)
(*     unfold b64_div, Bdiv in *; *)
(*     case_eq f1; intros * f1_eq; *)
(*       try rewrite f1_eq in *; simpl in *; *)
(*         case_eq f2; intros * f2_eq; *)
(*           try rewrite f2_eq in *; simpl in *; *)
(*             try congruence; try auto. *)
(* Qed. *)
161

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
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
228
229
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
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
(* (* Hypothesis Environment not finite --> remove finiteness assumption *)
(* no_overflow and no_underflow need to check only for expression in question. *)
(* As is they are ~ false, move the checks into eval function! *)
(* *) *)
(* Theorem eval_gives_sem e E v: *)
(*   (forall v, *)
(*      exists q, E v = Some q -> is_finite 53 1024 q = true) -> *)
(*   eval_exp_float e E = Some v -> *)
(*   no_overflow_in_eval -> *)
(*   no_underflow_in_eval -> *)
(*   is_finite 53 1024 v = true -> *)
(*   eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e) (B2R 53 1024 v). *)
(* Proof. *)
(*   revert v; induction e; *)
(*     intros v1 eval_succeeds eval_not_overf eval_not_underf finite_res; *)
(*     simpl in eval_succeeds. *)
(*   - econstructor. *)
(*     unfold toREnv. *)
(*     rewrite eval_succeeds; auto. *)
(*   - inversion eval_succeeds; subst. *)
(*     assert (perturb (B2R 53 1024 v1) 0 = B2R 53 1024 v1) as const_def. *)
(*     + unfold perturb. *)
(*       rewrite Rmult_plus_distr_l. *)
(*       rewrite Rmult_1_r. *)
(*       rewrite Rmult_0_r. *)
(*       rewrite Rplus_0_r. *)
(*       reflexivity. *)
(*     + simpl. *)
(*       setoid_rewrite <- const_def at 2. *)
(*       eapply Const_dist. *)
(*       rewrite Rabs_R0. *)
(*       eapply mEps_geq_zero. *)
(*   - inversion eval_succeeds. *)
(*   - case_eq (eval_exp_float e1 E); intros * eval_float_e1; *)
(*       case_eq (eval_exp_float e2 E); intros * eval_float_e2; *)
(*         rewrite eval_float_e1, eval_float_e2 in *; *)
(*         try congruence. *)
(*     specialize (IHe1 b0 (refl_equal (Some b0)) eval_not_overf eval_not_underf); *)
(*       specialize (IHe2 b1 (refl_equal (Some b1)) eval_not_overf eval_not_underf). *)
(*     specialize (eval_not_overf b e1 e2 E b0 b1 eval_float_e1 eval_float_e2). *)
(*     specialize (eval_not_underf b e1 e2 E b0 b1 eval_float_e1 eval_float_e2). *)
(*     unfold no_overflow, no_underflow in *. *)
(*     simpl. *)
(*     pose proof (fexp_correct 53 1024 eq_refl) as fexp_corr. *)
(*     pose proof (relative_error_N_ex radix2 *)
(*                                     (Fcore_FLT.FLT_exp (3 -1024 - 53) 53) *)
(*                                     (-1022)%Z 53%Z) *)
(*       as relative_error_binop. *)
(*     assert (forall k : Z, (-1022 < k)%Z -> (53 <= k - Fcore_FLT.FLT_exp (3 - 1024 - 53) 53 k)%Z) *)
(*       as exp_valid. *)
(*     + intros k k_pos. *)
(*       unfold Fcore_FLT.FLT_exp; simpl. *)
(*       destruct (Z.max_spec_le (k - 53) (-1074)); omega. *)
(*     + specialize (relative_error_binop exp_valid (fun x0 : Z => negb (Zeven x0))). *)
(*       clear exp_valid fexp_corr. *)
(*       destruct b. *)
(*       * pose proof (Bplus_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1) *)
(*           as addition_correct. *)
(*         assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true) *)
(*           as [finite_b0 finite_b1] *)
(*             by (apply b64_plus_preserves_finite; inversion eval_succeeds; subst; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e1) (B2R 53 1024 b0)) *)
(*           as eval_e1 *)
(*             by (apply IHe1; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e2) (B2R 53 1024 b1)) *)
(*           as eval_e2 *)
(*             by (apply IHe2; auto). *)
(*         clear IHe1 IHe2. *)
(*         unfold Binop_to_Rop in eval_not_overf. *)
(*         specialize (addition_correct finite_b0 finite_b1). *)
(*         rewrite eval_not_overf in addition_correct. *)
(*         destruct addition_correct as [add_is_round _]. *)
(*         inversion eval_succeeds; subst. *)
(*         unfold b64_plus. *)
(*         rewrite add_is_round. *)
(*         assert (exists delta,(Rabs(delta) <= Q2R machineEpsilon)%R /\ *)
(*                         Fcore_generic_fmt.round radix2 *)
(*                                                 (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53) *)
(*                                                 (round_mode mode_NE) *)
(*                                                 (B2R 53 1024 b0 + B2R 53 1024 b1) = *)
(*                         perturb (B2R 53 1024 b0 + B2R 53 1024 b1) delta) *)
(*                as delta_exists. *)
(*         { specialize (relative_error_binop (B2R 53 1024 b0 + B2R 53 1024 b1)%R). *)
(*           destruct relative_error_binop as [eps [abs_less round_eps]]. *)
(*           - auto. (* This is the proof that we do not handle denormals! *) *)
(*           - unfold round_mode; rewrite round_eps. *)
(*             unfold perturb; exists eps. *)
(*             split; try auto. *)
(*             eapply Rle_trans; eauto. *)
(*             unfold machineEpsilon. *)
(*             unfold Q2R. *)
(*             unfold Qnum, Qden. repeat rewrite <- Z2R_IZR. *)
(*             unfold Z2R. unfold P2R. simpl. lra. } *)
(*         { destruct delta_exists as [delta [delta_valid eval_perturb]]. *)
(*           rewrite eval_perturb. *)
(*           econstructor; eauto. *)
(*           intros; congruence. } *)
(*       * assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true) *)
(*           as [finite_b0 finite_b1] *)
(*             by (apply b64_minus_preserves_finite; inversion eval_succeeds; subst; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e1) (B2R 53 1024 b0)) *)
(*           as eval_e1 *)
(*             by (apply IHe1; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e2) (B2R 53 1024 b1)) *)
(*           as eval_e2 *)
(*             by (apply IHe2; auto). *)
(*         clear IHe1 IHe2. *)
(*         pose proof (Bminus_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1 finite_b0 finite_b1) *)
(*           as minus_correct. *)
(*         unfold Binop_to_Rop in eval_not_overf. *)
(*         rewrite eval_not_overf in minus_correct. *)
(*         destruct minus_correct as [minus_is_round _]. *)
(*         inversion eval_succeeds; subst. *)
(*         unfold b64_minus. *)
(*         rewrite minus_is_round. *)
(*         assert (exists delta,(Rabs(delta) <= Q2R machineEpsilon)%R /\ *)
(*                         Fcore_generic_fmt.round radix2 *)
(*                                                 (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53) *)
(*                                                 (round_mode mode_NE) *)
(*                                                 (B2R 53 1024 b0 - B2R 53 1024 b1) = *)
(*                         perturb (B2R 53 1024 b0 - B2R 53 1024 b1) delta) *)
(*                as delta_exists. *)
(*         { specialize (relative_error_binop (B2R 53 1024 b0 - B2R 53 1024 b1)%R). *)
(*           destruct relative_error_binop as [eps [abs_less round_eps]]. *)
(*           - auto. (* This is the proof that we do not handle denormals! *) *)
(*           - unfold round_mode; rewrite round_eps. *)
(*             unfold perturb; exists eps. *)
(*             split; try auto. *)
(*             eapply Rle_trans; eauto. *)
(*             unfold machineEpsilon. *)
(*             unfold Q2R. *)
(*             unfold Qnum, Qden. *)
(*             repeat rewrite <- Z2R_IZR. *)
(*             unfold Z2R; unfold P2R; simpl. lra. } *)
(*         { destruct delta_exists as [delta [delta_valid eval_perturb]]. *)
(*           rewrite eval_perturb. *)
(*           econstructor; eauto. *)
(*           intros; congruence. } *)
(*       * assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true) *)
(*           as [finite_b0 finite_b1] *)
(*             by (apply b64_mult_preserves_finite; inversion eval_succeeds; subst; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e1) (B2R 53 1024 b0)) *)
(*           as eval_e1 *)
(*             by (apply IHe1; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e2) (B2R 53 1024 b1)) *)
(*           as eval_e2 *)
(*             by (apply IHe2; auto). *)
(*         clear IHe1 IHe2. *)
(*         pose proof (Bmult_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1) *)
(*           as mult_correct. *)
(*         unfold Binop_to_Rop in eval_not_overf. *)
(*         rewrite eval_not_overf in mult_correct. *)
(*         destruct mult_correct as [mult_round _]. *)
(*         inversion eval_succeeds; subst. *)
(*         unfold b64_mult. *)
(*         rewrite mult_round. *)
(*         assert (exists delta,(Rabs(delta) <= Q2R machineEpsilon)%R /\ *)
(*                         Fcore_generic_fmt.round radix2 *)
(*                                                 (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53) *)
(*                                                 (round_mode mode_NE) *)
(*                                                 (B2R 53 1024 b0 * B2R 53 1024 b1) = *)
(*                         perturb (B2R 53 1024 b0 * B2R 53 1024 b1) delta) *)
(*                as delta_exists. *)
(*         { specialize (relative_error_binop (B2R 53 1024 b0 * B2R 53 1024 b1)%R). *)
(*           destruct relative_error_binop  as [eps [abs_less round_eps]]. *)
(*           - auto. (* This is the proof that we do not handle denormals! *) *)
(*           - unfold round_mode; rewrite round_eps. *)
(*             unfold perturb; exists eps. *)
(*             split; try auto. *)
(*             eapply Rle_trans; eauto. *)
(*             unfold machineEpsilon. *)
(*             unfold Q2R. *)
(*             unfold Qnum, Qden. *)
(*             repeat rewrite <- Z2R_IZR. *)
(*             unfold Z2R; unfold P2R; simpl. lra. } *)
(*         { destruct delta_exists as [delta [delta_valid eval_perturb]]. *)
(*           rewrite eval_perturb. *)
(*           econstructor; eauto. *)
(*           intros; congruence. } *)
(*       * assert (valid_div b1 = true) *)
(*           as div_valid *)
(*             by (destruct (valid_div b1); congruence). *)
(*         rewrite div_valid in eval_succeeds. *)
(*         assert (forall sgn, b1 = B754_infinity 53 1024 sgn -> False) *)
(*           by (intros; unfold valid_div in *; destruct b1; congruence). *)
(*         assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true) *)
(*           as [finite_b0 finite_b1] *)
(*             by (apply b64_div_preserves_finite; inversion eval_succeeds; subst; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e1) (B2R 53 1024 b0)) *)
(*           as eval_e1 *)
(*             by (apply IHe1; auto). *)
(*         assert (eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e2) (B2R 53 1024 b1)) *)
(*           as eval_e2 *)
(*             by (apply IHe2; auto). *)
(*         clear IHe1 IHe2. *)
(*         pose proof (Bdiv_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1) *)
(*           as div_correct. *)
(*         unfold Binop_to_Rop in eval_not_overf. *)
(*         rewrite eval_not_overf in div_correct. *)
(*         assert (B2R 53 1024 b1 <> 0%R). *)
(*         { hnf. intros. unfold B2R in H0. unfold valid_div in div_valid. destruct b1; simpl in H; try congruence. *)
(*           unfold Fcore_defs.F2R in H0. *)
(*           simpl in H0. *)
(*           pose proof (bpow_gt_0 radix2 e). *)
(*           rewrite Z2R_cond_Zopp in H0. *)
(*           unfold cond_Ropp in H0. *)
(*           destruct b. *)
(*           - pose proof (Zgt_pos_0 m). *)
(*             rewrite Z.gt_lt_iff in H2. *)
(*             apply Z2R_lt in H2. *)
(*             apply Rmult_integral in H0. *)
(*             unfold Z2R in *. *)
(*             destruct H0; try lra. *)
(*           - pose proof (Zgt_pos_0 m). *)
(*             rewrite Z.gt_lt_iff in H2. *)
(*             apply Z2R_lt in H2. *)
(*             apply Rmult_integral in H0. *)
(*             unfold Z2R in *; destruct H0; lra. } *)
(*         destruct div_correct as [round_correct _]; try auto. *)
(*         inversion eval_succeeds; subst. *)
(*         unfold b64_div. *)
(*         rewrite round_correct. *)
(*         assert (exists delta,(Rabs(delta) <= Q2R machineEpsilon)%R /\ *)
(*                         Fcore_generic_fmt.round radix2 *)
(*                                                 (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53) *)
(*                                                 (round_mode mode_NE) *)
(*                                                 (B2R 53 1024 b0 / B2R 53 1024 b1) = *)
(*                         perturb (B2R 53 1024 b0 / B2R 53 1024 b1) delta) *)
(*                as delta_exists. *)
(*         { specialize (relative_error_binop (B2R 53 1024 b0 / B2R 53 1024 b1)%R). *)
(*           destruct relative_error_binop  as [eps [abs_less round_eps]]. *)
(*           - auto. (* This is the proof that we do not handle denormals! *) *)
(*           - unfold round_mode; rewrite round_eps. *)
(*             unfold perturb; exists eps. *)
(*             split; try auto. *)
(*             eapply Rle_trans; eauto. *)
(*             unfold machineEpsilon. *)
(*             unfold Q2R. *)
(*             unfold Qnum, Qden. *)
(*             repeat rewrite <- Z2R_IZR. *)
(*             unfold Z2R; unfold P2R; simpl. lra. } *)
(*         { destruct delta_exists as [delta [delta_valid eval_perturb]]. *)
(*           rewrite eval_perturb. *)
(*           econstructor; eauto. } *)
(* Qed. *)