IEEE_connection.v 18.4 KB
Newer Older
1
Require Import Coq.Reals.Reals Coq.QArith.QArith Coq.QArith.Qabs Coq.micromega.Psatz
2 3 4 5 6
        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. *)