interval_arith.v 10.1 KB
Newer Older
1
Require Import Coq.Reals.Reals.
2
Require Import Interval.Interval_tactic.
3
Require Import Flover.Infra.abbrevs.
4 5 6 7

Definition min4 (a:R) (b:R) (c:R) (d:R) := Rmin a (Rmin b (Rmin c d)).
Definition max4 (a:R) (b:R) (c:R) (d:R) := Rmax a (Rmax b (Rmax c d)).

Heiko Becker's avatar
Heiko Becker committed
8
(*
Heiko Becker's avatar
Heiko Becker committed
9
Definition intv_upd_err (op:R->R->R) (I1:intv) (e1:err) (I2:intv) (e2:err) :=
10 11 12 13 14
  (
    min4 (op (IVlo I1 - e1)%R (IVlo I2 - e2)%R)
         (op (IVlo I1 -e1)%R (IVhi I2 + e2)%R)
         (op (IVhi I1 + e1)%R (IVlo I2 - e2)%R)
         (op (IVhi I1 + e1)%R (IVlo I2 - e2)%R),
Heiko Becker's avatar
Heiko Becker committed
15 16
    max4 (op (IVlo I1 - e1)%R (IVlo I2 - e2)%R)
         (op (IVlo I1 - e1)%R (IVhi I2 + e2)%R)
17 18 19
         (op (IVhi I1 + e1)%R (IVlo I2 - e2)%R)
         (op (IVhi I1 + e1)%R (IVlo I2 - e2)%R)
  ).
Heiko Becker's avatar
Heiko Becker committed
20
 *)
21

Heiko Becker's avatar
Heiko Becker committed
22 23 24 25 26
Definition intv_upd (op:R->R->R) (I1:intv) (I2:intv) :=
  (
    min4 (op (IVlo I1) (IVlo I2))
         (op (IVlo I1) (IVhi I2))
         (op (IVhi I1) (IVlo I2))
Heiko Becker's avatar
Heiko Becker committed
27
         (op (IVhi I1) (IVhi I2)),
Heiko Becker's avatar
Heiko Becker committed
28 29 30
    max4 (op (IVlo I1) (IVlo I2))
         (op (IVlo I1) (IVhi I2))
         (op (IVhi I1) (IVlo I2))
Heiko Becker's avatar
Heiko Becker committed
31
         (op (IVhi I1) (IVhi I2))
Heiko Becker's avatar
Heiko Becker committed
32 33
  ).

Heiko Becker's avatar
Heiko Becker committed
34 35 36
Definition intv_widen_err (Iv:intv) (e:err) := ((IVlo Iv - e)%R, (IVhi Iv + e)%R).

Definition intv_neg (Iv:intv) := ((- IVhi Iv)%R, (- IVlo Iv)%R).
37

Heiko Becker's avatar
Heiko Becker committed
38 39
Definition intv_add (I1:intv) (I2:intv) :=
  intv_upd Rplus I1 I2.
40

Heiko Becker's avatar
Heiko Becker committed
41 42 43 44 45 46
Definition intv_sub (I1:intv) (I2:intv) :=
  intv_add I1 (intv_neg I2).

Definition intv_add_err (I1:intv) (e1:err) (I2:intv) (e2:err) :=
  intv_add (intv_widen_err I1 e1) (intv_widen_err I2 e2).

Heiko Becker's avatar
Heiko Becker committed
47
Definition intv_sub_err (I1:intv) (e1:err) (I2:intv) (e2:err) :=
Heiko Becker's avatar
Heiko Becker committed
48
  intv_sub (intv_widen_err I1 e1) (intv_widen_err I2 e2).
Heiko Becker's avatar
Heiko Becker committed
49

50 51 52
(*
  FIXME: This needs to be beautified
*)
Heiko Becker's avatar
Heiko Becker committed
53
Definition intv_mult_err (I1:intv) (e1:err) (I2:intv) (e2:err) :=
54 55 56
  intv_add
    (intv_add
       (intv_add
Heiko Becker's avatar
Heiko Becker committed
57 58 59 60 61 62 63 64 65
          (intv_upd Rmult I1 I2)
          (intv_upd Rmult I1 (e1,e1)))
       (intv_upd Rmult I2 (e2,e2)))
    (Rmult e1 e2, Rmult e1 e2).

Definition intv_mult (I1:intv) (I2:intv) :=
  intv_upd Rmult I1 I2.

Definition intv_div_err (I1:intv) (e1:err) (I2:intv) (e2:err) :=
Heiko Becker's avatar
Heiko Becker committed
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
  intv_mult_err I1 e1 ( (/ IVlo (I2))%R, (/ IVhi (I2))%R) e2.

Definition contained (Iv:intv) (v:R) := (IVlo Iv <= v)%R /\ (v <= IVhi Iv)%R.

Lemma IntvNegMono (Iv:intv) (a:R) :
  contained Iv a -> contained (intv_neg Iv) (- a).
Proof.
  unfold contained; destruct Iv as [ivlo ivhi]; simpl; intros [ivlo_le_a a_le_ivhi].
  split; apply Ropp_ge_le_contravar; apply Rle_ge; auto.
Qed.

Lemma IntvAddMono (I1:intv) (I2:intv) (a:R) (b:R):
  contained I1 a ->
  contained I2 b ->
  contained (intv_add I1 I2) (a + b).
Proof.
  unfold contained.
  unfold intv_add, intv_upd, min4, max4.
  intros contained_I1 contained_I2.
  destruct contained_I1 as [lo_leq_a a_leq_hi].
  destruct contained_I2 as [lo_leq_b b_leq_hi].
  destruct I1 as [ivlo1 ivhi1].
  destruct I2 as [ivlo2 ivhi2].
  unfold IVlo, IVhi.
  simpl; split.
  (** Lower Bound **)
  - assert (ivlo1 + ivlo2 <= a + b)%R by (apply Rplus_le_compat; auto).
    assert
      ((Rmin (ivlo1 + ivlo2)
             (Rmin (ivlo1 + ivhi2) (Rmin (ivhi1 + ivlo2) (ivhi1 + ivhi2))) <= (ivlo1 + ivlo2)))%R
      by apply Rmin_l.
    apply (Rle_trans _ (ivlo1 + ivlo2) _); auto.
    (** Upper Bound **)
  - assert (a + b <= ivhi1 + ivhi2)%R by (apply Rplus_le_compat; auto).
    assert
      (ivhi1 + ivhi2 <= (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)))%R by apply Rmax_r.
    assert (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2) <= (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2))))%R
      by apply Rmax_r.
    assert
      (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)) <=
       (Rmax (ivlo1 + ivlo2)
             (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)))))%R
      by apply Rmax_r.
    assert (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2) <=
            (Rmax (ivlo1 + ivlo2)
                  (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)))))%R
      by (apply (Rle_trans _  (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2))) _); auto).
    assert (ivhi1 + ivhi2 <=
            (Rmax (ivlo1 + ivlo2)
                  (Rmax (ivlo1 + ivhi2) (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)))))%R
      by (apply (Rle_trans _  (Rmax (ivhi1 + ivlo2) (ivhi1 + ivhi2)) _); auto).
    apply (Rle_trans _ (ivhi1 + ivhi2) _); auto.
Qed.

Corollary IntvSubMono (I1:intv) (I2:intv) (a:R) (b:R) :
  contained I1 a ->
  contained I2 b ->
  contained (intv_sub I1 I2) (a - b).
Proof.
  unfold intv_sub.
  intros contained_1 contained_I2.
  assert (a - b = a + (- b))%R.
  field_simplify; reflexivity.
  rewrite H.
  apply IntvAddMono; auto.
  apply IntvNegMono; auto.
Qed.

(**
Useful Lemmata:
Rinv_mult_simpl:
  forall r1 r2 r3 : R,
  r1 <> 0%R -> (r1 * / r2 * (r3 * / r1))%R = (r3 * / r2)%R
Rmult_lt_compat_r:
  forall r r1 r2 : R, (0 < r)%R -> (r1 < r2)%R -> (r1 * r < r2 * r)%R
Rmult_gt_compat_r:
  forall r r1 r2 : R, (r > 0)%R -> (r1 > r2)%R -> (r1 * r > r2 * r)%R
Rmult_gt_compat_l:
  forall r r1 r2 : R, (r > 0)%R -> (r1 > r2)%R -> (r * r1 > r * r2)%R
Rmult_le_compat_l:
  forall r r1 r2 : R, (0 <= r)%R -> (r1 <= r2)%R -> (r * r1 <= r * r2)%R
Rmult_le_compat_r:
  forall r r1 r2 : R, (0 <= r)%R -> (r1 <= r2)%R -> (r1 * r <= r2 * r)%R
Rmult_ge_compat_l:
  forall r r1 r2 : R, (r >= 0)%R -> (r1 >= r2)%R -> (r * r1 >= r * r2)%R
Rmult_ge_compat_r:
  forall r r1 r2 : R, (r >= 0)%R -> (r1 >= r2)%R -> (r1 * r >= r2 * r)%R
Rmult_le_compat:
  forall r1 r2 r3 r4 : R,
  (0 <= r1)%R ->
  (0 <= r3)%R -> (r1 <= r2)%R -> (r3 <= r4)%R -> (r1 * r3 <= r2 * r4)%R
Rmult_ge_compat:
  forall r1 r2 r3 r4 : R,
  (r2 >= 0)%R ->
  (r4 >= 0)%R -> (r1 >= r2)%R -> (r3 >= r4)%R -> (r1 * r3 >= r2 * r4)%R
Rmult_gt_0_lt_compat:
  forall r1 r2 r3 r4 : R,
  (r3 > 0)%R ->
  (r2 > 0)%R -> (r1 < r2)%R -> (r3 < r4)%R -> (r1 * r3 < r2 * r4)%R
Rmult_le_0_lt_compat:
  forall r1 r2 r3 r4 : R,
  (0 <= r1)%R ->
  (0 <= r3)%R -> (r1 < r2)%R -> (r3 < r4)%R -> (r1 * r3 < r2 * r4)%R
Rmult_lt_0_compat:
  forall r1 r2 : R, (0 < r1)%R -> (0 < r2)%R -> (0 < r1 * r2)%R
Rmult_gt_0_compat:
  forall r1 r2 : R, (r1 > 0)%R -> (r2 > 0)%R -> (r1 * r2 > 0)%R
Rmult_le_compat_neg_l:
  forall r r1 r2 : R, (r <= 0)%R -> (r1 <= r2)%R -> (r * r2 <= r * r1)%R
Rmult_le_ge_compat_neg_l:
  forall r r1 r2 : R, (r <= 0)%R -> (r1 <= r2)%R -> (r * r1 >= r * r2)%R
Rmult_lt_gt_compat_neg_l:
  forall r r1 r2 : R, (r < 0)%R -> (r1 < r2)%R -> (r * r1 > r * r2)%R
**)

181
(**
Heiko Becker's avatar
Heiko Becker committed
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
Lemma IntvMultMono (I1:intv) (I2:intv) (a:R) (b:R) :
  contained I1 a ->
  contained I2 b ->
  contained (intv_mult I1 I2) (a * b).
Proof.
  unfold contained, intv_mult, intv_upd, IVlo, IVhi;
  destruct I1 as [ivlo1 ivhi1]; destruct I2 as [ivlo2 ivhi2]; simpl.
  intros [ lo_leq_a a_leq_hi] [lo_leq_b b_leq_hi].
  assert (ivlo1 <= ivhi1)%R by (apply (Rle_trans _ a _); auto).
  assert (ivlo2 <= ivhi2)%R by (apply (Rle_trans _ b _); auto).
  destruct (Rle_dec ivhi1 0); destruct (Rle_dec ivhi2 0);
  [(* Both less do no further case distinction *)
    assert (ivlo1 <= 0)%R by (apply (Rle_trans _ (ivhi1) _); auto);
    assert (ivlo2 <= 0)%R by (apply (Rle_trans _ (ivhi2) _); auto)
  | (* Fst one less, destruct second one *)
  assert (ivlo1 <= 0)%R by (apply (Rle_trans _ (ivhi1) _); auto);
  destruct (Rle_dec ivlo2 0)
  | (* Opposite case *)
  assert (ivlo2 <= 0)%R by (apply (Rle_trans _ (ivhi2) _); auto);
  destruct (Rle_dec ivlo1 0)
  | (* Both greater -> both destruct *)
  destruct (Rle_dec ivlo2 0); destruct (Rle_dec ivlo1 0) ].
  -
  (*
  First subcase, all intervals negative --> ivhi1 * ivhi2 is maximum
  Rmult_le_compat_neg_l b a ivhi1
  Rmult_le_compat_neg_l ivlo1 b ivhi2
   *)
    assert (b <= 0)%R as b_leq_0 by (apply (Rle_trans _ (ivhi2) _); auto).
    split.
212 213 214 215 216
    + pose proof (Rmult_le_compat_neg_l b a ivhi1 b_leq_0 a_leq_hi).
      pose proof (Rmult_le_compat_neg_l ivhi1 b ivhi2 r b_leq_hi).
      rewrite Rmult_comm in H3.
      pose proof (Rle_trans _ _ _ H4 H3).
      rewrite (Rmult_comm a b).
Heiko Becker's avatar
Heiko Becker committed
217
      unfold min4.
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
      assert
        (Rmin (ivlo1 * ivlo2)
              (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) <= ivhi1 * ivhi2)%R.
      * assert (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) <= ivhi1 * ivhi2)%R.
        {
          assert (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2) <= ivhi1 * ivhi2)%R by apply Rmin_r.
          apply (Rle_trans _ (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) _); auto.
          apply Rmin_r.
        }
        {
          apply (Rle_trans _ (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) _); auto.
          apply Rmin_r.
        }
      * apply (Rle_trans _ (ivhi1 * ivhi2) _); auto.
    + unfold max4.
      apply Rge_le.
      pose proof (Rmult_le_ge_compat_neg_l ivlo1 ivlo2 b H1 lo_leq_b).
      pose proof (Rmult_le_ge_compat_neg_l b ivlo1 a b_leq_0 lo_leq_a).
      rewrite Rmult_comm in H4.
      pose proof (Rge_trans _ _ _ H3 H4).
      rewrite (Rmult_comm a b).
      apply (Rge_trans _ (ivlo1 * ivlo2) _); auto.
      apply Rle_ge, Rmax_l.
  - assert (a <= 0)%R as a_leq_0 by (apply (Rle_trans _ (ivhi1) _); auto).
    split.
    + apply Rnot_le_gt in n.
      apply Rgt_ge in n.
      pose proof (Rmult_le_compat_neg_l a b ivhi2).
      pose proof (Rmult_le_compat_neg_l ivhi1).
      rewrite Rmult_comm in H3.
      pose proof (Rle_trans _ _ _ H4 H3).
      rewrite (Rmult_comm a b).
      unfold min4.
      assert
        (Rmin (ivlo1 * ivlo2)
              (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) <= ivhi1 * ivhi2)%R.
      * assert (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) <= ivhi1 * ivhi2)%R.
        {
          assert (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2) <= ivhi1 * ivhi2)%R by apply Rmin_r.
          apply (Rle_trans _ (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) _); auto.
          apply Rmin_r.
        }
        {
          apply (Rle_trans _ (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) _); auto.
          apply Rmin_r.
        }
      * apply (Rle_trans _ (ivhi1 * ivhi2) _); auto.
    + unfold max4.
      apply Rge_le.
      pose proof (Rmult_le_ge_compat_neg_l ivlo1 ivlo2 b H1 lo_leq_b).
      pose proof (Rmult_le_ge_compat_neg_l b ivlo1 a b_leq_0 lo_leq_a).
      rewrite Rmult_comm in H4.
      pose proof (Rge_trans _ _ _ H3 H4).
      rewrite (Rmult_comm a b).
      apply (Rge_trans _ (ivlo1 * ivlo2) _); auto.
      apply Rle_ge, Rmax_l.
Heiko Becker's avatar
Heiko Becker committed
274 275 276 277 278 279 280 281
  -  .
  -  .
  -  .
  -  .
  -  .
  -  .
  -  .
**)