IEEE_connection.v 16.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
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.Prop.Fprop_relative Flocq.Core.Fcore_defs Flocq.Calc.Fcalc_ops Flocq.Calc.Fcalc_div. *)
Require Import Flocq.Appli.Fappli_IEEE_bits Flocq.Appli.Fappli_IEEE Flocq.Core.Fcore_Raux Flocq.Prop.Fprop_relative.

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
  |Const c => Some c
  |Var _ x => E x
  |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 => if (valid_div f2) then Some (b64_div mode_NE f1 f2) else None
      end
     |_ , _ => None
   end
   | _ => None
  end.

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)
                               end.

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.

Theorem eval_gives_sem e E v:
  eval_exp_float e E = Some v ->
113 114
  (forall e E v,
      eval_exp_float e E = Some v ->
115 116 117 118 119 120
      (Rlt_bool
        (Rabs
           (Fcore_generic_fmt.round radix2
                                    (Fcore_FLT.FLT_exp (3 - 1024 - 53) 53)
                                    (round_mode mode_NE)
                                    (B2R 53 1024 v)))
121 122
        (bpow radix2 1024) = false)
      -> False) ->
123 124 125 126
  is_finite 53 1024 v = true ->
  eval_exp (Q2R machineEpsilon) (toREnv E) (toRExp e) (B2R 53 1024 v).
Proof.
  revert v; induction e;
127 128
    intros v1 eval_succeeds no_overflow finite_res;
    simpl in eval_succeeds.
129 130 131
  - econstructor.
    unfold toREnv.
    rewrite eval_succeeds; auto.
132 133
  - inversion eval_succeeds; subst.
    assert (perturb (B2R 53 1024 v1) 0 = B2R 53 1024 v1) as const_def.
134 135 136 137 138 139
    + unfold perturb.
      rewrite Rmult_plus_distr_l.
      rewrite Rmult_1_r.
      rewrite Rmult_0_r.
      rewrite Rplus_0_r.
      reflexivity.
140 141
    + simpl.
      setoid_rewrite <- const_def at 2.
142 143 144 145
      eapply Const_dist.
      rewrite Rabs_R0.
      eapply mEps_geq_zero.
  - inversion eval_succeeds.
146 147 148 149 150 151
  - 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)) no_overflow);
      specialize (IHe2 b1 (refl_equal (Some b1)) no_overflow).
152 153
    specialize (no_overflow (Binop b e1 e2) E).
    simpl.
154 155 156 157
    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)
158
      as relative_error_binop.
159 160 161 162 163 164 165
    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.
166 167 168 169 170 171 172
      destruct b.
      * pose proof (Bplus_correct 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1).
        unfold b64_plus in *.
        unfold Bplus in *.
        finite_b0 finite_b1).

        assert (is_finite 53 1024 b0 = true /\ is_finite 53 1024 b1 = true)
173 174 175 176 177 178 179 180 181
          as [finite_b0 finite_b1]
            by (apply b64_plus_preserves_finite; 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.
182 183 184 185 186 187 188 189 190 191
        simpl in *.
        remember (179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474124377767893424865485276302219601246094119453082952085005768838150682342462881473913110540827237163350510684586298239947245938479716304835356329624224137216)%R as cst.
        unfold b64_plus in *.
        destruct (
            Rlt_bool
         (Rabs
            (Fcore_generic_fmt.round radix2 (Fcore_FLT.FLT_exp (-1074) 53)
               (Fcore_generic_fmt.Znearest (fun x : Z => negb (Zeven x))) (B2R 53 1024 b0 + B2R 53 1024 b1))) cst).
        Focus 2.
        exfalso; apply no_overflow; auto. ; try congruence; auto).
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 407 408 409
        rewrite H0 in *.
        destruct H.
        destruct H1.
        unfold b64_plus.
        rewrite H.
        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).
        { specialize (relative_error_binop (B2R 53 1024 b0 + B2R 53 1024 b1)%R).
          destruct relative_error_binop.
          - unfold bpow, radix2.
            (* This is the proof that we do not handle denormals! *) admit.
          - destruct H3.
            unfold round_mode; simpl in *.
            rewrite H4.
            unfold perturb.
            exists x.
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
            assert (2 ^ 53 = 9007199254740992)%positive by (cbv; auto).
            rewrite H5.
            lra. }
        { destruct H3 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; auto).
        specialize (IHe1 b0).
        specialize (IHe2 b1).
        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).
        assert (Rlt_bool
                  (Rabs
                     (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)))
                  (bpow radix2 1024) = true)
          by admit (* TODO: Overflow assumption *).
        rewrite H0 in *.
        destruct H.
        destruct H1.
        unfold b64_minus.
        rewrite H.
        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).
        { specialize (relative_error_binop (B2R 53 1024 b0 - B2R 53 1024 b1)%R).
          destruct relative_error_binop.
          - unfold bpow, radix2.
            (* This is the proof that we do not handle denormals! *) admit.
          - destruct H3.
            unfold round_mode; simpl in *.
            rewrite H4.
            unfold perturb.
            exists x.
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
            assert (2 ^ 53 = 9007199254740992)%positive by (cbv; auto).
            rewrite H5.
            lra. }
        { destruct H3 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; auto).
        specialize (IHe1 b0).
        specialize (IHe2 b1).
        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).
        assert (Rlt_bool
                  (Rabs
                     (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)))
                  (bpow radix2 1024) = true)
          by admit (* TODO: Overflow assumption *).
        rewrite H0 in *.
        destruct H.
        destruct H1.
        (*assert (is_nan 53 1024 (Bmult 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1) = false).
        { unfold is_nan, is_finite, b64_mult in *.
          destruct (Bmult 53 1024 eq_refl eq_refl binop_nan_pl64 mode_NE b0 b1); auto. } *)
        unfold b64_mult.
        rewrite H.
        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).
        { specialize (relative_error_binop (B2R 53 1024 b0 * B2R 53 1024 b1)%R).
          destruct relative_error_binop.
          - unfold bpow, radix2.
            (* This is the proof that we do not handle denormals! *) admit.
          - destruct H3.
            unfold round_mode; simpl in *.
            rewrite H4.
            unfold perturb.
            exists x.
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
            assert (2 ^ 53 = 9007199254740992)%positive by (cbv; auto).
            rewrite H5.
            lra. }
        { destruct H3 as [delta [delta_valid eval_perturb]].
          rewrite eval_perturb.
          econstructor; eauto.
          intros; congruence. }
      * assert (valid_div b1 = true) by (destruct (valid_div b1); congruence).
        rewrite H in H0.
        inversion H0; subst.
        assert (forall sgn, b1 = B754_infinity 53 1024 sgn -> False)
          by (intros; unfold valid_div in H; 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; auto).
        specialize (IHe1 b0).
        specialize (IHe2 b1).
        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).
        assert (Rlt_bool
                  (Rabs
                     (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)))
                  (bpow radix2 1024) = true)
          by admit (* TODO: Overflow assumption *).
        assert (B2R 53 1024 b1 <> 0%R).
        { hnf. intros. unfold B2R in H4. unfold valid_div in H. destruct b1; simpl in H; try congruence.
          unfold Fcore_defs.F2R in H4.
          simpl in H4.
          pose proof (bpow_gt_0 radix2 e).
          rewrite Z2R_cond_Zopp in H4.
          unfold cond_Ropp in H4.
          destruct b.
          - pose proof (Zgt_pos_0 m).
            rewrite Z.gt_lt_iff in H6.
            apply Z2R_lt in H6.
            apply Rmult_integral in H4.
            unfold Z2R in *.
            destruct H4; try lra.
          - pose proof (Zgt_pos_0 m).
            rewrite Z.gt_lt_iff in H6.
            apply Z2R_lt in H6.
            apply Rmult_integral in H4.
            unfold Z2R in *; destruct H4; lra. }
        rewrite H3 in H2.
        destruct H2; try auto.
        destruct H5.
        unfold b64_div.
        rewrite H2.
        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).
        { specialize (relative_error_binop (B2R 53 1024 b0 / B2R 53 1024 b1)%R).
          destruct relative_error_binop.
          - unfold bpow, radix2.
            (* This is the proof that we do not handle denormals! *) admit.
          - destruct H7.
            unfold round_mode; simpl in *.
                        rewrite H8.
            unfold perturb.
            exists x.
            split; try auto.
            eapply Rle_trans; eauto.
            unfold machineEpsilon.
            unfold Q2R.
            unfold Qnum, Qden.
            assert (2 ^ 53 = 9007199254740992)%positive by (cbv; auto).
            rewrite H9.
            lra. }
        { destruct H7 as [delta [delta_valid eval_perturb]].
          rewrite eval_perturb.
          econstructor; eauto. }
Admitted.