diff options
| -rw-r--r-- | theories/Floats/FloatAxioms.v | 66 | ||||
| -rw-r--r-- | theories/Floats/Floats.v | 3 | ||||
| -rw-r--r-- | theories/Floats/PrimFloat.v | 95 | ||||
| -rw-r--r-- | theories/Floats/SpecFloat.v | 337 |
4 files changed, 501 insertions, 0 deletions
diff --git a/theories/Floats/FloatAxioms.v b/theories/Floats/FloatAxioms.v new file mode 100644 index 0000000000..4d74edddab --- /dev/null +++ b/theories/Floats/FloatAxioms.v @@ -0,0 +1,66 @@ +Require Import ZArith Int63 SpecFloat PrimFloat. + +(* Properties of the Binary64 IEEE 754 format *) +Definition prec := 53%Z. +Definition emax := 1024%Z. + +Notation valid_binary := (valid_binary prec emax). + +Definition SF64mul := SFmul prec emax. +Definition SF64add := SFadd prec emax. +Definition SF64sub := SFsub prec emax. +Definition SF64div := SFdiv prec emax. +Definition SF64sqrt := SFsqrt prec emax. + +Definition Prim2SF f := + if is_nan f then S754_nan + else if is_zero f then S754_zero (get_sign f) + else if is_infinity f then S754_infinity (get_sign f) + else + let (r, exp) := frexp f in + let e := (exp - prec)%Z in + let (shr, e') := shr_fexp prec emax [| normfr_mantissa r |]%int63 e loc_Exact in + match shr_m shr with + | Zpos p => S754_finite (get_sign f) p e' + | Zneg _ | Z0 => S754_zero false (* must never occur *) + end. + +Definition SF2Prim ef := + match ef with + | S754_nan => nan + | S754_zero false => zero + | S754_zero true => neg_zero + | S754_infinity false => infinity + | S754_infinity true => neg_infinity + | S754_finite s m e => + let pm := of_int63 (of_Z (Zpos m)) in + let f := ldexp pm e in + if s then (-f)%float else f + end. + +Axiom Prim2SF_valid : forall x, valid_binary (Prim2SF x) = true. +Axiom SF2Prim_Prim2SF : forall x, SF2Prim (Prim2SF x) = x. +Axiom Prim2SF_SF2Prim : forall x, valid_binary x = true -> Prim2SF (SF2Prim x) = x. + +Theorem Prim2SF_inj : forall x y, Prim2SF x = Prim2SF y -> x = y. + intros. rewrite <- SF2Prim_Prim2SF. symmetry. rewrite <- SF2Prim_Prim2SF. now rewrite H. +Qed. + +Theorem SF2Prim_inj : forall x y, SF2Prim x = SF2Prim y -> valid_binary x = true -> valid_binary y = true -> x = y. + intros. rewrite <- Prim2SF_SF2Prim by assumption. symmetry. rewrite <- Prim2SF_SF2Prim by assumption. rewrite H. reflexivity. +Qed. + +Axiom opp_spec : forall x, Prim2SF (-x)%float = SFopp (Prim2SF x). +Axiom abs_spec : forall x, Prim2SF (abs x) = SFabs (Prim2SF x). +Axiom compare_spec : forall x y, (x ?= y)%float = SFcompare (Prim2SF x) (Prim2SF y). +Axiom mul_spec : forall x y, Prim2SF (x * y)%float = SF64mul (Prim2SF x) (Prim2SF y). +Axiom add_spec : forall x y, Prim2SF (x + y)%float = SF64add (Prim2SF x) (Prim2SF y). +Axiom sub_spec : forall x y, Prim2SF (x - y)%float = SF64sub (Prim2SF x) (Prim2SF y). +Axiom div_spec : forall x y, Prim2SF (x / y)%float = SF64div (Prim2SF x) (Prim2SF y). +Axiom sqrt_spec : forall x, Prim2SF (sqrt x) = SF64sqrt (Prim2SF x). + +Axiom of_int63_spec : forall n, Prim2SF (of_int63 n) = binary_normalize prec emax (to_Z n) 0%Z false. +Axiom normfr_mantissa_spec : forall f, to_Z (normfr_mantissa f) = Z.of_N (SFnormfr_mantissa prec (Prim2SF f)). + +Axiom frexp_spec : forall f, let (m,e) := frexp f in (Prim2SF m, e) = SFfrexp prec emax (Prim2SF f). +Axiom ldexp_spec : forall f e, Prim2SF (ldexp f e) = SFldexp prec emax (Prim2SF f) e. diff --git a/theories/Floats/Floats.v b/theories/Floats/Floats.v new file mode 100644 index 0000000000..818de9ffb6 --- /dev/null +++ b/theories/Floats/Floats.v @@ -0,0 +1,3 @@ +Require Export SpecFloat. +Require Export PrimFloat. +Require Export FloatAxioms. diff --git a/theories/Floats/PrimFloat.v b/theories/Floats/PrimFloat.v new file mode 100644 index 0000000000..4cc0a9c4d5 --- /dev/null +++ b/theories/Floats/PrimFloat.v @@ -0,0 +1,95 @@ +Require Import ZArith. +Require Import Int63. + +Primitive float := #float64_type. + +Declare Scope float_scope. +Delimit Scope float_scope with float. +Bind Scope float_scope with float. + +Primitive opp := #float64_opp. +Notation "- x" := (opp x) : float_scope. + +Primitive abs := #float64_abs. + +Register option as kernel.ind_option. +Primitive compare := #float64_compare. +Notation "x ?= y" := (compare x y) (at level 70, no associativity) : float_scope. + +Primitive mul := #float64_mul. +Notation "x * y" := (mul x y) : float_scope. + +Primitive add := #float64_add. +Notation "x + y" := (add x y) : float_scope. +Primitive sub := #float64_sub. +Notation "x - y" := (sub x y) : float_scope. + +Primitive div := #float64_div. +Notation "x / y" := (div x y) : float_scope. + +Primitive sqrt := #float64_sqrt. + +(* Convert a primitive integer into a float value. + * The value is rounded if necessary. *) +Primitive of_int63 := #float64_of_int63. + +(* If the input is a float value with an absolute value + * inside [0.5; 1.) then return its mantissa as a primitive integer. + * The mantissa will be a 53-bit integer with its most significant bit set to 1. + * Else return zero. + * The sign bit is always ignored. *) +Primitive normfr_mantissa := #float64_normfr_mantissa. + +(* Exponent manipulation functions *) +Definition shift := (1022 + 52)%int63. +Primitive frshiftexp := #float64_frshiftexp. +Primitive ldshiftexp := #float64_ldshiftexp. + +Definition frexp f := + let (m, se) := frshiftexp f in + (m, ([| se |] - [| shift |])%Z%int63). + +Definition ldexp f e := ldshiftexp f (of_Z e + shift). + +Local Open Scope float_scope. + +(* Special values *) +Definition infinity := Eval compute in div (of_int63 1) (of_int63 0). +Definition neg_infinity := Eval compute in opp infinity. +Definition nan := Eval compute in div (of_int63 0) (of_int63 0). + +Definition one := Eval compute in (of_int63 1). +Definition zero := Eval compute in (of_int63 0). +Definition neg_zero := Eval compute in (-zero). +Definition two := Eval compute in (of_int63 2). + +Definition is_nan f := + match f ?= f with + | None => true + | _ => false + end. + +Definition is_zero f := + match f ?= zero with + | Some Eq => true (* note: 0 == -0 with floats *) + | _ => false + end. + +Definition is_infinity f := + match f ?= infinity with + | Some Eq => true + | Some Lt => match f ?= neg_infinity with + | Some Eq => true + | _ => false + end + | _ => false + end. + +Definition get_sign f := (* + => false, - => true *) + let f := if is_zero f then one / f else f in + match f ?= zero with + | Some Gt => false + | _ => true + end. + +Definition is_finite (x : float) := negb (is_nan x || is_infinity x). diff --git a/theories/Floats/SpecFloat.v b/theories/Floats/SpecFloat.v new file mode 100644 index 0000000000..60a70fc230 --- /dev/null +++ b/theories/Floats/SpecFloat.v @@ -0,0 +1,337 @@ +Require Import ZArith. + +Variant spec_float := + | S754_zero (s : bool) + | S754_infinity (s : bool) + | S754_nan + | S754_finite (s : bool) (m : positive) (e : Z). + +Section FloatOps. + Variable prec emax : Z. + + Definition emin := (3-emax-prec)%Z. + Definition fexp e := Z.max (e - prec) emin. + + Section Zdigits2. + Fixpoint digits2_pos (n : positive) : positive := + match n with + | xH => xH + | xO p => Pos.succ (digits2_pos p) + | xI p => Pos.succ (digits2_pos p) + end. + + Definition Zdigits2 n := + match n with + | Z0 => n + | Zpos p => Zpos (digits2_pos p) + | Zneg p => Zpos (digits2_pos p) + end. + End Zdigits2. + + Section ValidBinary. + Definition canonical_mantissa m e := + Zeq_bool (fexp (Zpos (digits2_pos m) + e)) e. + + Definition bounded m e := + andb (canonical_mantissa m e) (Zle_bool e (emax - prec)). + + Definition valid_binary x := + match x with + | S754_finite _ m e => bounded m e + | _ => true + end. + End ValidBinary. + + Section Iter. + Context {A : Type}. + Variable (f : A -> A). + + Fixpoint iter_pos (n : positive) (x : A) {struct n} : A := + match n with + | xI n' => iter_pos n' (iter_pos n' (f x)) + | xO n' => iter_pos n' (iter_pos n' x) + | xH => f x + end. + End Iter. + + Section Rounding. + Inductive location := loc_Exact | loc_Inexact : comparison -> location. + + Record shr_record := { shr_m : Z ; shr_r : bool ; shr_s : bool }. + + Definition shr_1 mrs := + let '(Build_shr_record m r s) := mrs in + let s := orb r s in + match m with + | Z0 => Build_shr_record Z0 false s + | Zpos xH => Build_shr_record Z0 true s + | Zpos (xO p) => Build_shr_record (Zpos p) false s + | Zpos (xI p) => Build_shr_record (Zpos p) true s + | Zneg xH => Build_shr_record Z0 true s + | Zneg (xO p) => Build_shr_record (Zneg p) false s + | Zneg (xI p) => Build_shr_record (Zneg p) true s + end. + + Definition loc_of_shr_record mrs := + match mrs with + | Build_shr_record _ false false => loc_Exact + | Build_shr_record _ false true => loc_Inexact Lt + | Build_shr_record _ true false => loc_Inexact Eq + | Build_shr_record _ true true => loc_Inexact Gt + end. + + Definition shr_record_of_loc m l := + match l with + | loc_Exact => Build_shr_record m false false + | loc_Inexact Lt => Build_shr_record m false true + | loc_Inexact Eq => Build_shr_record m true false + | loc_Inexact Gt => Build_shr_record m true true + end. + + Definition shr mrs e n := + match n with + | Zpos p => (iter_pos shr_1 p mrs, (e + n)%Z) + | _ => (mrs, e) + end. + + Definition shr_fexp m e l := + shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e). + + Definition round_nearest_even mx lx := + match lx with + | loc_Exact => mx + | loc_Inexact Lt => mx + | loc_Inexact Eq => if Z.even mx then mx else (mx + 1)%Z + | loc_Inexact Gt => (mx + 1)%Z + end. + + Definition binary_round_aux sx mx ex lx := + let '(mrs', e') := shr_fexp mx ex lx in + let '(mrs'', e'') := shr_fexp (round_nearest_even (shr_m mrs') (loc_of_shr_record mrs')) e' loc_Exact in + match shr_m mrs'' with + | Z0 => S754_zero sx + | Zpos m => if Zle_bool e'' (emax - prec) then S754_finite sx m e'' else S754_infinity sx + | _ => S754_nan + end. + + Definition shl_align mx ex ex' := + match (ex' - ex)%Z with + | Zneg d => (shift_pos d mx, ex') + | _ => (mx, ex) + end. + + Definition binary_round sx mx ex := + let '(mz, ez) := shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex))in + binary_round_aux sx (Zpos mz) ez loc_Exact. + + Definition binary_normalize m e szero := + match m with + | Z0 => S754_zero szero + | Zpos m => binary_round false m e + | Zneg m => binary_round true m e + end. + End Rounding. + + (* Define operations *) + Definition SFopp x := + match x with + | S754_nan => S754_nan + | S754_infinity sx => S754_infinity (negb sx) + | S754_finite sx mx ex => S754_finite (negb sx) mx ex + | S754_zero sx => S754_zero (negb sx) + end. + + Definition SFabs x := + match x with + | S754_nan => S754_nan + | S754_infinity sx => S754_infinity false + | S754_finite sx mx ex => S754_finite false mx ex + | S754_zero sx => S754_zero false + end. + + Definition SFcompare f1 f2 := + match f1, f2 with + | S754_nan , _ | _, S754_nan => None + | S754_infinity s1, S754_infinity s2 => + Some match s1, s2 with + | true, true => Eq + | false, false => Eq + | true, false => Lt + | false, true => Gt + end + | S754_infinity s, _ => Some (if s then Lt else Gt) + | _, S754_infinity s => Some (if s then Gt else Lt) + | S754_finite s _ _, S754_zero _ => Some (if s then Lt else Gt) + | S754_zero _, S754_finite s _ _ => Some (if s then Gt else Lt) + | S754_zero _, S754_zero _ => Some Eq + | S754_finite s1 m1 e1, S754_finite s2 m2 e2 => + Some match s1, s2 with + | true, false => Lt + | false, true => Gt + | false, false => + match Z.compare e1 e2 with + | Lt => Lt + | Gt => Gt + | Eq => Pcompare m1 m2 Eq + end + | true, true => + match Z.compare e1 e2 with + | Lt => Gt + | Gt => Lt + | Eq => CompOpp (Pcompare m1 m2 Eq) + end + end + end. + + Definition SFmul x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => S754_infinity (xorb sx sy) + | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy) + | S754_finite sx _ _, S754_infinity sy => S754_infinity (xorb sx sy) + | S754_infinity _, S754_zero _ => S754_nan + | S754_zero _, S754_infinity _ => S754_nan + | S754_finite sx _ _, S754_zero sy => S754_zero (xorb sx sy) + | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy) + | S754_zero sx, S754_zero sy => S754_zero (xorb sx sy) + | S754_finite sx mx ex, S754_finite sy my ey => + binary_round_aux (xorb sx sy) (Zpos (mx * my)) (ex + ey) loc_Exact + end. + + Definition cond_Zopp (b : bool) m := if b then Z.opp m else m. + + Definition SFadd x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => + if Bool.eqb sx sy then x else S754_nan + | S754_infinity _, _ => x + | _, S754_infinity _ => y + | S754_zero sx, S754_zero sy => + if Bool.eqb sx sy then x else + S754_zero false + | S754_zero _, _ => y + | _, S754_zero _ => x + | S754_finite sx mx ex, S754_finite sy my ey => + let ez := Z.min ex ey in + binary_normalize (Zplus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez))))) + ez false + end. + + Definition SFsub x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => + if Bool.eqb sx (negb sy) then x else S754_nan + | S754_infinity _, _ => x + | _, S754_infinity sy => S754_infinity (negb sy) + | S754_zero sx, S754_zero sy => + if Bool.eqb sx (negb sy) then x else + S754_zero false + | S754_zero _, S754_finite sy my ey => S754_finite (negb sy) my ey + | _, S754_zero _ => x + | S754_finite sx mx ex, S754_finite sy my ey => + let ez := Z.min ex ey in + binary_normalize (Zminus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez))))) + ez false + end. + + Definition new_location_even nb_steps k := + if Zeq_bool k 0 then loc_Exact + else loc_Inexact (Z.compare (2 * k) nb_steps). + + Definition new_location_odd nb_steps k := + if Zeq_bool k 0 then loc_Exact + else + loc_Inexact + match Z.compare (2 * k + 1) nb_steps with + | Lt => Lt + | Eq => Lt + | Gt => Gt + end. + + Definition new_location nb_steps := + if Z.even nb_steps then new_location_even nb_steps else new_location_odd nb_steps. + + Definition SFdiv_core_binary m1 e1 m2 e2 := + let d1 := Zdigits2 m1 in + let d2 := Zdigits2 m2 in + let e' := Z.min (fexp (d1 + e1 - (d2 + e2))) (e1 - e2) in + let s := (e1 - e2 - e')%Z in + let m' := + match s with + | Zpos _ => Z.shiftl m1 s + | Z0 => m1 + | Zneg _ => Z0 + end in + let '(q, r) := Z.div_eucl m' m2 in + (q, e', new_location m2 r). + + Definition SFdiv x y := + match x, y with + | S754_nan, _ | _, S754_nan => S754_nan + | S754_infinity sx, S754_infinity sy => S754_nan + | S754_infinity sx, S754_finite sy _ _ => S754_infinity (xorb sx sy) + | S754_finite sx _ _, S754_infinity sy => S754_zero (xorb sx sy) + | S754_infinity sx, S754_zero sy => S754_infinity (xorb sx sy) + | S754_zero sx, S754_infinity sy => S754_zero (xorb sx sy) + | S754_finite sx _ _, S754_zero sy => S754_infinity (xorb sx sy) + | S754_zero sx, S754_finite sy _ _ => S754_zero (xorb sx sy) + | S754_zero sx, S754_zero sy => S754_nan + | S754_finite sx mx ex, S754_finite sy my ey => + let '(mz, ez, lz) := SFdiv_core_binary (Zpos mx) ex (Zpos my) ey in + binary_round_aux (xorb sx sy) mz ez lz + end. + + Definition SFsqrt_core_binary m e := + let d := Zdigits2 m in + let e' := Z.min (fexp (Z.div2 (d + e + 1))) (Z.div2 e) in + let s := (e - 2 * e')%Z in + let m' := + match s with + | Zpos p => Z.shiftl m s + | Z0 => m + | Zneg _ => Z0 + end in + let (q, r) := Z.sqrtrem m' in + let l := + if Zeq_bool r 0 then loc_Exact + else loc_Inexact (if Zle_bool r q then Lt else Gt) in + (q, e', l). + + Definition SFsqrt x := + match x with + | S754_nan => S754_nan + | S754_infinity false => x + | S754_infinity true => S754_nan + | S754_finite true _ _ => S754_nan + | S754_zero _ => x + | S754_finite sx mx ex => + let '(mz, ez, lz) := SFsqrt_core_binary (Zpos mx) ex in + binary_round_aux false mz ez lz + end. + + Definition SFnormfr_mantissa f := + match f with + | S754_finite _ mx ex => + if Z.eqb ex (-prec) then Npos mx else 0%N + | _ => 0%N + end. + + Definition SFldexp f e := + match f with + | S754_finite sx mx ex => binary_round sx mx (ex+e) + | _ => f + end. + + Definition SFfrexp f := + match f with + | S754_finite sx mx ex => + if (Z.to_pos prec <=? digits2_pos mx)%positive then + (S754_finite sx mx (-prec), (ex+prec)%Z) + else + let d := (prec - Z.pos (digits2_pos mx))%Z in + (S754_finite sx (shift_pos (Z.to_pos d) mx) (-prec), (ex+prec-d)%Z) + | _ => (f, emin%Z) + end. +End FloatOps. |