SimpleDoppler.v 15.7 KB
Newer Older
1
Require Import Coq.Reals.Reals Interval.Interval_tactic Coq.micromega.Psatz.
2
Require Import Coq.Setoids.Setoid.
3 4 5
Require Import Daisy.AbsoluteError Daisy.Commands Daisy.IntervalArith
        Daisy.Expressions Daisy.ErrorBounds
        Daisy.Infra.RealConstruction Daisy.Infra.Abbrevs Daisy.Infra.RealSimps.
6
(*
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
TODO: update according to:
[  Info  ]
[  Info  ] Starting specs preprocessing phase
[  Info  ] Finished specs preprocessing phase
[  Info  ]
[  Info  ]
[  Info  ] Starting range-error phase
[  Info  ] Machine epsilon 1.1102230246251565E-16
[  Info  ] 331.4: [331.4, 331.4],3.6792791036077685e-14
[  Info  ] u: [-100.0, 100.0],1.1102230246251565e-14
[  Info  ] (331.4 + u): [231.4, 431.4],9.579004256465851e-14
[  Info  ] Finished range-error phase
[  Info  ]
[  Info  ] Starting info phase
[  Info  ] doppler
[  Info  ] abs-error: 9.579004256465851e-14, range: [231.4, 431.4],
[  Info  ] rel-error: 4.139586973407887E-16
[  Info  ] Finished info phase
[  Info  ] time:
[  Info  ] info:      7 ms, rangeError:     67 ms, analysis:     11 ms, frontend:   2171 ms,



Heiko Becker's avatar
Heiko Becker committed
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
[  Info  ]
[  Info  ] Starting specs preprocessing phase
[  Info  ] Finished specs preprocessing phase
[  Info  ]
[  Info  ]
[  Info  ] Starting range-error phase
[  Info  ] Machine epsilon 1.1102230246251565E-16
[  Info  ] 331.4: [331.4, 331.4],7.358558207215537e-14
[  Info  ] u: [-100.0, 100.0],2.220446049250313e-14
[  Info  ] (331.4 + u): [231.4, 431.4],1.9158008512931703e-13
[  Info  ] Finished range-error phase
[  Info  ]
[  Info  ] Starting info phase
[  Info  ] doppler
[  Info  ] abs-error: 1.9158008512931703e-13, range: [231.4, 431.4],
[  Info  ] rel-error: 8.279173946815775E-16
[  Info  ] Finished info phase
[  Info  ] time:
[  Info  ] info:      9 ms, rangeError:    149 ms, analysis:     15 ms, frontend:   2709 ms,
49
 *)
Heiko Becker's avatar
Heiko Becker committed
50

Heiko Becker's avatar
Heiko Becker committed
51

52 53 54 55
Ltac prove_constant := unfold realFromNum, negativePower; interval.
Ltac rw_asm H Lem := rewrite Lem; rewrite Lem in H.


Heiko Becker's avatar
Heiko Becker committed
56 57 58 59 60 61 62 63
Lemma bounds_abs_maxAbs (a:R) (b:R) (x:R) :
  contained x (mkInterval a b) ->
  (Rabs x <= Rmax (Rabs a) (Rabs b))%R.
Proof.
  intros [a_le_x x_leq_b].
  apply RmaxAbs; auto.
Qed.

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
Lemma Rplus_minus_general:
  forall a b c,
    (a + b + (- a) + c = b + c)%R.
Proof.
  intros a b c.
  repeat rewrite Rplus_assoc.
  assert (b + (-a + c) = (b + c) + - a)%R.
  - rewrite Rplus_comm.
    rewrite Rplus_assoc.
    assert (c + b = b + c)%R by (rewrite Rplus_comm; auto).
    rewrite H.
    rewrite Rplus_comm; auto.
  - rewrite H.
    rewrite <- Rplus_assoc.
    rewrite (Rplus_assoc a (b + c)).
    rewrite <- Rsub_eq_Ropp_Rplus.
    rewrite Rplus_minus; auto.
Qed.

83 84 85 86 87
Definition u:nat := 1.
(** 1655/5 = 331; 0,4 = 2/5 **)
Definition cst1:R := 1657 / 5.

(** Define abbreviations **)
88
Definition varU:exp R := Param R u.
89 90 91 92
Definition valCst:exp R := Const cst1.
Definition valCstAddVarU:exp R := Binop Plus valCst varU.

(** Error values **)
93 94
Definition errCst1 := realFromNum 7358558207215537 15 14.
Definition errVaru := realFromNum 2220446049250313 15 14.
95
Definition lowerBoundAddUCst:R := 1157 / 5.
96
Definition upperBoundAddUCst:R := 2157 / 5.
97
Definition errAddUCst := realFromNum 19158008512931703 16 13.
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

(** The added assertion becomes the precondition for us **)
Definition precondition := fun env:nat -> R => (- 100 <= env u)%R /\ (env u <= 100)%R.

(** As stated, we need to define the absolute environment now as an inductive predicate **)
Inductive absEnv : abs_env :=
  theCst: absEnv valCst (mkInterval cst1 cst1) errCst1
|theVar: absEnv varU (mkInterval (- 100) (100)) errVaru
|theAddition: absEnv valCstAddVarU (mkInterval lowerBoundAddUCst upperBoundAddUCst) errAddUCst.

(** Show only the result for the expression as a test,
  since lifting to a program will only add additional predicate unfolding in this example
 **)
Goal forall (cenv:nat->R) (vR:R) (vF:R),
       precondition cenv ->
       eval_exp (0)%R cenv valCstAddVarU vR ->
114
       eval_exp machineEpsilon cenv valCstAddVarU vF ->
115
       AbsErrExp valCstAddVarU absEnv errAddUCst /\ (Rabs (vR - vF) <= errAddUCst)%R.
116
Proof.
117 118 119 120 121 122 123 124 125
  intros cenv vR vF precond eval_real eval_float.
  split.
  (* Proving Absolute Error correct *)
  - unfold valCstAddVarU.
    apply (AbsErrBinop _ _ (mkInterval cst1 cst1) errCst1
                       _ (mkInterval (- 100) (100)) errVaru
                       (mkInterval lowerBoundAddUCst upperBoundAddUCst) errAddUCst);
      [ constructor | constructor | constructor | | | ].
    + unfold valCst.
126
      apply (AbsErrConst cst1 (mkInterval cst1 cst1) errCst1); [constructor | ].
127
      unfold isSoundErr; simpl.
128
      unfold errCst1, cst1, machineEpsilon.
129
      assert (1657 / 5 >= 0)%R by prove_constant.
Heiko Becker's avatar
Heiko Becker committed
130 131 132 133
      unfold Rabs;
      destruct Rcase_abs as [lt_plus | ge_plus];
      [ exfalso; apply Rlt_not_le in lt_plus; apply lt_plus; apply Rge_le in H; auto | ].
      rewrite Rmax_left; [ | apply Req_le; auto].
Heiko Becker's avatar
Heiko Becker committed
134
      unfold realFromNum, negativePower.
Heiko Becker's avatar
Heiko Becker committed
135
      interval.
136
    + apply (AbsErrParam u (mkInterval (- 100) (100)) errVaru); [constructor | ].
137
      unfold isSoundErr; simpl.
138
      unfold machineEpsilon, errVaru.
139
      unfold realFromNum.
Heiko Becker's avatar
Heiko Becker committed
140
      unfold negativePower.
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
      assert (Rabs (-100) = 100%R) by (unfold Rabs; destruct Rcase_abs; lra).
      rewrite H.
      assert (Rabs 100 = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
      rewrite H0.
      rewrite Rmax_left; [ | apply Req_le; auto].
      interval.
    + unfold isSoundErrBin, valCst; split; try split; simpl; [unfold min4| unfold max4| ].
      * assert (Rmin (cst1 + - 100) (cst1 + 100) = (cst1 + - 100))%R by (unfold cst1; apply Rmin_left; interval).
        rewrite H.
        rewrite Rmin_comm in H.
        rewrite H.
        rewrite Rmin_left; [ | apply Req_le; auto].
        unfold cst1, lowerBoundAddUCst.
        field_simplify.
        apply Req_le; auto. (* FIXME: Why does interval not work here? *)
      * assert (Rmax (cst1 + - 100) (cst1 + 100) = (cst1 + 100))%R by (unfold cst1; apply Rmax_right; interval).
        assert (Rmax (cst1 + 100) (cst1 + 100) = cst1 + 100)%R by (apply Rmax_left; apply Req_le; auto). (* TODO: lemma *)
        rewrite H.
        rewrite H0.
        rewrite H.
        unfold cst1, upperBoundAddUCst.
        field_simplify.
        apply Req_le; auto.
      * unfold isSoundErr; simpl.
        unfold lowerBoundAddUCst, upperBoundAddUCst, errAddUCst.
166
        unfold machineEpsilon.
167 168 169 170
        assert (0 <= (1157/5))%R by interval.
        assert (0 <= (2157/5))%R by interval.
        repeat rewrite Rabs_pos_eq; auto.
        rewrite Rmax_right; [ | interval].
171
        prove_constant.
Heiko Becker's avatar
Notes  
Heiko Becker committed
172
  (* Semantic judgement  Rabs_triang_gen, Rabs_triang! triangle property? *)
173 174 175
  - inversion eval_real as
        [ | | | op e1 e2 v1 v2 delta delta_leq_0 eval_cst1 eval_param_u [op_eq e1_eq e2_eq] vR_eq];
      subst.
176
    rewrite perturb_0_val in eval_real; auto.
177 178 179 180 181 182
    inversion eval_cst1; subst.
    rewrite perturb_0_val in eval_real; auto.
    inversion eval_param_u; subst.
    repeat rewrite perturb_0_val; auto.
    rewrite perturb_0_val in eval_real , eval_param_u, eval_cst1; auto.
    clear delta delta0 delta1 delta_leq_0 H0 H1.
183
    inversion eval_float; subst.
184 185
    inversion H4; subst.
    inversion H5; subst.
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
    pose proof (const_abs_err_bounded _ _ _ _ eval_cst1 H4).
    pose proof (param_abs_err_bounded _ _ _ _ eval_param_u H5).
    assert (eval_exp machineEpsilon (updEnv 2 (perturb (cenv u) delta1) (updEnv 1 (perturb cst1 delta0) cenv))
    (Binop Plus (Var R 1) (Var R 2)) (perturb (eval_binop Plus (perturb cst1 delta0) (perturb (cenv u) delta1)) delta)) by admit.
    pose proof (add_abs_err_bounded _ _ _ _ _ _ _ _ _ _ _ eval_cst1 H4 eval_param_u H5 eval_real H6 H H3).
    eapply Rle_trans.
    apply H7.
    unfold perturb.
    repeat rewrite Rmult_plus_distr_l.
    repeat rewrite Rmult_1_r.
    eapply Rle_trans.
    eapply Rplus_le_compat_l.
    eapply Rmult_le_compat_r; [unfold machineEpsilon; prove_constant | ].
    eapply Rabs_triang.

    eapply Rle_trans.
    eapply Rplus_le_compat_l.
    eapply Rmult_le_compat_r;  [unfold machineEpsilon; prove_constant | ].
    eapply Rplus_le_compat.
    eapply Rabs_triang.
    eapply Rabs_triang.

    (* Now prove interval bounds *)
    unfold precondition in precond.
    destruct precond.
    assert (contained (cenv u) (mkInterval (-100) (100))) as uContained by (split; auto).
    assert (contained (cst1) (mkInterval cst1 cst1)) as cst1Contained by (apply validPointInterval).
    assert (contained (cenv u + cst1) (addInterval  (mkInterval (-100) (100)) (mkInterval cst1 cst1)))
      as containedAdd
        by (apply additionIsValid; auto).
    unfold contained in *; simpl in *.
    destruct uContained, cst1Contained, containedAdd.
    assert (Rabs cst1 = cst1) by (assert (0 <= cst1)%R by (unfold cst1; interval); rewrite Rabs_pos_eq; auto).
    repeat rewrite Rabs_mult.
    repeat rewrite H16.
    assert (Rabs (cenv u) <= Rmax (Rabs (-100)) (Rabs (100)))%R by (apply RmaxAbs; auto).
    assert (Rabs (-100) = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
    assert (Rabs 100 = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
    rewrite H18, H19 in H17.
    rewrite Rmax_left in H17; [ | lra].
    assert (forall eps:R, 0 <= eps -> Rabs (cenv u) * eps <= 100 * eps)%R by (intros; apply Rmult_le_compat_r; auto).

    eapply Rle_trans.
    eapply Rplus_le_compat.
    eapply Rplus_le_compat_l.
    eapply H20; auto.
    unfold machineEpsilon; prove_constant.
    eapply Rmult_le_compat_r.
    unfold machineEpsilon; prove_constant.
    eapply Rplus_le_compat.
    eapply Rplus_le_compat_l.
    eapply Rmult_le_compat_l; auto.
    unfold cst1; prove_constant.
    apply H0.
    eapply Rplus_le_compat.
    eapply H17.
    eapply Rle_trans.
    eapply H20.
    eapply Rabs_pos.
    eapply Rmult_le_compat_l.
    prove_constant.
    apply H1.
    unfold cst1, machineEpsilon, errAddUCst; prove_constant.
Heiko Becker's avatar
Heiko Becker committed
249
Admitted.
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

    (*
    assert (Rabs (cenv u) * machineEpsilon <= 100 * machineEpsilon)%R.
        {
          eapply H20; auto.
          unfold machineEpsilon; prove_constant.
        }
        {
          (* bound manipulation *)
          eapply Rle_trans.
          eapply Rplus_le_compat.
          apply Req_le; auto.
          eapply Rplus_le_compat.
          apply H19.
          unfold machineEpsilon; prove_constant.
          eapply Rplus_le_compat.
          eapply Rmult_le_compat_l.
          unfold cst1; prove_constant.
          apply H2.
          eapply Rplus_le_compat.
          apply H20.
          eapply Rplus_le_compat.
          eapply Rmult_le_compat; try auto using Rabs_pos.
          apply H16.
          apply H2.
          apply H22.
          unfold cst1, errAddUCst, machineEpsilon; prove_constant.
        }
Qed.


    unfold
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
    unfold eval_binop.
    unfold perturb in *; simpl in *.
    rewrite Rabs_err_simpl in *.
    repeat rewrite Rmult_plus_distr_l.
    repeat rewrite Rmult_1_r.
    repeat rewrite Rmult_plus_distr_r.
    rewrite Rsub_eq_Ropp_Rplus.
    repeat rewrite Ropp_plus_distr.
    repeat rewrite <- Rplus_assoc.
    rewrite Rplus_assoc at 1.
    repeat rewrite Rplus_minus_general.
    rewrite <- Rplus_assoc.
    assert (
        Rabs
          (- (cst1 * delta0) + - (cenv u * delta1) + - (cst1 * delta) + - (cst1 * delta0 * delta) + - (cenv u * delta) +
           - (cenv u * delta1 * delta)) <=
         Rabs (cst1 * delta0) + Rabs (cenv u * delta1) + Rabs (cst1 * delta) + Rabs (cst1 * delta0 * delta) +
         Rabs (cenv u * delta) + Rabs (cenv u * delta1 * delta)
        )%R.
  - eapply Rle_trans.
    eapply Rabs_triang.
    assert (Rabs (- (cst1 * delta0) + - (cenv u * delta1) + - (cst1 * delta) + - (cst1 * delta0 * delta) + - (cenv u * delta)) <= Rabs (- (cst1 * delta0)) + Rabs (- (cenv u * delta1)) + Rabs (- (cst1 * delta)) + Rabs (- (cst1 * delta0 * delta)) + Rabs (- (cenv u * delta)))%R.
    + eapply Rle_trans.
      eapply Rabs_triang.
      assert (Rabs (- (cst1 * delta0) + - (cenv u * delta1) + - (cst1 * delta) + - (cst1 * delta0 * delta)) <= Rabs (- (cst1 * delta0)) + Rabs (- (cenv u * delta1)) + Rabs (- (cst1 * delta)) + Rabs (- (cst1 * delta0 * delta)))%R.
      * eapply Rle_trans.
        eapply Rabs_triang.
        assert (Rabs (- (cst1 * delta0) + - (cenv u * delta1) + - (cst1 * delta)) <= Rabs (- (cst1 * delta0)) + Rabs (- (cenv u * delta1)) + Rabs (- (cst1 * delta)))%R.
        {
          eapply Rle_trans.
          eapply Rabs_triang.
          pose proof (Rabs_triang (- (cst1 * delta0)) (- (cenv u * delta1))).
          eapply Rplus_le_compat; [auto | apply Req_le; auto].
        }
        {
          eapply Rplus_le_compat; [auto | apply Req_le; auto].
        }
      * eapply Rplus_le_compat; [auto | apply Req_le; auto].
    + eapply Rplus_le_compat; [auto | apply Req_le; auto].
321
      repeat rewrite Rabs_Ropp in H6; auto.
322
      rewrite Rabs_Ropp; auto.
323
    + eapply Rle_trans. apply H6.
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345
      repeat rewrite Rplus_assoc.
      eapply Rle_trans.
      apply Rplus_le_compat. apply H.
      apply Req_le; auto.
      eapply Rle_trans.
      eapply Rplus_le_compat.
      eapply Req_le; auto.
      eapply Rplus_le_compat.
      apply H3.
      apply Req_le; auto.
      rewrite (Rabs_mult (cst1 * delta0) _).
      rewrite (Rabs_mult (cenv u * delta1) _).
      rewrite (Rabs_mult cst1 delta).
      rewrite (Rabs_mult (cenv u) delta).

      (* Now prove interval bounds *)
      unfold precondition in precond.
      destruct precond.
      assert (contained (cenv u) (mkInterval (-100) (100))) as uContained by (split; auto).
      assert (contained (cst1) (mkInterval cst1 cst1)) as cst1Contained by (apply validPointInterval).
      assert (contained (cenv u + cst1) (addInterval  (mkInterval (-100) (100)) (mkInterval cst1 cst1)))
        as containedAdd
346
          by (apply additionIsValid; auto).
347 348
      unfold contained in *; simpl in *.
      destruct uContained, cst1Contained, containedAdd.
Heiko Becker's avatar
Heiko Becker committed
349 350
      assert (Rabs cst1 = cst1) by (assert (0 <= cst1)%R by (unfold cst1; interval); rewrite Rabs_pos_eq; auto).
      repeat rewrite Rabs_mult.
351
      repeat rewrite H15.
Heiko Becker's avatar
Heiko Becker committed
352 353 354
      assert (Rabs (cenv u) <= Rmax (Rabs (-100)) (Rabs (100)))%R by (apply RmaxAbs; auto).
      assert (Rabs (-100) = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
      assert (Rabs 100 = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
355 356 357
      rewrite H17, H18 in H16.
      rewrite Rmax_left in H16; [ | lra].
      assert (forall eps:R, 0 <= eps -> Rabs (cenv u) * eps <= 100 * eps)%R by (intros; apply Rmult_le_compat_r; auto).
Heiko Becker's avatar
Heiko Becker committed
358 359 360
      assert (cst1 * Rabs delta0 * Rabs delta <= cst1 * machineEpsilon * machineEpsilon)%R.
      * assert (cst1 * Rabs delta0 <= cst1 * machineEpsilon)%R by (apply Rmult_le_compat_l; [unfold cst1, realFromNum, negativePower; interval | auto]).
        repeat rewrite Rmult_assoc.
361
        apply Rmult_le_compat_l; [unfold cst1; prove_constant | ].
Heiko Becker's avatar
Heiko Becker committed
362 363 364 365 366 367 368 369 370 371 372
        apply Rmult_le_compat; auto using Rabs_pos.
      * assert (Rabs (cenv u) * Rabs delta1 * Rabs delta <= Rabs (cenv u) * machineEpsilon * machineEpsilon)%R.
        repeat rewrite Rmult_assoc.
        apply Rmult_le_compat_l; try auto using Rabs_pos.
        apply Rmult_le_compat; auto using Rabs_pos.
        assert (Rabs (cenv u) * Rabs delta1 * Rabs delta <= 100 * machineEpsilon * machineEpsilon)%R.
        {
          eapply Rle_trans.
          apply H21.
          rewrite Rmult_assoc.
          rewrite (Rmult_assoc 100 _).
373
          apply Rmult_le_compat_r; [unfold machineEpsilon; prove_constant | auto ].
Heiko Becker's avatar
Heiko Becker committed
374 375 376 377 378 379 380 381
        }
        {
          (* bound manipulation *)
          eapply Rle_trans.
          eapply Rplus_le_compat.
          apply Req_le; auto.
          eapply Rplus_le_compat.
          apply H19.
382
          unfold machineEpsilon; prove_constant.
Heiko Becker's avatar
Heiko Becker committed
383 384
          eapply Rplus_le_compat.
          eapply Rmult_le_compat_l.
385
          unfold cst1; prove_constant.
Heiko Becker's avatar
Heiko Becker committed
386 387 388 389 390
          apply H2.
          eapply Rplus_le_compat.
          apply H20.
          eapply Rplus_le_compat.
          eapply Rmult_le_compat; try auto using Rabs_pos.
391
          apply H16.
Heiko Becker's avatar
Heiko Becker committed
392 393
          apply H2.
          apply H22.
394
          unfold cst1, errAddUCst, machineEpsilon; prove_constant.
Heiko Becker's avatar
Heiko Becker committed
395
        }
396
Qed. *)