SimpleMultiplication.v 12 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 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 157 158 159 160 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
Require Import Coq.Reals.Reals Interval.Interval_tactic Coq.micromega.Psatz.
Require Import Coq.Setoids.Setoid.
Require Import Daisy.AbsoluteError Daisy.Commands Daisy.IntervalArith
        Daisy.Expressions Daisy.ErrorBounds
        Daisy.Infra.RealConstruction Daisy.Infra.Abbrevs Daisy.Infra.RealSimps.

(*
[  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  ] 100.: [100.0, 100.0],0.
[  Info  ] u: [-100.0, 100.0],2.220446049250313e-14
[  Info  ] (100. * u): [-10000.0, 10000.0],4.440892098500627e-12
[  Info  ] Finished range-error phase
[  Info  ]
[  Info  ] Starting info phase
[  Info  ] doppler
[  Info  ] abs-error: 4.440892098500627e-12, range: [-10000.0, 10000.0],
[  Info  ] rel-error: NaN
[  Info  ] Finished info phase
[  Info  ] time:
[  Info  ] info:      6 ms, rangeError:     77 ms, analysis:     13 ms, frontend:   2400 ms,
 *)


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

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.

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.

Definition u:nat := 1.
(** 1655/5 = 331; 0,4 = 2/5 **)
Definition cst1:R := 100.

(** Define abbreviations **)
Definition varU:exp R := Param R u.
Definition valCst:exp R := Const cst1.
Definition valCstMultVarU:exp R := Binop Mult valCst varU.

(** Error values **)
Definition errCst1 := realFromNum 0 1 1.
Definition errVaru := realFromNum 2220446049250313 15 14.
Definition lowerBoundMultUCst:R := - realFromNum 10000 0 0.
Definition upperBoundMultUCst:R := realFromNum 10000 0 0.
Definition errMultUCst := realFromNum 4440892098500627 15 12.

(** 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 valCstMultVarU (mkInterval lowerBoundMultUCst upperBoundMultUCst) errMultUCst.


(** 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 valCstMultVarU vR ->
       eval_exp machineEpsilon cenv valCstMultVarU vF ->
       AbsErrExp valCstMultVarU absEnv errMultUCst /\
       (Rabs (vR - vF) <= errMultUCst)%R.
Proof.
  intros cenv vR vF precond eval_real eval_float.
  rewrite Rabs_simplify.
  split.
  (* Proving Absolute Error correct *)
  - unfold valCstMultVarU.
    apply (AbsErrBinop _ _ (mkInterval cst1 cst1) errCst1
                       _ (mkInterval (- 100) (100)) errVaru
                       (mkInterval lowerBoundMultUCst upperBoundMultUCst) errMultUCst);
      [ constructor | constructor | constructor | | | ].
    + unfold valCst.
      apply (AbsErrConst cst1 (mkInterval cst1 cst1) errCst1); [constructor | ].
      unfold isSoundErr; simpl.
      unfold errCst1, cst1, machineEpsilon.
      unfold realFromNum, negativePower.
      rewrite Rmax_left; [ |apply Req_le; auto].
      assert (Rabs 100 = 100)%R by (unfold Rabs; destruct Rcase_abs; lra).
      rewrite H.
      (* TODO: What about vlaues that are actually floats ? *) admit.
    + apply (AbsErrParam u (mkInterval (- 100) (100)) errVaru); [constructor | ].
      unfold isSoundErr; simpl.
      unfold machineEpsilon, errVaru.
      unfold realFromNum.
      unfold negativePower.
      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, lowerBoundMultUCst.
        prove_constant. (* Here it works... *)
      * 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, upperBoundMultUCst.
        field_simplify.
        prove_constant. (* And here too *)
      * unfold isSoundErr; simpl.
        unfold lowerBoundMultUCst, upperBoundMultUCst, errMultUCst.
        unfold machineEpsilon.
        assert (- realFromNum 10000 0 0 <= 0)%R by prove_constant.
        assert (0 <= realFromNum 10000 0 0) %R by prove_constant.
        rewrite Rabs_left1; auto.
        rewrite Rabs_pos_eq; auto.
        rewrite Ropp_involutive.
        rewrite Rmax_left; [ | apply Req_le; auto].
        prove_constant.

  (* Semantic judgement  Rabs_triang_gen, Rabs_triang! triangle property? *)
  - 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.
    rewrite perturb_0_val in eval_real; auto.
    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.
    inversion eval_float; subst.
    inversion H4; subst.
    inversion H5; subst.
    pose proof (const_abs_err_bounded _ _ _ _ eval_cst1 H4).
    pose proof (param_abs_err_bounded _ _ _ _ eval_param_u H5).
    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 (forall a b, a + - a + b = b)%R by (intros; lra).
    erewrite H6.
    repeat rewrite <- Rsub_eq_Ropp_Rplus.
    (** Lots of subtractions and we have positive numbers **)
    right.
    split.
    + eapply Rge_trans.
      eapply Fcore_Raux.Rabs_le_inv in H0; eapply Fcore_Raux.Rabs_le_inv in H2.
    (* TODO *)

    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].
      repeat rewrite Rabs_Ropp in H6; auto.
      rewrite Rabs_Ropp; auto.
    + eapply Rle_trans. apply H6.
      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
          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 H15.
      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 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).
      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.
        apply Rmult_le_compat_l; [unfold cst1; prove_constant | ].
        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 _).
          apply Rmult_le_compat_r; [unfold machineEpsilon; prove_constant | auto ].
        }
        {
          (* 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. *)