Library Flocq.IEEE754.Binary

This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2010-2018 Sylvie Boldo
Copyright (C) 2010-2018 Guillaume Melquiond
This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the COPYING file for more details.

IEEE-754 arithmetic

Require Import Core Digits Round Bracket Operations Div Sqrt Relative.

Section AnyRadix.

Inductive full_float :=
  | F754_zero (s : bool)
  | F754_infinity (s : bool)
  | F754_nan (s : bool) (m : positive)
  | F754_finite (s : bool) (m : positive) (e : Z).

Definition FF2R beta x :=
  match x with
  | F754_finite s m eF2R (Float beta (cond_Zopp s (Zpos m)) e)
  | _ ⇒ 0%R
  end.

End AnyRadix.

Section Binary.


prec is the number of bits of the mantissa including the implicit one; emax is the exponent of the infinities. For instance, binary32 is defined by prec = 24 and emax = 128.
Variable prec emax : Z.
Context (prec_gt_0_ : Prec_gt_0 prec).
Hypothesis Hmax : (prec < emax)%Z.

Let emin := (3 - emax - prec)%Z.
Let fexp := FLT_exp emin prec.
Instance fexp_correct : Valid_exp fexp := FLT_exp_valid emin prec.
Instance fexp_monotone : Monotone_exp fexp := FLT_exp_monotone emin prec.

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 nan_pl pl :=
  Zlt_bool (Zpos (digits2_pos pl)) prec.

Definition valid_binary x :=
  match x with
  | F754_finite _ m ebounded m e
  | F754_nan _ plnan_pl pl
  | _true
  end.

Basic type used for representing binary FP numbers. Note that there is exactly one such object per FP datum.

Inductive binary_float :=
  | B754_zero (s : bool)
  | B754_infinity (s : bool)
  | B754_nan (s : bool) (pl : positive) :
    nan_pl pl = true binary_float
  | B754_finite (s : bool) (m : positive) (e : Z) :
    bounded m e = true binary_float.

Definition FF2B x :=
  match x as x return valid_binary x = true binary_float with
  | F754_finite s m eB754_finite s m e
  | F754_infinity sfun _B754_infinity s
  | F754_zero sfun _B754_zero s
  | F754_nan b plfun HB754_nan b pl H
  end.

Definition B2FF x :=
  match x with
  | B754_finite s m e _F754_finite s m e
  | B754_infinity sF754_infinity s
  | B754_zero sF754_zero s
  | B754_nan b pl _F754_nan b pl
  end.

Definition B2R f :=
  match f with
  | B754_finite s m e _F2R (Float radix2 (cond_Zopp s (Zpos m)) e)
  | _ ⇒ 0%R
  end.

Theorem FF2R_B2FF :
   x,
  FF2R radix2 (B2FF x) = B2R x.

Theorem B2FF_FF2B :
   x Hx,
  B2FF (FF2B x Hx) = x.

Theorem valid_binary_B2FF :
   x,
  valid_binary (B2FF x) = true.

Theorem FF2B_B2FF :
   x H,
  FF2B (B2FF x) H = x.

Theorem FF2B_B2FF_valid :
   x,
  FF2B (B2FF x) (valid_binary_B2FF x) = x.

Theorem B2R_FF2B :
   x Hx,
  B2R (FF2B x Hx) = FF2R radix2 x.

Theorem match_FF2B :
   {T} fz fi fn ff x Hx,
  match FF2B x Hx return T with
  | B754_zero sxfz sx
  | B754_infinity sxfi sx
  | B754_nan b p _fn b p
  | B754_finite sx mx ex _ff sx mx ex
  end =
  match x with
  | F754_zero sxfz sx
  | F754_infinity sxfi sx
  | F754_nan b pfn b p
  | F754_finite sx mx exff sx mx ex
  end.

Theorem canonical_canonical_mantissa :
   (sx : bool) mx ex,
  canonical_mantissa mx ex = true
  canonical radix2 fexp (Float radix2 (cond_Zopp sx (Zpos mx)) ex).

Theorem generic_format_B2R :
   x,
  generic_format radix2 fexp (B2R x).

Theorem FLT_format_B2R :
   x,
  FLT_format radix2 emin prec (B2R x).

Theorem B2FF_inj :
   x y : binary_float,
  B2FF x = B2FF y
  x = y.

Definition is_finite_strict f :=
  match f with
  | B754_finite _ _ _ _true
  | _false
  end.

Theorem B2R_inj:
   x y : binary_float,
  is_finite_strict x = true
  is_finite_strict y = true
  B2R x = B2R y
  x = y.

Definition Bsign x :=
  match x with
  | B754_nan s _ _s
  | B754_zero ss
  | B754_infinity ss
  | B754_finite s _ _ _s
  end.

Definition sign_FF x :=
  match x with
  | F754_nan s _s
  | F754_zero ss
  | F754_infinity ss
  | F754_finite s _ _s
  end.

Theorem Bsign_FF2B :
   x H,
  Bsign (FF2B x H) = sign_FF x.

Definition is_finite f :=
  match f with
  | B754_finite _ _ _ _true
  | B754_zero _true
  | _false
  end.

Definition is_finite_FF f :=
  match f with
  | F754_finite _ _ _true
  | F754_zero _true
  | _false
  end.

Theorem is_finite_FF2B :
   x Hx,
  is_finite (FF2B x Hx) = is_finite_FF x.

Theorem is_finite_FF_B2FF :
   x,
  is_finite_FF (B2FF x) = is_finite x.

Theorem B2R_Bsign_inj:
   x y : binary_float,
    is_finite x = true
    is_finite y = true
    B2R x = B2R y
    Bsign x = Bsign y
    x = y.

Definition is_nan f :=
  match f with
  | B754_nan _ _ _true
  | _false
  end.

Definition is_nan_FF f :=
  match f with
  | F754_nan _ _true
  | _false
  end.

Theorem is_nan_FF2B :
   x Hx,
  is_nan (FF2B x Hx) = is_nan_FF x.

Theorem is_nan_FF_B2FF :
   x,
  is_nan_FF (B2FF x) = is_nan x.

Definition build_nan (x : { x | is_nan x = true }) :=
  proj1_sig x.

Theorem B2R_build_nan :
   x, B2R (build_nan x) = 0%R.

Theorem is_finite_build_nan :
   x, is_finite (build_nan x) = false.

Theorem is_nan_build_nan :
   x, is_nan (build_nan x) = true.

Opposite

Definition Bopp opp_nan x :=
  match x with
  | B754_nan _ _ _build_nan (opp_nan x)
  | B754_infinity sxB754_infinity (negb sx)
  | B754_finite sx mx ex HxB754_finite (negb sx) mx ex Hx
  | B754_zero sxB754_zero (negb sx)
  end.

Theorem Bopp_involutive :
   opp_nan x,
  is_nan x = false
  Bopp opp_nan (Bopp opp_nan x) = x.

Theorem B2R_Bopp :
   opp_nan x,
  B2R (Bopp opp_nan x) = (- B2R x)%R.

Theorem is_finite_Bopp :
   opp_nan x,
  is_finite (Bopp opp_nan x) = is_finite x.

Absolute value

Definition Babs abs_nan (x : binary_float) : binary_float :=
  match x with
  | B754_nan _ _ _build_nan (abs_nan x)
  | B754_infinity sxB754_infinity false
  | B754_finite sx mx ex HxB754_finite false mx ex Hx
  | B754_zero sxB754_zero false
  end.

Theorem B2R_Babs :
   abs_nan x,
  B2R (Babs abs_nan x) = Rabs (B2R x).

Theorem is_finite_Babs :
   abs_nan x,
  is_finite (Babs abs_nan x) = is_finite x.

Theorem Bsign_Babs :
   abs_nan x,
  is_nan x = false
  Bsign (Babs abs_nan x) = false.

Theorem Babs_idempotent :
   abs_nan (x: binary_float),
  is_nan x = false
  Babs abs_nan (Babs abs_nan x) = Babs abs_nan x.

Theorem Babs_Bopp :
   abs_nan opp_nan x,
  is_nan x = false
  Babs abs_nan (Bopp opp_nan x) = Babs abs_nan x.

Comparison
Some c means ordered as per c; None means unordered.

Definition Bcompare (f1 f2 : binary_float) : option comparison :=
  match f1, f2 with
  | B754_nan _ _ _,_ | _,B754_nan _ _ _None
  | B754_infinity true, B754_infinity true
  | B754_infinity false, B754_infinity falseSome Eq
  | B754_infinity true, _Some Lt
  | B754_infinity false, _Some Gt
  | _, B754_infinity trueSome Gt
  | _, B754_infinity falseSome Lt
  | B754_finite true _ _ _, B754_zero _Some Lt
  | B754_finite false _ _ _, B754_zero _Some Gt
  | B754_zero _, B754_finite true _ _ _Some Gt
  | B754_zero _, B754_finite false _ _ _Some Lt
  | B754_zero _, B754_zero _Some Eq
  | B754_finite s1 m1 e1 _, B754_finite s2 m2 e2 _
    match s1, s2 with
    | true, falseSome Lt
    | false, trueSome Gt
    | false, false
      match Zcompare e1 e2 with
      | LtSome Lt
      | GtSome Gt
      | EqSome (Pcompare m1 m2 Eq)
      end
    | true, true
      match Zcompare e1 e2 with
      | LtSome Gt
      | GtSome Lt
      | EqSome (CompOpp (Pcompare m1 m2 Eq))
      end
    end
  end.

Theorem Bcompare_correct :
   f1 f2,
  is_finite f1 = true is_finite f2 = true
  Bcompare f1 f2 = Some (Rcompare (B2R f1) (B2R f2)).
  Ltac apply_Rcompare :=
    match goal with
      | [ |- Some Lt = Some (Rcompare _ _) ] ⇒ f_equal; symmetry; apply Rcompare_Lt
      | [ |- Some Eq = Some (Rcompare _ _) ] ⇒ f_equal; symmetry; apply Rcompare_Eq
      | [ |- Some Gt = Some (Rcompare _ _) ] ⇒ f_equal; symmetry; apply Rcompare_Gt
    end.

Theorem Bcompare_swap :
   x y,
  Bcompare y x = match Bcompare x y with Some cSome (CompOpp c) | NoneNone end.

Theorem bounded_lt_emax :
   mx ex,
  bounded mx ex = true
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R.

Theorem bounded_ge_emin :
   mx ex,
  bounded mx ex = true
  (bpow radix2 emin F2R (Float radix2 (Zpos mx) ex))%R.

Theorem abs_B2R_lt_emax :
   x,
  (Rabs (B2R x) < bpow radix2 emax)%R.

Theorem bounded_canonical_lt_emax :
   mx ex,
  canonical radix2 fexp (Float radix2 (Zpos mx) ex)
  (F2R (Float radix2 (Zpos mx) ex) < bpow radix2 emax)%R
  bounded mx ex = true.

Truncation

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
  | Z0Build_shr_record Z0 false s
  | Zpos xHBuild_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 xHBuild_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 falseloc_Exact
  | Build_shr_record _ false trueloc_Inexact Lt
  | Build_shr_record _ true falseloc_Inexact Eq
  | Build_shr_record _ true trueloc_Inexact Gt
  end.

Definition shr_record_of_loc m l :=
  match l with
  | loc_ExactBuild_shr_record m false false
  | loc_Inexact LtBuild_shr_record m false true
  | loc_Inexact EqBuild_shr_record m true false
  | loc_Inexact GtBuild_shr_record m true true
  end.

Theorem shr_m_shr_record_of_loc :
   m l,
  shr_m (shr_record_of_loc m l) = m.

Theorem loc_of_shr_record_of_loc :
   m l,
  loc_of_shr_record (shr_record_of_loc m l) = l.

Definition shr mrs e n :=
  match n with
  | Zpos p(iter_pos shr_1 p mrs, (e + n)%Z)
  | _(mrs, e)
  end.

Lemma inbetween_shr_1 :
   x mrs e,
  (0 shr_m mrs)%Z
  inbetween_float radix2 (shr_m mrs) e x (loc_of_shr_record mrs)
  inbetween_float radix2 (shr_m (shr_1 mrs)) (e + 1) x (loc_of_shr_record (shr_1 mrs)).

Theorem inbetween_shr :
   x m e l n,
  (0 m)%Z
  inbetween_float radix2 m e x l
  let '(mrs, e') := shr (shr_record_of_loc m l) e n in
  inbetween_float radix2 (shr_m mrs) e' x (loc_of_shr_record mrs).

Definition shr_fexp m e l :=
  shr (shr_record_of_loc m l) e (fexp (Zdigits2 m + e) - e).

Theorem shr_truncate :
   m e l,
  (0 m)%Z
  shr_fexp m e l =
  let '(m', e', l') := truncate radix2 fexp (m, e, l) in (shr_record_of_loc m' l', e').

Rounding modes

Inductive mode := mode_NE | mode_ZR | mode_DN | mode_UP | mode_NA.

Definition round_mode m :=
  match m with
  | mode_NEZnearestE
  | mode_ZRZtrunc
  | mode_DNZfloor
  | mode_UPZceil
  | mode_NAZnearestA
  end.

Definition choice_mode m sx mx lx :=
  match m with
  | mode_NEcond_incr (round_N (negb (Z.even mx)) lx) mx
  | mode_ZRmx
  | mode_DNcond_incr (round_sign_DN sx lx) mx
  | mode_UPcond_incr (round_sign_UP sx lx) mx
  | mode_NAcond_incr (round_N true lx) mx
  end.

Global Instance valid_rnd_round_mode : m, Valid_rnd (round_mode m).

Definition overflow_to_inf m s :=
  match m with
  | mode_NEtrue
  | mode_NAtrue
  | mode_ZRfalse
  | mode_UPnegb s
  | mode_DNs
  end.

Definition binary_overflow m s :=
  if overflow_to_inf m s then F754_infinity s
  else F754_finite s (match (Zpower 2 prec - 1)%Z with Zpos pp | _xH end) (emax - prec).

Definition binary_round_aux mode sx mx ex lx :=
  let '(mrs', e') := shr_fexp mx ex lx in
  let '(mrs'', e'') := shr_fexp (choice_mode mode sx (shr_m mrs') (loc_of_shr_record mrs')) e' loc_Exact in
  match shr_m mrs'' with
  | Z0F754_zero sx
  | Zpos mif Zle_bool e'' (emax - prec) then F754_finite sx m e'' else binary_overflow mode sx
  | _F754_nan false xH
  end.

Theorem binary_round_aux_correct' :
   mode x mx ex lx,
  (x 0)%R
  inbetween_float radix2 mx ex (Rabs x) lx
  (ex cexp radix2 fexp x)%Z
  let z := binary_round_aux mode (Rlt_bool x 0) mx ex lx in
  valid_binary z = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode mode) x)) (bpow radix2 emax) then
    FF2R radix2 z = round radix2 fexp (round_mode mode) x
    is_finite_FF z = true sign_FF z = Rlt_bool x 0
  else
    z = binary_overflow mode (Rlt_bool x 0).

Theorem binary_round_aux_correct :
   mode x mx ex lx,
  inbetween_float radix2 (Zpos mx) ex (Rabs x) lx
  (ex fexp (Zdigits radix2 (Zpos mx) + ex))%Z
  let z := binary_round_aux mode (Rlt_bool x 0) (Zpos mx) ex lx in
  valid_binary z = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode mode) x)) (bpow radix2 emax) then
    FF2R radix2 z = round radix2 fexp (round_mode mode) x
    is_finite_FF z = true sign_FF z = Rlt_bool x 0
  else
    z = binary_overflow mode (Rlt_bool x 0).

Multiplication

Lemma Bmult_correct_aux :
   m sx mx ex (Hx : bounded mx ex = true) sy my ey (Hy : bounded my ey = true),
  let x := F2R (Float radix2 (cond_Zopp sx (Zpos mx)) ex) in
  let y := F2R (Float radix2 (cond_Zopp sy (Zpos my)) ey) in
  let z := binary_round_aux m (xorb sx sy) (Zpos (mx × my)) (ex + ey) loc_Exact in
  valid_binary z = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x × y))) (bpow radix2 emax) then
    FF2R radix2 z = round radix2 fexp (round_mode m) (x × y)
    is_finite_FF z = true sign_FF z = xorb sx sy
  else
    z = binary_overflow m (xorb sx sy).

Definition Bmult mult_nan m x y :=
  match x, y with
  | B754_nan _ _ _, _ | _, B754_nan _ _ _build_nan (mult_nan x y)
  | B754_infinity sx, B754_infinity syB754_infinity (xorb sx sy)
  | B754_infinity sx, B754_finite sy _ _ _B754_infinity (xorb sx sy)
  | B754_finite sx _ _ _, B754_infinity syB754_infinity (xorb sx sy)
  | B754_infinity _, B754_zero _build_nan (mult_nan x y)
  | B754_zero _, B754_infinity _build_nan (mult_nan x y)
  | B754_finite sx _ _ _, B754_zero syB754_zero (xorb sx sy)
  | B754_zero sx, B754_finite sy _ _ _B754_zero (xorb sx sy)
  | B754_zero sx, B754_zero syB754_zero (xorb sx sy)
  | B754_finite sx mx ex Hx, B754_finite sy my ey Hy
    FF2B _ (proj1 (Bmult_correct_aux m sx mx ex Hx sy my ey Hy))
  end.

Theorem Bmult_correct :
   mult_nan m x y,
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x × B2R y))) (bpow radix2 emax) then
    B2R (Bmult mult_nan m x y) = round radix2 fexp (round_mode m) (B2R x × B2R y)
    is_finite (Bmult mult_nan m x y) = andb (is_finite x) (is_finite y)
    (is_nan (Bmult mult_nan m x y) = false
      Bsign (Bmult mult_nan m x y) = xorb (Bsign x) (Bsign y))
  else
    B2FF (Bmult mult_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)).

Normalization and rounding

Definition shl_align mx ex ex' :=
  match (ex' - ex)%Z with
  | Zneg d(shift_pos d mx, ex')
  | _(mx, ex)
  end.

Theorem shl_align_correct :
   mx ex ex',
  let (mx', ex'') := shl_align mx ex ex' in
  F2R (Float radix2 (Zpos mx) ex) = F2R (Float radix2 (Zpos mx') ex'')
  (ex'' ex')%Z.

Theorem snd_shl_align :
   mx ex ex',
  (ex' ex)%Z
  snd (shl_align mx ex ex') = ex'.

Definition shl_align_fexp mx ex :=
  shl_align mx ex (fexp (Zpos (digits2_pos mx) + ex)).

Theorem shl_align_fexp_correct :
   mx ex,
  let (mx', ex') := shl_align_fexp mx ex in
  F2R (Float radix2 (Zpos mx) ex) = F2R (Float radix2 (Zpos mx') ex')
  (ex' fexp (Zdigits radix2 (Zpos mx') + ex'))%Z.

Definition binary_round m sx mx ex :=
  let '(mz, ez) := shl_align_fexp mx ex in binary_round_aux m sx (Zpos mz) ez loc_Exact.

Theorem binary_round_correct :
   m sx mx ex,
  let z := binary_round m sx mx ex in
  valid_binary z = true
  let x := F2R (Float radix2 (cond_Zopp sx (Zpos mx)) ex) in
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) x)) (bpow radix2 emax) then
    FF2R radix2 z = round radix2 fexp (round_mode m) x
    is_finite_FF z = true
    sign_FF z = sx
  else
    z = binary_overflow m sx.

Definition binary_normalize mode m e szero :=
  match m with
  | Z0B754_zero szero
  | Zpos mFF2B _ (proj1 (binary_round_correct mode false m e))
  | Zneg mFF2B _ (proj1 (binary_round_correct mode true m e))
  end.

Theorem binary_normalize_correct :
   m mx ex szero,
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (F2R (Float radix2 mx ex)))) (bpow radix2 emax) then
    B2R (binary_normalize m mx ex szero) = round radix2 fexp (round_mode m) (F2R (Float radix2 mx ex))
    is_finite (binary_normalize m mx ex szero) = true
    Bsign (binary_normalize m mx ex szero) =
      match Rcompare (F2R (Float radix2 mx ex)) 0 with
        | Eqszero
        | Lttrue
        | Gtfalse
      end
  else
    B2FF (binary_normalize m mx ex szero) = binary_overflow m (Rlt_bool (F2R (Float radix2 mx ex)) 0).

Addition

Definition Bplus plus_nan m x y :=
  match x, y with
  | B754_nan _ _ _, _ | _, B754_nan _ _ _build_nan (plus_nan x y)
  | B754_infinity sx, B754_infinity sy
    if Bool.eqb sx sy then x else build_nan (plus_nan x y)
  | B754_infinity _, _x
  | _, B754_infinity _y
  | B754_zero sx, B754_zero sy
    if Bool.eqb sx sy then x else
    match m with mode_DNB754_zero true | _B754_zero false end
  | B754_zero _, _y
  | _, B754_zero _x
  | B754_finite sx mx ex Hx, B754_finite sy my ey Hy
    let ez := Zmin ex ey in
    binary_normalize m (Zplus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez)))))
      ez (match m with mode_DNtrue | _false end)
  end.

Theorem Bplus_correct :
   plus_nan m x y,
  is_finite x = true
  is_finite y = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x + B2R y))) (bpow radix2 emax) then
    B2R (Bplus plus_nan m x y) = round radix2 fexp (round_mode m) (B2R x + B2R y)
    is_finite (Bplus plus_nan m x y) = true
    Bsign (Bplus plus_nan m x y) =
      match Rcompare (B2R x + B2R y) 0 with
        | Eqmatch m with mode_DNorb (Bsign x) (Bsign y)
                                 | _andb (Bsign x) (Bsign y) end
        | Lttrue
        | Gtfalse
      end
  else
    (B2FF (Bplus plus_nan m x y) = binary_overflow m (Bsign x) Bsign x = Bsign y).

Subtraction

Definition Bminus minus_nan m x y :=
  match x, y with
  | B754_nan _ _ _, _ | _, B754_nan _ _ _build_nan (minus_nan x y)
  | B754_infinity sx, B754_infinity sy
    if Bool.eqb sx (negb sy) then x else build_nan (minus_nan x y)
  | B754_infinity _, _x
  | _, B754_infinity syB754_infinity (negb sy)
  | B754_zero sx, B754_zero sy
    if Bool.eqb sx (negb sy) then x else
    match m with mode_DNB754_zero true | _B754_zero false end
  | B754_zero _, B754_finite sy my ey HyB754_finite (negb sy) my ey Hy
  | _, B754_zero _x
  | B754_finite sx mx ex Hx, B754_finite sy my ey Hy
    let ez := Zmin ex ey in
    binary_normalize m (Zminus (cond_Zopp sx (Zpos (fst (shl_align mx ex ez)))) (cond_Zopp sy (Zpos (fst (shl_align my ey ez)))))
      ez (match m with mode_DNtrue | _false end)
  end.

Theorem Bminus_correct :
   minus_nan m x y,
  is_finite x = true
  is_finite y = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x - B2R y))) (bpow radix2 emax) then
    B2R (Bminus minus_nan m x y) = round radix2 fexp (round_mode m) (B2R x - B2R y)
    is_finite (Bminus minus_nan m x y) = true
    Bsign (Bminus minus_nan m x y) =
      match Rcompare (B2R x - B2R y) 0 with
        | Eqmatch m with mode_DNorb (Bsign x) (negb (Bsign y))
                                 | _andb (Bsign x) (negb (Bsign y)) end
        | Lttrue
        | Gtfalse
      end
  else
    (B2FF (Bminus minus_nan m x y) = binary_overflow m (Bsign x) Bsign x = negb (Bsign y)).

Division

Definition Fdiv_core_binary m1 e1 m2 e2 :=
  let d1 := Zdigits2 m1 in
  let d2 := Zdigits2 m2 in
  let e' := Zmin (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
    | Z0m1
    | Zneg _Z0
    end in
  let '(q, r) := Zfast_div_eucl m' m2 in
  (q, e', new_location m2 r loc_Exact).

Lemma Bdiv_correct_aux :
   m sx mx ex sy my ey,
  let x := F2R (Float radix2 (cond_Zopp sx (Zpos mx)) ex) in
  let y := F2R (Float radix2 (cond_Zopp sy (Zpos my)) ey) in
  let z :=
    let '(mz, ez, lz) := Fdiv_core_binary (Zpos mx) ex (Zpos my) ey in
    binary_round_aux m (xorb sx sy) mz ez lz in
  valid_binary z = true
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (x / y))) (bpow radix2 emax) then
    FF2R radix2 z = round radix2 fexp (round_mode m) (x / y)
    is_finite_FF z = true sign_FF z = xorb sx sy
  else
    z = binary_overflow m (xorb sx sy).

Definition Bdiv div_nan m x y :=
  match x, y with
  | B754_nan _ _ _, _ | _, B754_nan _ _ _build_nan (div_nan x y)
  | B754_infinity sx, B754_infinity sybuild_nan (div_nan x y)
  | B754_infinity sx, B754_finite sy _ _ _B754_infinity (xorb sx sy)
  | B754_finite sx _ _ _, B754_infinity syB754_zero (xorb sx sy)
  | B754_infinity sx, B754_zero syB754_infinity (xorb sx sy)
  | B754_zero sx, B754_infinity syB754_zero (xorb sx sy)
  | B754_finite sx _ _ _, B754_zero syB754_infinity (xorb sx sy)
  | B754_zero sx, B754_finite sy _ _ _B754_zero (xorb sx sy)
  | B754_zero sx, B754_zero sybuild_nan (div_nan x y)
  | B754_finite sx mx ex _, B754_finite sy my ey _
    FF2B _ (proj1 (Bdiv_correct_aux m sx mx ex sy my ey))
  end.

Theorem Bdiv_correct :
   div_nan m x y,
  B2R y 0%R
  if Rlt_bool (Rabs (round radix2 fexp (round_mode m) (B2R x / B2R y))) (bpow radix2 emax) then
    B2R (Bdiv div_nan m x y) = round radix2 fexp (round_mode m) (B2R x / B2R y)
    is_finite (Bdiv div_nan m x y) = is_finite x
    (is_nan (Bdiv div_nan m x y) = false
      Bsign (Bdiv div_nan m x y) = xorb (Bsign x) (Bsign y))
  else
    B2FF (Bdiv div_nan m x y) = binary_overflow m (xorb (Bsign x) (Bsign y)).

Square root

Definition Fsqrt_core_binary m e :=
  let d := Zdigits2 m in
  let e' := Zmin (fexp (Z.div2 (d + e + 1))) (Z.div2 e) in
  let s := (e - 2 × e')%Z in
  let m' :=
    match s with
    | Zpos pZ.shiftl m s
    | Z0m
    | 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).

Lemma Bsqrt_correct_aux :
   m mx ex (Hx : bounded mx ex = true),
  let x := F2R (Float radix2 (Zpos mx) ex) in
  let z :=
    let '(mz, ez, lz) := Fsqrt_core_binary (Zpos mx) ex in
    binary_round_aux m false mz ez lz in
  valid_binary z = true
  FF2R radix2 z = round radix2 fexp (round_mode m) (sqrt x)
  is_finite_FF z = true sign_FF z = false.

Definition Bsqrt sqrt_nan m x :=
  match x with
  | B754_nan sx plx _build_nan (sqrt_nan x)
  | B754_infinity falsex
  | B754_infinity truebuild_nan (sqrt_nan x)
  | B754_finite true _ _ _build_nan (sqrt_nan x)
  | B754_zero _x
  | B754_finite sx mx ex Hx
    FF2B _ (proj1 (Bsqrt_correct_aux m mx ex Hx))
  end.

Theorem Bsqrt_correct :
   sqrt_nan m x,
  B2R (Bsqrt sqrt_nan m x) = round radix2 fexp (round_mode m) (sqrt (B2R x))
  is_finite (Bsqrt sqrt_nan m x) = match x with B754_zero _true | B754_finite false _ _ _true | _false end
  (is_nan (Bsqrt sqrt_nan m x) = false Bsign (Bsqrt sqrt_nan m x) = Bsign x).

End Binary.