IEEE_connection.v 16.2 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

Definition valid_div a b (f:binary_float a b):=
  match f with
  |B754_finite _ _ _ _ _ _ => true
  |_ => false
  end.

Fixpoint eval_exp_float (e:exp (binary_float 53 1024)) (E:nat -> option (binary_float 53 1024)):=
  match e with
16
  |Const c => if (is_finite 53 1024 c) then Some c else None
17 18 19
  |Var _ x => E x
  |Binop b e1 e2 =>
   match eval_exp_float e1 E, eval_exp_float e2 E with
20 21 22
   |Some f1, Some f2 =>
    if (is_finite 53 1024 f1 && is_finite 53 1024 f2)
    then
23 24 25 26 27 28
      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 => if (valid_div f2) then Some (b64_div mode_NE f1 f2) else None
      end
29 30
    else
      None
31 32 33 34 35
     |_ , _ => None
   end
   | _ => None
  end.

36 37 38 39 40 41 42 43
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

44 45 46 47 48 49 50 51 52 53 54
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.

Definition toREnv E :env:= fun x:nat => match E x with
                               |None => None
                               |Some c => Some (B2R 53 1024 c)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
                                   end.

Definition Binop_to_Rop (b:binop) :=
  match b with
  |Plus => Rplus
  |Sub => Rminus
  |Mult => Rmult
  |Div => Rdiv
  end.

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)
72
                                 the_val))
73 74
     (bpow radix2 1024) = true).

75 76 77 78 79 80
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.

81 82 83
(* 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
84 85 86 87 88 89 90
  ((bpow radix2 (-1022)) <= (Rabs the_val))%R.

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.
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156

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.


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.

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.

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.

157 158 159 160
(* 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!
*)
161
Theorem eval_gives_sem e E v:
162 163
  (forall v,
     exists q, E v = Some q -> is_finite 53 1024 q = true) ->
164
  eval_exp_float e E = Some v ->
165 166
  no_overflow_in_eval ->
  no_underflow_in_eval ->
167 168 169 170
  is_finite 53 1024 v = true ->
  eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e) (B2R 53 1024 v).
Proof.
  revert v; induction e;
171
    intros v1 eval_succeeds eval_not_overf eval_not_underf finite_res;
172
    simpl in eval_succeeds.
173 174 175
  - econstructor.
    unfold toREnv.
    rewrite eval_succeeds; auto.
176 177
  - inversion eval_succeeds; subst.
    assert (perturb (B2R 53 1024 v1) 0 = B2R 53 1024 v1) as const_def.
178 179 180 181 182 183
    + unfold perturb.
      rewrite Rmult_plus_distr_l.
      rewrite Rmult_1_r.
      rewrite Rmult_0_r.
      rewrite Rplus_0_r.
      reflexivity.
184 185
    + simpl.
      setoid_rewrite <- const_def at 2.
186 187 188 189
      eapply Const_dist.
      rewrite Rabs_R0.
      eapply mEps_geq_zero.
  - inversion eval_succeeds.
190 191 192 193
  - 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.
194 195 196 197 198
    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 *.
199
    simpl.
200 201 202 203
    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)
204
      as relative_error_binop.
205 206 207 208 209 210 211
    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.
212
      destruct b.
213 214
      * pose proof (Bplus_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1)
          as addition_correct.
215
        assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true)
216
          as [finite_b0 finite_b1]
217
            by (apply b64_plus_preserves_finite; inversion eval_succeeds; subst; auto).
218 219 220 221 222 223 224
        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.
225 226 227 228 229
        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.
230
        unfold b64_plus.
231
        rewrite add_is_round.
232 233 234 235 236
        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) =
237 238
                        perturb (B2R 53 1024 b0 + B2R 53 1024 b1) delta)
               as delta_exists.
239
        { specialize (relative_error_binop (B2R 53 1024 b0 + B2R 53 1024 b1)%R).
240 241 242 243
          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.
244 245 246 247
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
248 249 250
            unfold Qnum, Qden. repeat rewrite <- Z2R_IZR.
            unfold Z2R. unfold P2R. simpl. lra. }
        { destruct delta_exists as [delta [delta_valid eval_perturb]].
251 252 253 254 255
          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]
256
            by (apply b64_minus_preserves_finite; inversion eval_succeeds; subst; auto).
257 258 259 260 261 262 263
        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.
264 265 266 267 268 269
        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.
270
        unfold b64_minus.
271
        rewrite minus_is_round.
272 273 274 275 276
        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) =
277 278
                        perturb (B2R 53 1024 b0 - B2R 53 1024 b1) delta)
               as delta_exists.
279
        { specialize (relative_error_binop (B2R 53 1024 b0 - B2R 53 1024 b1)%R).
280 281 282 283
          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.
284 285 286 287 288
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
289 290 291
            repeat rewrite <- Z2R_IZR.
            unfold Z2R; unfold P2R; simpl. lra. }
        { destruct delta_exists as [delta [delta_valid eval_perturb]].
292 293 294 295 296
          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]
297
            by (apply b64_mult_preserves_finite; inversion eval_succeeds; subst; auto).
298 299 300 301 302 303 304
        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.
305 306 307 308 309 310
        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.
311
        unfold b64_mult.
312
        rewrite mult_round.
313 314 315 316 317
        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) =
318 319
                        perturb (B2R 53 1024 b0 * B2R 53 1024 b1) delta)
               as delta_exists.
320
        { specialize (relative_error_binop (B2R 53 1024 b0 * B2R 53 1024 b1)%R).
321 322 323 324
          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.
325 326 327 328 329
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
330 331 332
            repeat rewrite <- Z2R_IZR.
            unfold Z2R; unfold P2R; simpl. lra. }
        { destruct delta_exists as [delta [delta_valid eval_perturb]].
333 334 335
          rewrite eval_perturb.
          econstructor; eauto.
          intros; congruence. }
336 337 338 339
      * assert (valid_div b1 = true)
          as div_valid
            by (destruct (valid_div b1); congruence).
        rewrite div_valid in eval_succeeds.
340
        assert (forall sgn, b1 = B754_infinity 53 1024 sgn -> False)
341
          by (intros; unfold valid_div in *; destruct b1; congruence).
342 343
        assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true)
          as [finite_b0 finite_b1]
344
            by (apply b64_div_preserves_finite; inversion eval_succeeds; subst; auto).
345 346 347 348 349 350 351
        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.
352 353 354 355
        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.
356
        assert (B2R 53 1024 b1 <> 0%R).
357 358 359
        { 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.
360
          pose proof (bpow_gt_0 radix2 e).
361 362
          rewrite Z2R_cond_Zopp in H0.
          unfold cond_Ropp in H0.
363 364
          destruct b.
          - pose proof (Zgt_pos_0 m).
365 366 367
            rewrite Z.gt_lt_iff in H2.
            apply Z2R_lt in H2.
            apply Rmult_integral in H0.
368
            unfold Z2R in *.
369
            destruct H0; try lra.
370
          - pose proof (Zgt_pos_0 m).
371 372 373 374 375 376
            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.
377
        unfold b64_div.
378
        rewrite round_correct.
379 380 381 382 383
        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) =
384 385
                        perturb (B2R 53 1024 b0 / B2R 53 1024 b1) delta)
               as delta_exists.
386
        { specialize (relative_error_binop (B2R 53 1024 b0 / B2R 53 1024 b1)%R).
387 388 389 390
          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.
391 392 393 394 395
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
396 397 398
            repeat rewrite <- Z2R_IZR.
            unfold Z2R; unfold P2R; simpl. lra. }
        { destruct delta_exists as [delta [delta_valid eval_perturb]].
399 400
          rewrite eval_perturb.
          econstructor; eauto. }
401
Qed.