Push to branch with affine arith changes

coq/Expressions.v

@@ -164,40 +164,48 @@ using a perturbation of the real valued computation by (1 + delta), where
 |delta| <= machine epsilon.
 **)
 Open Scope R_scope.
-Inductive eval_expr (E:env) (Gamma: nat -> option mType)
+
+Inductive eval_expr (E:env)
+          (Gamma: nat -> option mType)
+          (fBits:expr R -> option nat)
   :(expr R) -> R -> mType -> Prop :=
 | Var_load m x v:
     Gamma x = Some m ->
     E x = Some v ->
-    eval_expr E Gamma (Var R x) v m
+    eval_expr E Gamma fBits (Var R x) v m
 | Const_dist m n delta:
     Rabs delta <= mTypeToR m ->
-    eval_expr E Gamma (Const m n) (perturb n m delta) m
+    eval_expr E Gamma fBits (Const m n) (perturb n m delta) m
 | Unop_neg m f1 v1:
-    eval_expr E Gamma f1 v1 m ->
-    eval_expr E Gamma (Unop Neg f1) (evalUnop Neg v1) m
+    eval_expr E Gamma fBits f1 v1 m ->
+    eval_expr E Gamma fBits (Unop Neg f1) (evalUnop Neg v1) m
 | Unop_inv m f1 v1 delta:
     Rabs delta <= mTypeToR m ->
-    eval_expr E Gamma  f1 v1 m ->
+    eval_expr E Gamma fBits  f1 v1 m ->
     (~ v1 = 0)%R  ->
-    eval_expr E Gamma (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) m
+    eval_expr E Gamma fBits (Unop Inv f1) (perturb (evalUnop Inv v1) m delta) m
 | Downcast_dist m m1 f1 v1 delta:
     isMorePrecise m1 m = true ->
     Rabs delta <= mTypeToR m ->
-    eval_expr E Gamma f1 v1 m1 ->
-    eval_expr E Gamma (Downcast m f1) (perturb v1 m delta) m
-| Binop_dist m1 m2 op f1 f2 v1 v2 delta:
-    Rabs delta <= mTypeToR (join m1 m2) ->
-    eval_expr E Gamma f1 v1 m1 ->
-    eval_expr E Gamma f2 v2 m2 ->
+    eval_expr E Gamma fBits f1 v1 m1 ->
+    eval_expr E Gamma fBits (Downcast m f1) (perturb v1 m delta) m
+| Binop_dist m1 m2 op f1 f2 v1 v2 fBit delta m:
+    Rabs delta <= mTypeToR m ->
+    eval_expr E Gamma fBits f1 v1 m1 ->
+    eval_expr E Gamma fBits f2 v2 m2 ->
     ((op = Div) -> (~ v2 = 0)%R) ->
-    eval_expr E Gamma (Binop op f1 f2) (perturb (evalBinop op v1 v2) (join m1 m2) delta) (join m1 m2)
-| Fma_dist m1 m2 m3 f1 f2 f3 v1 v2 v3 delta:
-    Rabs delta <= mTypeToR (join3 m1 m2 m3) ->
-    eval_expr E Gamma f1 v1 m1 ->
-    eval_expr E Gamma f2 v2 m2 ->
-    eval_expr E Gamma f3 v3 m3 ->
-    eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3).
+    (isFixedPoint m1 -> isFixedPoint m2 -> fBits (Binop op f1 f2) = Some fBit) ->
+    join m1 m2 fBit = Some m ->
+    eval_expr E Gamma fBits (Binop op f1 f2) (perturb (evalBinop op v1 v2) m delta) m
+| Fma_dist m1 m2 m3 m f1 f2 f3 v1 v2 v3 fBit delta:
+    Rabs delta <= mTypeToR m ->
+    eval_expr E Gamma fBits f1 v1 m1 ->
+    eval_expr E Gamma fBits f2 v2 m2 ->
+    eval_expr E Gamma fBits f3 v3 m3 ->
+    (isFixedPoint m1 -> isFixedPoint m2 ->
+     isFixedPoint m3 -> fBits (Fma f1 f2 f3) = Some fBit) ->
+    join3 m1 m2 m3 fBit = Some m ->
+    eval_expr E Gamma fBits (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m.
 
 Close Scope R_scope.
 
@@ -206,78 +214,81 @@ Hint Constructors eval_expr.
 (** *)
 (*   Show some simpler (more general) rule lemmata *)
 (* **)
-Lemma Const_dist' m n delta v m' E Gamma:
+Lemma Const_dist' m n delta v m' E Gamma fBits:
   Rle (Rabs delta) (mTypeToR m') ->
   v = perturb n m delta ->
   m' = m ->
-  eval_expr E Gamma (Const m n) v m'.
+  eval_expr E Gamma fBits (Const m n) v m'.
 Proof.
   intros; subst; auto.
 Qed.
 
 Hint Resolve Const_dist'.
 
-Lemma Unop_neg' m f1 v1 v m' E Gamma:
-  eval_expr E Gamma f1 v1 m ->
+Lemma Unop_neg' m f1 v1 v m' E Gamma fBits:
+  eval_expr E Gamma fBits f1 v1 m ->
   v = evalUnop Neg v1 ->
   m' = m ->
-  eval_expr E Gamma (Unop Neg f1) v m'.
+  eval_expr E Gamma fBits (Unop Neg f1) v m'.
 Proof.
   intros; subst; auto.
 Qed.
 
 Hint Resolve Unop_neg'.
 
-Lemma Unop_inv' m f1 v1 delta v m' E Gamma:
+Lemma Unop_inv' m f1 v1 delta v m' E Gamma fBits:
   Rle (Rabs delta) (mTypeToR m') ->
-  eval_expr E Gamma  f1 v1 m ->
+  eval_expr E Gamma fBits  f1 v1 m ->
   (~ v1 = 0)%R  ->
   v = perturb (evalUnop Inv v1) m delta ->
   m' = m ->
-  eval_expr E Gamma (Unop Inv f1) v m'.
+  eval_expr E Gamma fBits (Unop Inv f1) v m'.
 Proof.
   intros; subst; auto.
 Qed.
 
 Hint Resolve Unop_inv'.
 
-Lemma Downcast_dist' m m1 f1 v1 delta v m' E Gamma:
+Lemma Downcast_dist' m m1 f1 v1 delta v m' E Gamma fBits:
   isMorePrecise m1 m = true ->
   Rle (Rabs delta) (mTypeToR m') ->
-  eval_expr E Gamma f1 v1 m1 ->
+  eval_expr E Gamma fBits f1 v1 m1 ->
   v = (perturb v1 m delta) ->
   m' = m ->
-  eval_expr E Gamma (Downcast m f1) v m'.
+  eval_expr E Gamma fBits (Downcast m f1) v m'.
 Proof.
   intros; subst; eauto.
 Qed.
 
 Hint Resolve Downcast_dist'.
 
-Lemma Binop_dist' m1 m2 op f1 f2 v1 v2 delta v m' E Gamma:
+Lemma Binop_dist' m1 m2 op f1 f2 v1 v2 delta v m' E Gamma fBits fBit:
   Rle (Rabs delta) (mTypeToR m') ->
-  eval_expr E Gamma f1 v1 m1 ->
-  eval_expr E Gamma f2 v2 m2 ->
+  eval_expr E Gamma fBits f1 v1 m1 ->
+  eval_expr E Gamma fBits f2 v2 m2 ->
   ((op = Div) -> (~ v2 = 0)%R) ->
   v = perturb (evalBinop op v1 v2) m' delta ->
-  m' = join m1 m2 ->
-  eval_expr E Gamma (Binop op f1 f2) v m'.
+  Some m' = join m1 m2 fBit ->
+  (isFixedPoint m1 -> isFixedPoint m2 -> fBits (Binop op f1 f2) = Some fBit) ->
+  eval_expr E Gamma fBits (Binop op f1 f2) v m'.
 Proof.
-  intros; subst; auto.
+  intros; subst; eauto.
 Qed.
 
 Hint Resolve Binop_dist'.
 
-Lemma Fma_dist' m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m' E Gamma:
+Lemma Fma_dist' m1 m2 m3 f1 f2 f3 v1 v2 v3 delta v m' E Gamma fBits fBit:
   Rle (Rabs delta) (mTypeToR m') ->
-  eval_expr E Gamma f1 v1 m1 ->
-  eval_expr E Gamma f2 v2 m2 ->
-  eval_expr E Gamma f3 v3 m3 ->
+  eval_expr E Gamma fBits f1 v1 m1 ->
+  eval_expr E Gamma fBits f2 v2 m2 ->
+  eval_expr E Gamma fBits f3 v3 m3 ->
   v = perturb (evalFma v1 v2 v3) m' delta ->
-  m' = join3 m1 m2 m3 ->
-  eval_expr E Gamma (Fma f1 f2 f3) v m'.
+  Some m' = join3 m1 m2 m3 fBit ->
+  (isFixedPoint m1 -> isFixedPoint m2 ->
+   isFixedPoint m3 -> fBits (Fma f1 f2 f3) = Some fBit) ->
+  eval_expr E Gamma fBits (Fma f1 f2 f3) v m'.
 Proof.
-  intros; subst; auto.
+  intros; subst; eauto.
 Qed.
 
 Hint Resolve Fma_dist'.
@@ -296,10 +307,10 @@ Fixpoint usedVars (V:Type) (e:expr V) :NatSet.t :=
   | _ => NatSet.empty
   end.
 
-Lemma toRMap_eval_REAL f v E Gamma m:
-  eval_expr E (toRMap Gamma) (toREval f) v m -> m = REAL.
+Lemma toRMap_eval_REAL f v E Gamma fBits m:
+  eval_expr E (toRMap Gamma) fBits (toREval f) v m -> m = REAL.
 Proof.
-  revert v E Gamma m.
+  revert v E Gamma fBits m.
   induction f; intros * eval_f; inversion eval_f; subst;
   repeat
     match goal with
@@ -314,7 +325,8 @@ Proof.
       by (eapply IHf1; eauto).
     assert (m2 = REAL)
       by (eapply IHf2; eauto);
-      subst; auto.
+      subst; eauto.
+    cbn in *; congruence.
   - assert (m1 = REAL)
       by (eapply IHf1; eauto).
     assert (m2 = REAL)
@@ -322,6 +334,7 @@ Proof.
     assert (m3 = REAL)
       by (eapply IHf3; eauto);
       subst; auto.
+    cbn in *; congruence.
 Qed.
 
 (**
@@ -338,10 +351,10 @@ Qed.
 (**
 Evaluation with 0 as machine epsilon is deterministic
 **)
-Lemma meps_0_deterministic (f:expr R) (E:env) Gamma:
+Lemma meps_0_deterministic (f:expr R) (E:env) Gamma fBits:
   forall v1 v2,
-  eval_expr E (toRMap Gamma) (toREval f) v1 REAL ->
-  eval_expr E (toRMap Gamma) (toREval f) v2 REAL ->
+  eval_expr E (toRMap Gamma) fBits (toREval f) v1 REAL ->
+  eval_expr E (toRMap Gamma) fBits (toREval f) v2 REAL ->
   v1 = v2.
 Proof.
   induction f;
@@ -385,41 +398,59 @@ For each evaluation of using an arbitrary epsilon, we can replace it by
 evaluating the subexprressions and then binding the result values to different
 variables in the Environment.
  **)
-Lemma binary_unfolding b f1 f2 E v1 v2 m1 m2 Gamma delta:
+Lemma binary_unfolding b f1 f2 E v1 v2 m1 m2 m Gamma fBits fBit delta:
   (b = Div -> ~(v2 = 0 )%R) ->
-  (Rabs delta <= (mTypeToR (join m1 m2)))%R ->
-  eval_expr E Gamma f1 v1 m1 ->
-  eval_expr E Gamma f2 v2 m2 ->
-  eval_expr E Gamma (Binop b f1 f2) (perturb (evalBinop b v1 v2) (join m1 m2) delta) (join m1 m2) ->
+  (isFixedPoint m1 -> isFixedPoint m2 -> fBits (Binop b f1 f2) = Some fBit) ->
+  (join m1 m2 fBit = Some m) ->
+  (Rabs delta <= mTypeToR m)%R ->
+  eval_expr E Gamma fBits f1 v1 m1 ->
+  eval_expr E Gamma fBits f2 v2 m2 ->
+  eval_expr E Gamma fBits (Binop b f1 f2) (perturb (evalBinop b v1 v2) m delta) m ->
   eval_expr (updEnv 2 v2 (updEnv 1 v1 emptyEnv))
            (updDefVars 2 m2 (updDefVars 1 m1 Gamma))
-             (Binop b (Var R 1) (Var R 2)) (perturb (evalBinop b v1 v2) (join m1 m2) delta) (join m1 m2).
+           (fun e =>
+              match e with
+              | Binop b (Var _ 1) (Var _ 2) => fBits (Binop b f1 f2)
+              | _ => fBits e
+           end)
+             (Binop b (Var R 1) (Var R 2)) (perturb (evalBinop b v1 v2) m delta) m.
 Proof.
   intros no_div_zero eval_f1 eval_f2 eval_float.
-  econstructor; try auto.
+  econstructor; try eauto.
+  - eapply Var_load; cbn; auto.
+  - eapply Var_load; cbn; auto.
 Qed.
 
-Lemma fma_unfolding f1 f2 f3 E v1 v2 v3 m1 m2 m3 Gamma delta:
-  (Rabs delta <= (mTypeToR (join3 m1 m2 m3)))%R ->
-  eval_expr E Gamma f1 v1 m1 ->
-  eval_expr E Gamma f2 v2 m2 ->
-  eval_expr E Gamma f3 v3 m3 ->
-  eval_expr E Gamma (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3) ->
+Lemma fma_unfolding f1 f2 f3 E v1 v2 v3 m1 m2 m3 m Gamma fBits fBit delta:
+  (isFixedPoint m1 -> isFixedPoint m2 ->
+   isFixedPoint m3 -> fBits (Fma f1 f2 f3) = Some fBit) ->
+  Some m = join3 m1 m2 m3 fBit ->
+  (Rabs delta <= mTypeToR m)%R ->
+  eval_expr E Gamma fBits f1 v1 m1 ->
+  eval_expr E Gamma fBits f2 v2 m2 ->
+  eval_expr E Gamma fBits f3 v3 m3 ->
+  eval_expr E Gamma fBits (Fma f1 f2 f3) (perturb (evalFma v1 v2 v3) m delta) m ->
   eval_expr (updEnv 3 v3 (updEnv 2 v2 (updEnv 1 v1 emptyEnv)))
            (updDefVars 3 m3 (updDefVars 2 m2 (updDefVars 1 m1 Gamma)))
-             (Fma (Var R 1) (Var R 2) (Var R 3)) (perturb (evalFma v1 v2 v3) (join3 m1 m2 m3) delta) (join3 m1 m2 m3).
+           (fun e =>
+              match e with
+              | Fma (Var _ 1) (Var _ 2) (Var _ 3) => fBits (Fma f1 f2 f3)
+              | _ => fBits e
+           end)
+             (Fma (Var R 1) (Var R 2) (Var R 3)) (perturb (evalFma v1 v2 v3) m delta) m.
 Proof.
-  econstructor; try auto.
+  econstructor; try eauto;
+    eapply Var_load; cbn; auto.
 Qed.
 
 Lemma eval_eq_env e:
-  forall E1 E2 Gamma v m,
+  forall E1 E2 Gamma fBits v m,
     (forall x, E1 x = E2 x) ->
-    eval_expr E1 Gamma e v m ->
-    eval_expr E2 Gamma e v m.
+    eval_expr E1 Gamma fBits e v m ->
+    eval_expr E2 Gamma fBits e v m.
 Proof.
   induction e; intros;
-    (match_pat (eval_expr _ _ _ _ _) (fun H => inversion H; subst; simpl in *));
+    (match_pat (eval_expr _ _ _ _ _ _) (fun H => inversion H; subst; simpl in *));
     try eauto.
   eapply Var_load; auto.
   rewrite <- (H n); auto.

coq/Infra/MachineType.v

@@ -19,7 +19,33 @@ From Flover
    the base field.
  **)
 Inductive mType: Type := REAL | M16 | M32 | M64
-                         | F (w:positive) (f:positive). (*| M128 | M256*)
+                         | F (w:positive) (f:nat). (*| M128 | M256*)
+
+Definition isFixedPoint m :Prop :=
+  match m with
+  |F _ _ => True
+  | _ => False
+  end.
+
+Definition isFixedPointB m :bool :=
+  match m with
+  |F _ _ => true
+  | _ => false
+  end.
+
+Lemma isFixedPoint_isFixedPointB m :
+  isFixedPoint m <-> isFixedPointB m = true.
+Proof.
+  destruct m; unfold isFixedPoint; simpl; split; try congruence; try auto.
+Qed.
+
+Lemma isFixedPoint_isFixedPointB_false m:
+  (~ isFixedPoint m) <-> isFixedPointB m = false.
+Proof.
+  destruct m; unfold isFixedPoint; simpl; split; try congruence; try auto.
+  intros.
+  exfalso; auto.
+Qed.
 
 (**
   Compute a machine epsilon from a  machine types in Q
@@ -30,7 +56,7 @@ Definition mTypeToR (m:mType) :R :=
   | M16 => 1 / 2^11
   | M32 => 1/ 2^24
   | M64 => 1/ 2^53
-  | F w f => 1/ 2^(Pos.to_nat f)
+  | F w f => 1/ 2^f
   (*
   (* the epsilons below match what is used internally in Daisy,
   although these value do not match the IEEE standard *)
@@ -44,7 +70,7 @@ Definition mTypeToQ (m:mType) :Q :=
   | M16 => (Qpower (2#1) (Zneg 11))
   | M32 => (Qpower (2#1) (Zneg 24))
   | M64 => (Qpower (2#1) (Zneg 53))
-  | F w f => Qpower (2#1) (Zneg f)
+  | F w f => Qpower (2#1) (- Z.of_nat f)
   (*
   (* the epsilons below match what is used internally in Daisy,
   although these value do not match the IEEE standard *)
@@ -77,24 +103,32 @@ Lemma mTypeToQ_mTypeToR m :
   Q2R (mTypeToQ m) = mTypeToR m.
 Proof.
   destruct m; cbn; try auto using Q2R0_is_0; try (unfold Q2R; simpl; lra).
+  pose proof (Qpower_opp (2#1) (Z.of_nat f)) as Qpower_eq.
+  apply Qeq_eqR in Qpower_eq.
+  rewrite Qpower_eq.
   rewrite Q2R_inv.
   - unfold Rdiv. rewrite Rmult_1_l.
     f_equal.
-    rewrite Qpower_decomp.
-    unfold Q2R; simpl.
-    rewrite Zpower_pos_powerRZ.
-    rewrite pow_powerRZ.
-    rewrite positive_nat_Z.
-    assert (IZR (Z.pos (1 ^ f)) = 1%R) as pow_1.
-    { unfold IZR.
-      f_equal.
-      rewrite Pos_pow_1_l; auto. }
-    rewrite pow_1.
-    rewrite Rinv_1, Rmult_1_r; auto.
+    destruct (Z.of_nat f) eqn:z_nat; try (destruct f; cbn in z_nat; congruence).
+    + cbn. destruct f; simpl in *; try congruence.
+      lra.
+    + unfold Qpower.
+      rewrite Qpower_decomp.
+      unfold Q2R; simpl.
+      rewrite Zpower_pos_powerRZ.
+      rewrite pow_powerRZ.
+      rewrite <- z_nat.
+      assert (IZR (Z.pos (1 ^ p)) = 1%R) as pow_1.
+      { unfold IZR.
+        f_equal.
+        rewrite Pos_pow_1_l; auto. }
+      rewrite pow_1.
+      rewrite Rinv_1, Rmult_1_r; auto.
   - hnf; intros power_0.
+    destruct f; simpl in *; try lra.
     rewrite Qpower_decomp in *.
     unfold Qeq in *; simpl in *.
-    pose proof (Zpow_facts.Zpower_pos_pos 2 f) as Zpower_pos.
+    pose proof (Zpow_facts.Zpower_pos_pos 2 (Pos.of_succ_nat f)) as Zpower_pos.
     assert (0 <2)%Z as pos2 by (omega).
     specialize (Zpower_pos pos2).
     rewrite Z.mul_1_r in *.
@@ -143,6 +177,8 @@ Lemma mTypeToQ_pos_Q m:
   0 <= mTypeToQ m.
 Proof.
   destruct m; simpl; cbn; try lra.
+  unfold Qpower.
+  destruct f; simpl in *; try lra.
   apply Qinv_le_0_compat.
   apply Qpower_pos_positive.
   lra.
@@ -163,7 +199,7 @@ Definition mTypeEq (m1:mType) (m2:mType) :=
   | M16, M16 => true
   | M32, M32 => true
   | M64, M64 => true
-  | F w1 f1, F w2 f2 => (w1 =? w2)%positive && (f1 =? f2)%positive
+  | F w1 f1, F w2 f2 => (w1 =? w2)%positive && (f1 =? f2)%nat
   (* | M128, M128 => true *)
   (* | M256, M256 => true *)
   | _, _ => false
@@ -173,7 +209,7 @@ Lemma mTypeEq_refl (m:mType):
   mTypeEq m m = true.
 Proof.
   intros. destruct m; try auto; simpl.
-  repeat rewrite Pos.eqb_refl; auto.
+  rewrite Pos.eqb_refl, Nat.eqb_refl; auto.
 Qed.
 
 Lemma mTypeEq_compat_eq(m1:mType) (m2:mType):
@@ -185,6 +221,7 @@ Proof.
         try congruence; try auto;
           try simpl in eq_case; try inversion eq_case.
   - andb_to_prop eq_case. f_equal; auto using Peqb_true_eq.
+    rewrite <- Nat.eqb_eq; auto.
   - inversion m2_case. apply mTypeEq_refl.
 Qed.
 
@@ -196,8 +233,8 @@ Proof.
     congruence.
   - case_eq m1; intros; case_eq m2; intros; subst; simpl; try congruence.
     destruct (w =? w0)%positive eqn:?; try auto.
-    destruct (f =? f0)%positive eqn:?; try auto.
-    rewrite Pos.eqb_eq in *; subst. congruence.
+    destruct (f =? f0)%nat eqn:?; try auto.
+    rewrite Pos.eqb_eq, Nat.eqb_eq in *; subst. congruence.
 Qed.
 
 Ltac type_conv :=
@@ -213,7 +250,7 @@ Proof.
   intros.
   destruct m1, m2; simpl; auto.
   rewrite Pos.eqb_sym; f_equal.
-  apply Pos.eqb_sym.
+  apply Nat.eqb_sym.
 Qed.
 
 (**
@@ -229,7 +266,7 @@ Qed.
 Definition isMorePrecise (m1:mType) (m2:mType) :=
   match m1, m2 with
   |REAL, _ => true
-  | F w1 f1, F w2 f2 => (w1 <=? w2)%positive (*&& (f1 <=? f2)%positive *)
+  | F w1 f1, F w2 f2 => (w2 <=? w1)%positive (*&& (f1 <=? f2)%positive *)
   | F w f, _ => false
   | _ , F w f => false
   | _, _ => Qle_bool (mTypeToQ m1) (mTypeToQ m2)
@@ -240,7 +277,7 @@ Definition isMorePrecise (m1:mType) (m2:mType) :=
 Definition morePrecise (m1:mType) (m2:mType) :=
   match m1, m2 with
   | REAL, _ => true
-  | F w1 f1, F w2 f2 => (w1 <=? w2)%positive (*&& (f1 <=? f2)%positive*)
+  | F w1 f1, F w2 f2 => (w2 <=? w1)%positive (*&& (f1 <=? f2)%positive*)
   | _ , F w f => false
   | F w f, _ => false
   | M16, M16 => true
@@ -251,20 +288,6 @@ Definition morePrecise (m1:mType) (m2:mType) :=
   | _, _ => false
   end.
 
-(* Lemma morePrecise_antisym m1 m2: *)
-(*   morePrecise m1 m2 = true -> *)
-(*   morePrecise m2 m1 = true -> *)
-(*   mTypeEq m1 m2 = true. *)
-(* Proof. *)
-(*   destruct m1; destruct m2; simpl; auto. *)
-(*   intros le_m1m2 le_m2m1. andb_to_prop le_m1m2. andb_to_prop le_m2m1. *)
-(*   apply Pos.leb_le in L; apply Pos.leb_le in R. *)
-(*   apply Pos.leb_le in L0; apply Pos.leb_le in R0. *)
-(*   split_bool; *)
-(*     apply Pos.eqb_eq; *)
-(*   apply Pos.le_antisym; auto. *)
-(* Qed. *)
-
 Lemma morePrecise_trans m1 m2 m3:
   morePrecise m1 m2 = true ->
   morePrecise m2 m3 = true ->
@@ -313,11 +336,20 @@ Qed.
   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) :=
-  if (morePrecise m1 m2) then m1 else m2.
+Definition join (m1:mType) (m2:mType) (fracBits:nat) :option mType:=
+  match m1, m2 with
+  | F w1 f1, F w2 f2 =>
+    if (w2 <=? w1)%positive
+    then Some (F w1 fracBits)
+    else Some (F w2 fracBits)
+  | F _ _, _ => None
+  | _ , F _ _ => None
+  | _ , _ => if (morePrecise m1 m2) then Some m1 else Some m2
+  end.
 
-Definition join3 (m1:mType) (m2:mType) (m3:mType) :=
-  join m1 (join m2 m3).
+Definition join3 (m1:mType) (m2:mType) (m3:mType) (fBits:nat):=
+  olet msub := (join m2 m3 fBits) in
+    join m1 msub fBits.
 
 (* Lemma REAL_join_is_REAL m1 m2: *)
 (*   join m1 m2 = REAL -> m1 = REAL /\ m2 = REAL. *)
@@ -327,7 +359,7 @@ Definition join3 (m1:mType) (m2:mType) (m3:mType) :=
 (*   destruct m1, m2; simpl in *; cbv in *; try congruence; try auto. *)
 (* Qed. *)
 
-Definition maxExponent (m:mType) :positive :=
+Definition maxExponent (m:mType) :nat :=
   match m with
   | REAL => 1
   | M16 => 15
@@ -336,7 +368,7 @@ Definition maxExponent (m:mType) :positive :=
   | F w f => f
   end.
 
-Definition minExponentPos (m:mType) :positive :=
+Definition minExponentPos (m:mType) :nat :=
   match m with
   | REAL => 1
   | M16 => 14
@@ -357,15 +389,15 @@ Fixed-Points:
 **)
 Definition maxValue (m:mType) :=
   match m with
-  |F w f => (((Z.pow_pos 2 (w -1))-1)#1) * Qinv ((Z.pow_pos 2 (maxExponent m))#1)
-  |_ => (Z.pow_pos 2 (maxExponent m)) # 1
+  |F w f => (((Z.pow_pos 2 (w -1))-1)#1) * Qinv ((Z.pow 2 (Z.of_nat (maxExponent m)))#1)
+  |_ => (Z.pow 2 (Z.of_nat (maxExponent m))) # 1
   end.
 
 (* Similarly: minimum values: we return 0 for fixed-point numbers here to make no fixed-point number be a denormal number ever*)
 Definition minValue_pos (m:mType) :=
   match m with
   |F w f => 0
-  | _ => Qinv ((Z.pow_pos 2 (minExponentPos m)) # 1)
+  | _ => Qinv ((Z.pow 2 (Z.of_nat (minExponentPos m))) # 1)
   end.
 
 (* (** Goldberg - Handbook of Floating-Point Arithmetic: (p.183) *)