Library FNumeric

Add LoadPath "IEEE754" as IEEE754.
Require Export IEEE754.IEEE.
Require Import IEEE754.Diadic.
Require Export ZArith.
Require Import Numeric.
Open Local Scope Z_scope.

@author Aleksy Schubert
Common interface for IEEE 754
Module Type FNUMERIC.
  Parameter t : Set.
  Parameter toIEEE : tabstract_IEEE.
  Parameter precision : float_type.
  Definition rmode := Rounding_nearest.
  Inductive isNaN : abstract_IEEEProp :=
     | IsNaNIs : x, isNaN (NaN x)
   .
  Parameter toDiadic : abstract_IEEEdiadic.
  Parameter fromDiadic : diadicabstract_IEEE.

  Parameter add : ttt.
  Parameter add_prop : i1 i2,
    toIEEE (add i1 i2) = fst (abstract_add precision rmode (toIEEE i1) (toIEEE i2)).

  Parameter sub : ttt.
  Parameter sub_prop : i1 i2,
    toIEEE (sub i1 i2) = fst (abstract_sub precision rmode (toIEEE i1) (toIEEE i2)).

  Parameter mul : ttt.
  Parameter mul_prop : i1 i2,
    toIEEE (mul i1 i2)= fst (abstract_mult precision rmode (toIEEE i1) (toIEEE i2)).

  Parameter div : ttt.
  Parameter div_prop : i1 i2,
    toIEEE (div i1 i2)= fst (abstract_div precision rmode (toIEEE i1) (toIEEE i2)).

  Parameter neg : tt.
  Parameter neg_prop : f,
    toIEEE (neg f) = abstract_opp (toIEEE f).

The function defines the operation of comparison operators. According to the Java Virtual Machine Specification it is:
If value1' is greater than value2', the int value 1 is pushed onto the operand stack. Otherwise, if value1' is equal to value2', the int value 0 is pushed onto the operand stack. Otherwise, if value1' is less than value2', the int value -1 is pushed onto the operand stack. Otherwise, at least one of value1' or value2' is NaN. The fcmpg instruction pushes the int value 1 onto the operand stack and the fcmpl instruction pushes the int value -1 onto the operand stack.
  Parameter cmpg : ttZ.
  Parameter cmpg_gt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Gt, Nothing)cmpg f1 f2 = 1.
  Parameter cmpg_eq : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Eq, Nothing)cmpg f1 f2 = 0.
  Parameter cmpg_lt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Lt, Nothing)cmpg f1 f2 = -1.
  Parameter cmpg_nan : f1 f2,
    (isNaN (toIEEE f1)) (isNaN (toIEEE f2))cmpg f1 f2 = 1.

  Parameter cmpl : ttZ.
  Parameter cmpl_gt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Gt, Nothing)cmpl f1 f2 = 1.
  Parameter cmpl_eq : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Eq, Nothing)cmpl f1 f2 = 0.
  Parameter cmpl_lt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Lt, Nothing)cmpl f1 f2 = -1.
  Parameter cmpl_nan : f1 f2,
    (isNaN (toIEEE f1)) (isNaN (toIEEE f2))cmpl f1 f2 = -1.

The function defines the operation of floating point remainder. According to the Java Virtual Machine Specification it is:
The result of an frem instruction is governed by these rules: If either value1' or value2' is NaN, the result is NaN. If neither value1' nor value2' is NaN, the sign of the result equals the sign of the dividend. If the dividend is an infinity or the divisor is a zero or both, the result is NaN. If the dividend is finite and the divisor is an infinity, the result equals the dividend. If the dividend is a zero and the divisor is finite, the result equals the dividend.
In the remaining cases, where neither operand is an infinity, a zero, or NaN, the floating-point remainder result from a dividend value1' and a divisor value2' is defined by the mathematical relation result = value1' - (value2' * q), where q is an integer that is negative only if value1' / value2' is negative and positive only if value1' / value2' is positive, and whose magnitude is as large as possible without exceeding the magnitude of the true mathematical quotient of value1' and value2'.
  Parameter abstract_rem : rounding_modeabstract_IEEEabstract_IEEEabstract_IEEE × exception.

  Parameter rem : ttt.
  Parameter rem_prop : i1 i2,
    toIEEE (rem i1 i2) = fst (abstract_rem rmode (toIEEE i1) (toIEEE i2)).

  Parameter zero: t.
  Parameter zero_prop : toIEEE zero = Zero true.

  Parameter one: t.
  Parameter one_prop : toIEEE one = Normal true (dig_length precision) (Z.to_pos 1).

  Parameter two: t.
  Parameter two_prop : toIEEE two = Normal true (1 + (dig_length precision)) (Z.to_pos 1).

  Definition GtIEEE : abstract_IEEEabstract_IEEEProp :=
    fun x y(abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Gt, Nothing).

  Definition GeIEEE : abstract_IEEEabstract_IEEEProp :=
    fun x y(abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Gt, Nothing)
      (abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Eq, Nothing).

The function defines the operation of floating point truncation to an integer. Its properties according to the Java Virtual Machine Specification are mentioned along the paramters below.
The first parameter is the floating point value to truncate, the second parameter is the number of bits of the resulting integer minus 1 (it corresponds do half_base in Numeric).

  Parameter truncate : tZZ.
If the value' is NaN, the result of the conversion is 0.
  Parameter truncate_nan :
    f1 n, ( x, toIEEE f1 = NaN x) → truncate f1 n = 0.

Otherwise, if the value' is not an infinity, it is rounded to an integer value V, rounding towards zero using IEEE 754 round towards zero mode. If this integer value V can be represented as an int, then the result is the int value V.
  Parameter truncate_range_positive :
    f half_base,
     let f1 := toIEEE f in
     let maxint := abstract_of_diadic precision Rounding_sup
                    (Diadic (2^half_base) 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
     let (truncp1, _) :=
           (abstract_add precision Rounding_sup trunc (toIEEE one)) in
       ( (GtIEEE maxint f1) (GeIEEE f1 (Zero true))) →
         ( (GtIEEE truncp1 f1) (GeIEEE f1 trunc)).

  Parameter truncate_range_negative :
    f half_base,
     let f1 := toIEEE f in
     let minint := abstract_of_diadic precision Rounding_sup
                    (Diadic (2^half_base) 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
     let (truncm1, _) :=
           (abstract_sub precision Rounding_inf trunc (toIEEE one)) in
       ( (GeIEEE f1 minint) (GeIEEE (Zero true) f1)) →
         ( (GtIEEE f1 truncm1) (GeIEEE trunc f1)).

Otherwise, either the value' must be too small (a negative value of large magnitude or negative infinity), and the result is the smallest representable value of type integer,
  Parameter truncate_too_small :
    f half_base,
     let f1 := toIEEE f in
     let minintZ := - 2^half_base in
     let minint := abstract_of_diadic precision Rounding_sup
                     (Diadic minintZ 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
       (GtIEEE minint f1 (f1 = Infty false)) →
        truncate f half_base = minintZ.

...or the value' must be too large (a positive value of large magnitude or positive infinity), and the result is the largest representable value of type integer.
  Parameter truncate_too_large :
    f half_base,
     let f1 := toIEEE f in
     let maxintZ := 2^half_base - 1 in
     let maxint := abstract_of_diadic precision Rounding_sup
                     (Diadic maxintZ 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
       (GeIEEE f1 maxint (f1 = Infty true)) →
        truncate f half_base = maxintZ.

  Parameter const : abstract_IEEEt.

  Parameter iconst : Zt.

End FNUMERIC.

Module Type FPRECISION.
  Parameter precision : float_type.
End FPRECISION.

Module FMake (S:FPRECISION) <: FNUMERIC with Definition precision:=S.precision.
  Definition precision := S.precision.
  Definition t := abstract_IEEE.
  Definition toIEEE : tabstract_IEEE := fun xx.
  Definition rmode := Rounding_nearest.
  Inductive isNaN : abstract_IEEEProp :=
     | IsNaNIs : x, isNaN (NaN x)
   .

This requires amendment since structure must be preserved
  Definition toDiadic (num : abstract_IEEE) := diadic_of_abstract precision num.

  Definition fromDiadic (num : diadic) := abstract_of_diadic precision rmode num.

TODO: gradual underflow TODO: value set conversion

  Definition add x y := fst (abstract_add precision rmode x y).
  Lemma add_prop : i1 i2,
    toIEEE (add i1 i2) = fst (abstract_add precision rmode (toIEEE i1) (toIEEE i2)).
  Proof.
   tauto.
  Qed.

  Definition sub x y := fst (abstract_sub precision rmode x y).
  Lemma sub_prop : i1 i2,
    toIEEE (sub i1 i2) = fst (abstract_sub precision rmode (toIEEE i1) (toIEEE i2)).
  Proof.
   tauto.
  Qed.

  Definition mul x y := fst (abstract_mult precision rmode x y).
  Lemma mul_prop : i1 i2,
    toIEEE (mul i1 i2)= fst (abstract_mult precision rmode (toIEEE i1) (toIEEE i2)).
  Proof.
   tauto.
  Qed.

  Definition div x y := fst (abstract_div precision rmode x y).
  Lemma div_prop : i1 i2,
    toIEEE (div i1 i2)= fst (abstract_div precision rmode (toIEEE i1) (toIEEE i2)).
  Proof.
   tauto.
  Qed.

  Definition neg x := abstract_opp x.
  Lemma neg_prop : f,
    toIEEE (neg f) = abstract_opp (toIEEE f).
  Proof.
   tauto.
  Qed.

Function cmpg x y : Z :=
  match x, y with
   | NaN z, _ ⇒ 1%Z
   | Normal _ _ _, NaN z ⇒ 1%Z
   | Subnormal _ _, NaN z ⇒ 1%Z
   | Zero _, NaN z ⇒ 1%Z
   | Infty _, NaN z ⇒ 1%Z
   | _ , _match abstract_compare precision (toIEEE x) (toIEEE y) tso_default with
                    | (Ordered Datatypes.Gt, Nothing) ⇒ 1%Z
                    | (Ordered Datatypes.Eq, Nothing) ⇒ 0%Z
                    | (Ordered Datatypes.Lt, Nothing)-1%Z
                    | _-2%Z
                   end
  end.

  Lemma cmpg_gt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Gt, Nothing)cmpg f1 f2 = 1.
  Proof.
    intuition.
    induction f1, f2;
       try (unfold cmpg; try rewrite H1; auto).
Qed.

  Lemma cmpg_eq : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Eq, Nothing)cmpg f1 f2 = 0.
  Proof.
    intuition.
    induction f1, f2;
             (unfold cmpg; try rewrite H1; try discriminate; auto).
    assert False.
    apply H.
    unfold toIEEE.
    auto using isNaN.
    intuition.
Qed.

  Lemma cmpg_lt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Lt, Nothing)cmpg f1 f2 = -1.
  Proof.
    intuition.
    induction f1, f2;
       (unfold cmpg; try rewrite H1; try discriminate; auto).
  Qed.

  Lemma cmpg_nan : f1 f2,
    (isNaN (toIEEE f1)) (isNaN (toIEEE f2))cmpg f1 f2 = 1.
  Proof.
    unfold toIEEE.
    intuition.
    inversion H0; unfold cmpg; trivial.
    inversion H0; unfold cmpg.
    induction f1; trivial.
  Qed.

Function cmpl x y : Z :=
  match x, y with
   | NaN z, _-1%Z
   | Normal _ _ _, NaN z-1%Z
   | Subnormal _ _, NaN z-1%Z
   | Zero _, NaN z-1%Z
   | Infty _, NaN z-1%Z
   | _ , _match abstract_compare precision (toIEEE x) (toIEEE y) tso_default with
                    | (Ordered Datatypes.Gt, Nothing) ⇒ 1%Z
                    | (Ordered Datatypes.Eq, Nothing) ⇒ 0%Z
                    | (Ordered Datatypes.Lt, Nothing)-1%Z
                    | _-2%Z
                   end
  end.

  Lemma cmpl_gt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Gt, Nothing)cmpl f1 f2 = 1.
  Proof.
    intuition.
    induction f1, f2;
       try (unfold cmpl; try rewrite H1; try discriminate; auto).
  Qed.

  Lemma cmpl_eq : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Eq, Nothing)cmpl f1 f2 = 0.
  Proof.
    intuition.
    induction f1, f2;
       (unfold cmpl; try rewrite H1; try discriminate; auto).
    assert False.
    apply H.
    unfold toIEEE.
    auto using isNaN.
    intuition.
 Qed.

  Lemma cmpl_lt : f1 f2,
    (¬isNaN (toIEEE f1)) → (¬isNaN (toIEEE f2)) →
    abstract_compare precision (toIEEE f1) (toIEEE f2) tso_default = (Ordered Datatypes.Lt, Nothing)cmpl f1 f2 = -1.
  Proof.
    intuition.
    induction f1, f2;
       (unfold cmpl; try rewrite H1; try discriminate; auto).
  Qed.

  Lemma cmpl_nan : f1 f2,
    (isNaN (toIEEE f1)) (isNaN (toIEEE f2))cmpl f1 f2 = -1.
  Proof.
    unfold toIEEE.
    intuition.
    inversion H0; unfold cmpl; trivial.
    inversion H0; unfold cmpl.
    induction f1; trivial.
  Qed.

Function shiftr_iter (howmany : positive) (what : positive) :=
  match howmany return positive with
   | xHwhat
   | _match what with
               | xHwhat
               | xO xshiftr_iter (Pos.pred howmany) x
               | xI xshiftr_iter (Pos.pred howmany) x
              end
  end.

Defines the integer number the magnitude of which is as large as possible without exceeding the magnitude of the given number.
Definition round_mag_down_int (num : abstract_IEEE) :=
  match num return abstract_IEEE with
  | Normal sign e m
       match (e - dig_length precision) with
        | Z0Normal sign e m
this is an integer number
        | Zpos xNormal sign e m
this is an integer number
        | Zneg xNormal sign e (shiftr_iter x m)
x-times divide by 2 the value of m
       end
  | Subnormal b posZero b
  | Zero bZero b
  | Infty bInfty b
  | NaN posNaN pos
 end.

Definition abstract_rem (m : rounding_mode) (x y : abstract_IEEE) :=
  match x, y return (abstract_IEEE × exception) with
  | NaN s1, _(NaN s1, Nothing)
  | Normal s1 _ _, NaN _(NaN quiet, Nothing)
  | Subnormal s1 _, NaN _(NaN quiet, Nothing)
  | Zero s1, NaN _(NaN quiet, Nothing)
  | Infty s1, NaN _(NaN quiet, Nothing)

  | Infty b1, _(NaN quiet, Invalid_operation)
  | Normal b1 _ _ , Zero _(NaN quiet, Invalid_operation)
  | Subnormal b1 _ , Zero _(NaN quiet, Invalid_operation)
  | Zero b1, Zero _(NaN quiet, Invalid_operation)

  | Zero b1, Infty _(Zero b1, Nothing)
  | Normal b1 b2 b3, Infty _(Normal b1 b2 b3, Nothing)
  | Subnormal b1 b2, Infty _(Subnormal b1 b2, Nothing)

  | Zero b1, Normal _ _ _(Zero b1, Nothing)
  | Zero b1, Subnormal _ _(Zero b1, Nothing)

  | _, _
      let q := round_mag_down_int (fst (abstract_div precision m x y)) in
      let diffr := fst (abstract_mult precision m y q) in
          abstract_add precision m x (abstract_opp diffr)
  end.

  Definition rem x y := fst (abstract_rem rmode x y).
  Lemma rem_prop : i1 i2,
    toIEEE (rem i1 i2) = fst (abstract_rem rmode (toIEEE i1) (toIEEE i2)).
  Proof.
    auto.
  Qed.

  Definition zero := Zero true.
  Lemma zero_prop : toIEEE zero = Zero true.
  Proof.
    auto.
  Qed.

  Definition one := Normal true (dig_length precision) (Z.to_pos 1).
  Lemma one_prop : toIEEE one = Normal true (dig_length precision) (Z.to_pos 1).
  Proof.
    auto.
  Qed.

  Definition two := Normal true (1 + (dig_length precision)) (Z.to_pos 1).
  Lemma two_prop : toIEEE two = Normal true (1 + (dig_length precision)) (Z.to_pos 1).
  Proof.
    auto.
  Qed.

  Definition GtIEEE : abstract_IEEEabstract_IEEEProp :=
    fun x y(abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Gt, Nothing).

  Definition GeIEEE : abstract_IEEEabstract_IEEEProp :=
    fun x y(abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Gt, Nothing)
      (abstract_compare precision x y tso_default) =
                  (Ordered Datatypes.Eq, Nothing).

This function supplies maxint or minint in case the number num exceeds the range defined by the half base hb.
   Definition fit_into (num : Z) (hb : Z) :=
     num.

   Definition sign_of (sign : bool) :=
     if sign then 1 else -1.

The function defines the operation of floating point truncation to an integer. Its properties according to the Java Virtual Machine Specification are mentioned along the paramters below.
The first parameter is the floating point value to truncate, the second parameter is the number of bits of the resulting integer minus 1 (it corresponds do half_base in Numeric.

  Definition truncate (num : abstract_IEEE) (hb : Z):=
  match num return Z with
  | Normal sign e m
       match (e - dig_length precision) with
        | Z0fit_into
           ((sign_of sign × Zpos m)* (2^(e - dig_length precision))) hb
                    
Normal sign e m is an integer number
        | Zpos xfit_into
            ((sign_of sign × Zpos m)* (2^(e - dig_length precision))) hb
                    
Normal sign e m is an integer number
        | Zneg xfit_into
           ((sign_of sign × Zpos (shiftr_iter x m))*
                           (2^(e - dig_length precision))) hb
                    
x-times divide by 2 the value of m to get integer
       end
  | Subnormal b pos ⇒ 0
  | Zero b ⇒ 0
  | Infty true ⇒ 2^hb -1
  | Infty false-2^hb
  | NaN pos ⇒ 0
 end.

If the value' is NaN, the result of the conversion is 0.
   Lemma truncate_nan :
    f1 n, ( x, toIEEE f1 = NaN x) → truncate f1 n = 0.
   Proof.
   intros.
   unfold truncate.
   unfold toIEEE.
   unfold toIEEE in H.
   destruct H.
   rewrite H.
   intuition.
   Qed.

Otherwise, if the value' is not an infinity, it is rounded to an integer value V, rounding towards zero using IEEE 754 round towards zero mode. If this integer value V can be represented as an int, then the result is the int value V.
  Lemma truncate_range_positive :
    f half_base,
     let f1 := toIEEE f in
     let maxint := abstract_of_diadic precision Rounding_sup
                    (Diadic (2^half_base) 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
     let (truncp1, _) :=
           (abstract_add precision Rounding_sup trunc (toIEEE one)) in
       ( (GtIEEE maxint f1) (GeIEEE f1 (Zero true))) →
         ( (GtIEEE truncp1 f1) (GeIEEE f1 trunc)).
  Proof.
  admit.
  Qed.

  Lemma truncate_range_negative :
    f half_base,
     let f1 := toIEEE f in
     let minint := abstract_of_diadic precision Rounding_sup
                    (Diadic (2^half_base) 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
     let (truncm1, _) :=
           (abstract_sub precision Rounding_inf trunc (toIEEE one)) in
       ( (GeIEEE f1 minint) (GeIEEE (Zero true) f1)) →
         ( (GtIEEE f1 truncm1) (GeIEEE trunc f1)).
  Proof.
  admit.
  Qed.

Otherwise, either the value' must be too small (a negative value of large magnitude or negative infinity), and the result is the smallest representable value of type integer,
  Lemma truncate_too_small :
    f half_base,
     let f1 := toIEEE f in
     let minintZ := - 2^half_base in
     let minint := abstract_of_diadic precision Rounding_sup
                     (Diadic minintZ 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
       (GtIEEE minint f1 (f1 = Infty false)) →
        truncate f half_base = minintZ.
  Proof.
  intros.
  unfold toIEEE in ×.
  destruct H.
  admit.
  unfold truncate.
  unfold f1 in ×.
  rewrite H.
  auto.
  Qed.

...or the value' must be too large (a positive value of large magnitude or positive infinity), and the result is the largest representable value of type integer.
  Lemma truncate_too_large :
    f half_base,
     let f1 := toIEEE f in
     let maxintZ := 2^half_base - 1 in
     let maxint := abstract_of_diadic precision Rounding_sup
                     (Diadic maxintZ 0) in
     let trunc := abstract_of_diadic precision Rounding_zero
                    (Diadic (truncate f half_base) 0) in
       (GeIEEE f1 maxint (f1 = Infty true)) →
        truncate f half_base = maxintZ.
  Proof.
  intros.
  unfold toIEEE in ×.
  destruct H.
  admit.
  unfold f1 in ×.
  rewrite H.
  unfold truncate.
  trivial.
  Qed.

  Definition const (num: abstract_IEEE) := num.

  Definition iconst (num: Z) :=
    abstract_of_diadic precision rmode (Diadic num 0).

End FMake.