6 files changed, 236 insertions, 26 deletions
diff --git a/flocq/Calc/Bracket.v b/flocq/Calc/Bracket.v
index 838cadfa..9ab55165 100644
--- a/flocq/Calc/Bracket.v
+++ b/flocq/Calc/Bracket.v
@@ -19,13 +19,14 @@ COPYING file for more details.
 
 (** * Locations: where a real number is positioned with respect to its rounded-down value in an arbitrary format. *)
 
-From Coq Require Import Lia.
-Require Import Raux Defs Float_prop.
-Require Import SpecFloatCompat.
+From Coq Require Import ZArith Reals Lia.
+From Coq Require SpecFloat.
 
-Notation location := location (only parsing).
-Notation loc_Exact := loc_Exact (only parsing).
-Notation loc_Inexact := loc_Inexact (only parsing).
+Require Import Zaux Raux Defs Float_prop.
+
+Notation location := SpecFloat.location (only parsing).
+Notation loc_Exact := SpecFloat.loc_Exact (only parsing).
+Notation loc_Inexact := SpecFloat.loc_Inexact (only parsing).
 
 Section Fcalc_bracket.
 
@@ -533,16 +534,35 @@ Qed.
 
 End Fcalc_bracket_step.
 
-Section Fcalc_bracket_scale.
+Section Bracket_plus.
 
-Lemma inbetween_mult_aux :
-  forall x d s,
-  ((x * s + d * s) / 2 = (x + d) / 2 * s)%R.
+Theorem inbetween_plus_compat :
+  forall x d u l t,
+  inbetween x d u l ->
+  inbetween (x + t) (d + t) (u + t) l.
 Proof.
-intros x d s.
-field.
+intros x d u l t [Hx|l' Hx Hl] ; constructor.
+now rewrite Hx.
+now split ; apply Rplus_lt_compat_r.
+replace ((x + t + (d + t)) / 2)%R with ((x + d) / 2 + t)%R by field.
+now rewrite Rcompare_plus_r.
 Qed.
 
+Theorem inbetween_plus_reg :
+  forall x d u l t,
+  inbetween (x + t) (d + t) (u + t) l ->
+  inbetween x d u l.
+Proof.
+intros x d u l t H.
+generalize (inbetween_plus_compat _ _ _ _ (Ropp t) H).
+assert (K: forall y, (y + t + -t = y)%R) by (intros y ; ring).
+now rewrite 3!K.
+Qed.
+
+End Bracket_plus.
+
+Section Fcalc_bracket_scale.
+
 Theorem inbetween_mult_compat :
   forall x d u l s,
   (0 < s)%R ->
@@ -552,7 +572,7 @@ Proof.
 intros x d u l s Hs [Hx|l' Hx Hl] ; constructor.
 now rewrite Hx.
 now split ; apply Rmult_lt_compat_r.
-rewrite inbetween_mult_aux.
+replace ((x * s + d * s) / 2)%R with ((x + d) / 2 * s)%R by field.
 now rewrite Rcompare_mult_r.
 Qed.
 
@@ -562,12 +582,11 @@ Theorem inbetween_mult_reg :
   inbetween (x * s) (d * s) (u * s) l ->
   inbetween x d u l.
 Proof.
-intros x d u l s Hs [Hx|l' Hx Hl] ; constructor.
-apply Rmult_eq_reg_r with (1 := Hx).
-now apply Rgt_not_eq.
-now split ; apply Rmult_lt_reg_r with s.
-rewrite <- Rcompare_mult_r with (1 := Hs).
-now rewrite inbetween_mult_aux in Hl.
+intros x d u l s Hs H.
+generalize (inbetween_mult_compat _ _ _ _ _ (Rinv_0_lt_compat s Hs) H).
+assert (K: forall y, (y * s * /s = y)%R).
+{ intros y. field. now apply Rgt_not_eq. }
+now rewrite 3!K.
 Qed.
 
 End Fcalc_bracket_scale.
diff --git a/flocq/Calc/Div.v b/flocq/Calc/Div.v
index 48e3bb51..88d99a1f 100644
--- a/flocq/Calc/Div.v
+++ b/flocq/Calc/Div.v
@@ -19,8 +19,9 @@ COPYING file for more details.
 
 (** * Helper function and theorem for computing the rounded quotient of two floating-point numbers. *)
 
-From Coq Require Import Lia.
-Require Import Raux Defs Generic_fmt Float_prop Digits Bracket.
+From Coq Require Import ZArith Reals Lia.
+
+Require Import Zaux Raux Defs Generic_fmt Float_prop Digits Bracket.
 
 Set Implicit Arguments.
 Set Strongly Strict Implicit.
diff --git a/flocq/Calc/Operations.v b/flocq/Calc/Operations.v
index ac93d412..bcc93f6a 100644
--- a/flocq/Calc/Operations.v
+++ b/flocq/Calc/Operations.v
@@ -19,8 +19,9 @@ COPYING file for more details.
 
 (** * Basic operations on floats: alignment, addition, multiplication *)
 
-From Coq Require Import Lia.
-Require Import Raux Defs Float_prop.
+From Coq Require Import ZArith Reals Lia.
+
+Require Import Zaux Raux Defs Float_prop.
 
 Set Implicit Arguments.
 Set Strongly Strict Implicit.
diff --git a/flocq/Calc/Plus.v b/flocq/Calc/Plus.v
new file mode 100644
index 00000000..bd338af8
--- /dev/null
+++ b/flocq/Calc/Plus.v
@@ -0,0 +1,171 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2021 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2021 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.
+*)
+
+(** * Helper function and theorem for computing the rounded sum of two floating-point numbers. *)
+
+From Coq Require Import ZArith Reals Lia.
+
+Require Import Core Bracket Operations Round.
+
+Section Plus.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+Variable fexp : Z -> Z.
+
+Context { monotone_exp : Monotone_exp fexp }.
+
+Definition Fplus_core m1 e1 m2 e2 e :=
+  let k := (e - e2)%Z in
+  let '(m2', _, l) :=
+    if Zlt_bool 0 k then truncate_aux beta (m2, e2, loc_Exact) k
+    else (m2 * Zpower beta (-k), e, loc_Exact)%Z in
+  let m1' := (m1 * Zpower beta (e1 - e))%Z in
+  (m1' + m2', l)%Z.
+
+Theorem Fplus_core_correct :
+  forall m1 e1 m2 e2 e,
+  (e <= e1)%Z ->
+  let '(m, l) := Fplus_core m1 e1 m2 e2 e in
+  inbetween_float beta m e (F2R (Float beta m1 e1) + F2R (Float beta m2 e2)) l.
+Proof.
+intros m1 e1 m2 e2 e He1.
+unfold Fplus_core.
+case Zlt_bool_spec ; intros He2.
+- unfold truncate_aux.
+  unfold inbetween_float, F2R. simpl.
+  rewrite 3!plus_IZR.
+  rewrite Rplus_assoc.
+  rewrite 2!Rmult_plus_distr_r.
+  replace (IZR (m1 * Zpower beta (e1 - e)) * bpow e)%R with (IZR m1 * bpow e1)%R.
+  2: {
+    rewrite mult_IZR, IZR_Zpower by lia.
+    rewrite Rmult_assoc, <- bpow_plus.
+    apply (f_equal (fun v => IZR m1 * bpow v)%R).
+    ring.
+  }
+  rewrite <- 3!(Rplus_comm _ (IZR m1 * bpow e1)).
+  apply inbetween_plus_compat.
+  set (k := (e - e2)%Z).
+  rewrite <- (plus_IZR _ 1).
+  replace e with (e2 + k)%Z by (unfold k ; ring).
+  apply inbetween_float_new_location.
+  exact He2.
+  now constructor 1.
+- constructor 1.
+  unfold F2R. simpl.
+  rewrite plus_IZR, Rmult_plus_distr_r.
+  rewrite 2!mult_IZR, 2!IZR_Zpower by lia.
+  rewrite 2!Rmult_assoc, <- 2!bpow_plus.
+  apply (f_equal2 (fun v w => IZR m1 * bpow v + IZR m2 * bpow w)%R) ; ring.
+Qed.
+
+Definition Fplus (f1 f2 : float beta) :=
+  let (m1, e1) := f1 in
+  let (m2, e2) := f2 in
+  if Zeq_bool m1 0 then
+    (m2, e2, loc_Exact)
+  else if Zeq_bool m2 0 then
+    (m1, e1, loc_Exact)
+  else
+  let p1 := (Zdigits beta m1 + e1)%Z in
+  let p2 := (Zdigits beta m2 + e2)%Z in
+  if Zle_bool 2 (Z.abs (p1 - p2)) then
+    let e := Z.min (Z.max e1 e2) (fexp (Z.max p1 p2 - 1)) in
+    let (m, l) :=
+      if Zlt_bool e1 e then
+        Fplus_core m2 e2 m1 e1 e
+      else
+        Fplus_core m1 e1 m2 e2 e in
+    (m, e, l)
+  else
+    let (m, e) := Fplus f1 f2 in
+    (m, e, loc_Exact).
+
+Theorem Fplus_correct :
+  forall x y,
+  let '(m, e, l) := Fplus x y in
+  (l = loc_Exact \/ e <= cexp beta fexp (F2R x + F2R y))%Z /\
+  inbetween_float beta m e (F2R x + F2R y) l.
+Proof.
+intros [m1 e1] [m2 e2].
+unfold Fplus.
+case Zeq_bool_spec ; intros Hm1.
+{ rewrite Hm1.
+  split.
+  now left.
+  rewrite F2R_0, Rplus_0_l.
+  now constructor 1. }
+case Zeq_bool_spec ; intros Hm2.
+{ rewrite Hm2.
+  split.
+  now left.
+  rewrite F2R_0, Rplus_0_r.
+  now constructor 1. }
+set (p1 := (Zdigits beta m1 + e1)%Z).
+set (p2 := (Zdigits beta m2 + e2)%Z).
+set (e := Z.min (Z.max e1 e2) (fexp (Z.max p1 p2 - 1))).
+case Zle_bool_spec ; intros Hp ; cycle 1.
+{ rewrite <- F2R_plus.
+  destruct Operations.Fplus as [m' e'].
+  split.
+  now left.
+  now constructor 1. }
+set (z := (F2R _ + F2R _)%R).
+assert (Hz: (e <= cexp beta fexp z)%Z).
+{ apply Z.le_trans with (fexp (Z.max p1 p2 - 1)).
+  apply Z.le_min_r.
+  unfold cexp.
+  apply monotone_exp.
+  unfold z.
+  revert Hp.
+  unfold p1, p2.
+  rewrite <- 2!mag_F2R_Zdigits by easy.
+  clear -Hm1 Hm2.
+  destruct (Zle_or_lt (mag beta (F2R (Float beta m1 e1))) (mag beta (F2R (Float beta m2 e2)))) as [Hp'|Hp'].
+  - rewrite Z.max_r by easy.
+    rewrite Z.abs_neq by (clear -Hp'; lia).
+    rewrite Rplus_comm.
+    intros H.
+    apply mag_plus_ge.
+    now apply F2R_neq_0.
+    clear -H ; lia.
+  - rewrite Z.max_l, Z.abs_eq by (clear -Hp'; lia).
+    intros H.
+    apply mag_plus_ge.
+    now apply F2R_neq_0.
+    clear -H ; lia. }
+case Zlt_bool_spec ; intros He.
+- assert (He': (e <= e2)%Z) by (clear -He ; lia).
+  generalize (Fplus_core_correct m2 e2 m1 e1 e He').
+  rewrite Rplus_comm.
+  fold z.
+  destruct Fplus_core as [m' l].
+  refine (fun H => conj _ H).
+  now right.
+- assert (He': (e <= e1)%Z) by (clear -He ; lia).
+  generalize (Fplus_core_correct m1 e1 m2 e2 e He').
+  fold z.
+  destruct Fplus_core as [m' l].
+  refine (fun H => conj _ H).
+  now right.
+Qed.
+
+End Plus.
diff --git a/flocq/Calc/Round.v b/flocq/Calc/Round.v
index 704a1ab2..b4693ed7 100644
--- a/flocq/Calc/Round.v
+++ b/flocq/Calc/Round.v
@@ -19,7 +19,8 @@ COPYING file for more details.
 
 (** * Helper function for computing the rounded value of a real number. *)
 
-From Coq Require Import Lia.
+From Coq Require Import ZArith Reals Lia.
+
 Require Import Core Digits Float_prop Bracket.
 
 Section Fcalc_round.
@@ -113,6 +114,12 @@ Qed.
 
 Definition cond_incr (b : bool) m := if b then (m + 1)%Z else m.
 
+Lemma le_cond_incr_le :
+  forall b m, (m <= cond_incr b m <= m + 1)%Z.
+Proof.
+unfold cond_incr; intros b; case b; lia.
+Qed.
+
 Theorem inbetween_float_round_sign :
   forall rnd choice,
   ( forall x m l, inbetween_int m (Rabs x) l ->
@@ -589,6 +596,16 @@ now case Rlt_bool.
 lia.
 Qed.
 
+Theorem inbetween_float_NA_sign :
+  forall x m l,
+  let e := cexp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp ZnearestA x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_N true l) m)) e).
+Proof.
+apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_N true l) m).
+exact inbetween_int_NA_sign.
+Qed.
+
 Definition truncate_aux t k :=
   let '(m, e, l) := t in
   let p := Zpower beta k in
diff --git a/flocq/Calc/Sqrt.v b/flocq/Calc/Sqrt.v
index 4d267d21..3c885bba 100644
--- a/flocq/Calc/Sqrt.v
+++ b/flocq/Calc/Sqrt.v
@@ -19,8 +19,9 @@ COPYING file for more details.
 
 (** * Helper functions and theorems for computing the rounded square root of a floating-point number. *)
 
-From Coq Require Import Lia.
-Require Import Raux Defs Digits Generic_fmt Float_prop Bracket.
+From Coq Require Import ZArith Reals Lia.
+
+Require Import Zaux Raux Defs Digits Generic_fmt Float_prop Bracket.
 
 Set Implicit Arguments.
 Set Strongly Strict Implicit.