MachineType.v 5.31 KB
Newer Older
1 2 3 4 5 6
(**
  Implementation of machine-precision as a datatype for mixed-precision computations

  @author: Raphael Monat
  @maintainer: Heiko Becker
 **)
7

8 9 10
Require Import Coq.QArith.QArith Coq.QArith.Qminmax Coq.QArith.Qabs
               Coq.QArith.Qreals Coq.Reals.Reals Coq.micromega.Psatz.
Require Import Daisy.Infra.RealRationalProps.
11
(**
12
   Define machine precision as datatype
13
 **)
Heiko Becker's avatar
Heiko Becker committed
14
Inductive mType: Type := M0 | M16 | M32 | M64. (*| M128 | M256*)
15

16 17 18 19
(**
  Injection of machine types into Q
**)
Definition mTypeToQ (e:mType) :Q :=
20 21
  match e with
  | M0 => 0
Heiko Becker's avatar
Heiko Becker committed
22
  | M16 => (Qpower (2#1) (Zneg 11))
23
  | M32 => (Qpower (2#1) (Zneg 24))
24
  | M64 => (Qpower (2#1) (Zneg 53))
Heiko Becker's avatar
Heiko Becker committed
25
  (*
26 27 28
  (* the epsilons below match what is used internally in daisy,
  although these value do not match the IEEE standard *)
  | M128 => (Qpower (2#1) (Zneg 105))
Heiko Becker's avatar
Heiko Becker committed
29
  | M256 => (Qpower (2#1) (Zneg 211)) *)
30
  end.
31

32 33 34 35
Arguments mTypeToQ _/.

Lemma mTypeToQ_pos_Q:
  forall e, 0 <= mTypeToQ e.
36 37
Proof.
  intro e.
38 39
  case_eq e; intro;
    simpl mTypeToQ; apply Qle_bool_iff; auto.
40
Qed.
41

42 43
Lemma mTypeToQ_pos_R :
  forall e, (0 <= Q2R (mTypeToQ e))%R.
44 45
Proof.
  intro.
46
  rewrite <- Q2R0_is_0.
47
  apply Qle_Rle.
48
  apply mTypeToQ_pos_Q.
49
Qed.
='s avatar
= committed
50

51
Definition mTypeEq (m1:mType) (m2:mType) :=
52 53
  match m1, m2 with
  | M0, M0 => true
Heiko Becker's avatar
Heiko Becker committed
54
  | M16, M16 => true
55 56
  | M32, M32 => true
  | M64, M64 => true
Heiko Becker's avatar
Heiko Becker committed
57 58
  (* | M128, M128 => true *)
  (* | M256, M256 => true *)
59 60 61
  | _, _ => false
  end.

62 63
Lemma mTypeEq_compat_eq(m1:mType) (m2:mType):
  mTypeEq m1 m2 = true <-> m1 = m2.
64
Proof.
65 66 67 68 69
  split; intros eq_case;
    case_eq m1; intros m1_case;
      case_eq m2; intros m2_case; subst;
        try congruence; try auto;
          cbv in eq_case; inversion eq_case.
='s avatar
= committed
70
Qed.
71

72 73
Lemma mTypeEq_refl (m:mType):
  mTypeEq m m = true.
74
Proof.
75
  intros. rewrite mTypeEq_compat_eq; auto.
76 77
Qed.

78 79
Lemma mTypeEq_compat_eq_false (m1:mType) (m2:mType):
  mTypeEq m1 m2 = false <-> ~(m1 = m2).
80
Proof.
81 82 83 84
  split; intros eq_case.
  - hnf; intros; subst. rewrite mTypeEq_refl in eq_case.
    congruence.
  - case_eq m1; intros; case_eq m2; intros; subst; cbv; congruence.
85
Qed.
86

87 88 89 90 91 92 93
Ltac type_conv :=
  repeat
    (match goal with
     | [H : mTypeEq _ _ = true |- _] => rewrite mTypeEq_compat_eq in H; subst
     | [H : mTypeEq _ _ = false |- _] => rewrite mTypeEq_compat_eq_false in H
     end).

94 95
Lemma mTypeEq_sym (m1:mType) (m2:mType):
  forall b, mTypeEq m1 m2 = b -> mTypeEq m2 m1 = b.
96
Proof.
97 98 99 100
  intros.
  destruct b;
    [rewrite mTypeEq_compat_eq in * | rewrite mTypeEq_compat_eq_false in *];
    auto.
101 102
Qed.

103 104 105 106 107 108 109 110
(**
  Check if machine precision m1 is more precise than machine precision m2.
  M0 is real-valued evaluation, hence the most precise.
  All others are compared by
    mTypeToQ m1 <= mTypeToQ m2
 **)
Definition isMorePrecise (m1:mType) (m2:mType) :=
  Qle_bool (mTypeToQ m1) (mTypeToQ m2).
111

112 113
Lemma isMorePrecise_refl (m:mType) :
  isMorePrecise m m = true.
114
Proof.
115 116
  unfold isMorePrecise; simpl.
  apply Qle_bool_iff; lra.
117 118
Qed.

119 120
Lemma M0_least_precision (m:mType) :
  isMorePrecise m M0 = true -> m = M0.
121
Proof.
122 123 124
  intro m_morePrecise.
  unfold isMorePrecise in *.
  destruct m; cbv in m_morePrecise; congruence.
125
Qed.
126

127 128
Lemma M0_lower_bound (m:mType) :
  isMorePrecise M0 m = true.
129
Proof.
130
  destruct m; unfold isMorePrecise; cbv; auto.
131 132
Qed.

133 134 135 136 137 138
(**
  Function computing the join of two precisions, this is the most precise type,
  in which evaluation has to be performed, e.g. addition of 32 and 64 bit floats
  has to happen in 64 bits
**)
Definition join (m1:mType) (m2:mType) :=
139 140
  if (isMorePrecise m1 m2) then m1 else m2.

Nikita Zyuzin's avatar
Nikita Zyuzin committed
141 142 143
Definition join3 (m1:mType) (m2:mType) (m3:mType) :=
  join m1 (join m2 m3).

144 145 146 147 148 149
(* Lemma M0_join_is_M0 m1 m2: *)
(*   join m1 m2 = M0 -> m1 = M0 /\ m2 = M0. *)
(* Proof. *)
(*   unfold join. *)
(*   intros. *)
(*   destruct m1, m2; simpl in *; cbv in *; try congruence; try auto. *)
Heiko Becker's avatar
Heiko Becker committed
150 151
(* Qed. *)

Heiko Becker's avatar
Heiko Becker committed
152
Definition maxExponent (m:mType) :positive :=
Heiko Becker's avatar
Heiko Becker committed
153
  match m with
Heiko Becker's avatar
Heiko Becker committed
154
  | M0 => 1
Heiko Becker's avatar
Heiko Becker committed
155 156 157 158 159
  | M16 => 15
  | M32 => 127
  | M64 => 1023
  end.

Heiko Becker's avatar
Heiko Becker committed
160
Definition minExponentPos (m:mType) :positive :=
Heiko Becker's avatar
Heiko Becker committed
161
  match m with
Heiko Becker's avatar
Heiko Becker committed
162
  | M0 => 1
Heiko Becker's avatar
Heiko Becker committed
163 164 165 166 167 168 169 170 171 172 173 174
  | M16 => 14
  | M32 => 126
  | M64 => 1022
  end.

(**
Goldberg - Handbook of Floating-Point Arithmetic: (p.183)
(𝛃 - 𝛃^(1 - p)) * 𝛃^(e_max)

which simplifies to 2^(e_max) for base 2
**)
Definition maxValue (m:mType) :=
Heiko Becker's avatar
Heiko Becker committed
175
 (Z.pow_pos 2 (maxExponent m)) # 1.
Heiko Becker's avatar
Heiko Becker committed
176 177

Definition minValue (m:mType) :=
Heiko Becker's avatar
Heiko Becker committed
178
  Qinv ((Z.pow_pos 2 (minExponentPos m)) # 1).
Heiko Becker's avatar
Heiko Becker committed
179 180 181 182 183

(** Goldberg - Handbook of Floating-Point Arithmetic: (p.183)
  𝛃^(e_min -p + 1) = 𝛃^(e_min -1) for base 2
**)
Definition minDenormalValue (m:mType) :=
Heiko Becker's avatar
Heiko Becker committed
184
  Qinv (Z.pow_pos 2 (minExponentPos m - 1) # 1).
Heiko Becker's avatar
Heiko Becker committed
185 186 187 188

Definition normal (v:Q) (m:mType) :=
  Qle_bool (minValue m) (Qabs v) && Qle_bool (Qabs v) (maxValue m).

189
Definition denormal (v:Q) (m:mType) :=
190
  Qle_bool (Qabs v) (minValue m) && negb (Qeq_bool v 0).
191 192 193 194 195

Definition Normal (v:R) (m:mType) :=
  (Q2R (minValue m) <= (Rabs v) /\ (Rabs v) <= Q2R (maxValue m))%R.

Definition Denormal (v:R) (m:mType) :=
Heiko Becker's avatar
Heiko Becker committed
196
  ((Rabs v) < Q2R (minValue m) /\ ~ (v = 0))%R.
Heiko Becker's avatar
Heiko Becker committed
197 198 199 200 201 202
(**
  Predicate that is true if and only if the given value v is a valid
  floating-point value according to the the type m.
  Since we use the 1 + 𝝳 abstraction, the value must either be
  in normal range or 0.
**)
203
Definition validFloatValue (v:R) (m:mType) :=
Heiko Becker's avatar
Heiko Becker committed
204
  match m with
205 206
  | M0 => True
  | _ => Normal v m \/ Denormal v m \/ v = 0%R
Heiko Becker's avatar
Heiko Becker committed
207 208 209 210 211 212
  end.

Definition validValue (v:Q) (m:mType) :=
  match m with
  | M0 => true
  | _ => Qle_bool (Qabs v) (maxValue m)
Nikita Zyuzin's avatar
Nikita Zyuzin committed
213
  end.