From 5312915c1b29929f82e1f8de80609a277584913f Mon Sep 17 00:00:00 2001
From: xleroy <xleroy@fca1b0fc-160b-0410-b1d3-a4f43f01ea2e>
Date: Thu, 28 Jun 2012 07:59:03 +0000
Subject: Use Flocq for floats

git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@1939 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e
---
 flocq/Calc/Fcalc_bracket.v |  688 ++++++++++++++++++++++++++++
 flocq/Calc/Fcalc_digits.v  |  382 ++++++++++++++++
 flocq/Calc/Fcalc_div.v     |  165 +++++++
 flocq/Calc/Fcalc_ops.v     |  161 +++++++
 flocq/Calc/Fcalc_round.v   | 1093 ++++++++++++++++++++++++++++++++++++++++++++
 flocq/Calc/Fcalc_sqrt.v    |  346 ++++++++++++++
 6 files changed, 2835 insertions(+)
 create mode 100644 flocq/Calc/Fcalc_bracket.v
 create mode 100644 flocq/Calc/Fcalc_digits.v
 create mode 100644 flocq/Calc/Fcalc_div.v
 create mode 100644 flocq/Calc/Fcalc_ops.v
 create mode 100644 flocq/Calc/Fcalc_round.v
 create mode 100644 flocq/Calc/Fcalc_sqrt.v

(limited to 'flocq/Calc')

diff --git a/flocq/Calc/Fcalc_bracket.v b/flocq/Calc/Fcalc_bracket.v
new file mode 100644
index 00000000..dd4bd97f
--- /dev/null
+++ b/flocq/Calc/Fcalc_bracket.v
@@ -0,0 +1,688 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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.
+*)
+
+(** * Locations: where a real number is positioned with respect to its rounded-down value in an arbitrary format. *)
+
+Require Import Fcore_Raux.
+Require Import Fcore_defs.
+Require Import Fcore_float_prop.
+
+Section Fcalc_bracket.
+
+Variable d u : R.
+Hypothesis Hdu : (d < u)%R.
+
+Inductive location := loc_Exact | loc_Inexact : comparison -> location.
+
+Variable x : R.
+
+Definition inbetween_loc :=
+  match Rcompare x d with
+  | Gt => loc_Inexact (Rcompare x ((d + u) / 2))
+  | _ => loc_Exact
+  end.
+
+(** Locates a real number with respect to the middle of two other numbers. *)
+
+Inductive inbetween : location -> Prop :=
+  | inbetween_Exact : x = d -> inbetween loc_Exact
+  | inbetween_Inexact l : (d < x < u)%R -> Rcompare x ((d + u) / 2)%R = l -> inbetween (loc_Inexact l).
+
+Theorem inbetween_spec :
+  (d <= x < u)%R -> inbetween inbetween_loc.
+Proof.
+intros Hx.
+unfold inbetween_loc.
+destruct (Rcompare_spec x d) as [H|H|H].
+now elim Rle_not_lt with (1 := proj1 Hx).
+now constructor.
+constructor.
+now split.
+easy.
+Qed.
+
+Theorem inbetween_unique :
+  forall l l',
+  inbetween l -> inbetween l' -> l = l'.
+Proof.
+intros l l' Hl Hl'.
+inversion_clear Hl ; inversion_clear Hl'.
+apply refl_equal.
+rewrite H in H0.
+elim Rlt_irrefl with (1 := proj1 H0).
+rewrite H1 in H.
+elim Rlt_irrefl with (1 := proj1 H).
+apply f_equal.
+now rewrite <- H0.
+Qed.
+
+Section Fcalc_bracket_any.
+
+Variable l : location.
+
+Theorem inbetween_bounds :
+  inbetween l ->
+  (d <= x < u)%R.
+Proof.
+intros [Hx|l' Hx Hl] ; clear l.
+rewrite Hx.
+split.
+apply Rle_refl.
+exact Hdu.
+now split ; try apply Rlt_le.
+Qed.
+
+Theorem inbetween_bounds_not_Eq :
+  inbetween l ->
+  l <> loc_Exact ->
+  (d < x < u)%R.
+Proof.
+intros [Hx|l' Hx Hl] H.
+now elim H.
+exact Hx.
+Qed.
+
+End Fcalc_bracket_any.
+
+Theorem inbetween_distance_inexact :
+  forall l,
+  inbetween (loc_Inexact l) ->
+  Rcompare (x - d) (u - x) = l.
+Proof.
+intros l Hl.
+inversion_clear Hl as [|l' Hl' Hx].
+now rewrite Rcompare_middle.
+Qed.
+
+Theorem inbetween_distance_inexact_abs :
+  forall l,
+  inbetween (loc_Inexact l) ->
+  Rcompare (Rabs (d - x)) (Rabs (u - x)) = l.
+Proof.
+intros l Hl.
+rewrite Rabs_left1.
+rewrite Rabs_pos_eq.
+rewrite Ropp_minus_distr.
+now apply inbetween_distance_inexact.
+apply Rle_0_minus.
+apply Rlt_le.
+apply (inbetween_bounds _ Hl).
+apply Rle_minus.
+apply (inbetween_bounds _ Hl).
+Qed.
+
+End Fcalc_bracket.
+
+Theorem inbetween_ex :
+  forall d u l,
+  (d < u)%R ->
+  exists x,
+  inbetween d u x l.
+Proof.
+intros d u [|l] Hdu.
+exists d.
+now constructor.
+exists (d + match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4 * (u - d))%R.
+constructor.
+(* *)
+set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end / 4)%R).
+assert (0 < v < 1)%R.
+split.
+unfold v, Rdiv.
+apply Rmult_lt_0_compat.
+case l.
+now apply (Z2R_lt 0 2).
+now apply (Z2R_lt 0 1).
+now apply (Z2R_lt 0 3).
+apply Rinv_0_lt_compat.
+now apply (Z2R_lt 0 4).
+unfold v, Rdiv.
+apply Rmult_lt_reg_r with 4%R.
+now apply (Z2R_lt 0 4).
+rewrite Rmult_assoc, Rinv_l.
+rewrite Rmult_1_r, Rmult_1_l.
+case l.
+now apply (Z2R_lt 2 4).
+now apply (Z2R_lt 1 4).
+now apply (Z2R_lt 3 4).
+apply Rgt_not_eq.
+now apply (Z2R_lt 0 4).
+split.
+apply Rplus_lt_reg_r with (d * (v - 1))%R.
+ring_simplify.
+rewrite Rmult_comm.
+now apply Rmult_lt_compat_l.
+apply Rplus_lt_reg_r with (-u * v)%R.
+replace (- u * v + (d + v * (u - d)))%R with (d * (1 - v))%R by ring.
+replace (- u * v + u)%R with (u * (1 - v))%R by ring.
+apply Rmult_lt_compat_r.
+apply Rplus_lt_reg_r with v.
+now ring_simplify.
+exact Hdu.
+(* *)
+set (v := (match l with Lt => 1 | Eq => 2 | Gt => 3 end)%R).
+rewrite <- (Rcompare_plus_r (- (d + u) / 2)).
+rewrite <- (Rcompare_mult_r 4).
+2: now apply (Z2R_lt 0 4).
+replace (((d + u) / 2 + - (d + u) / 2) * 4)%R with ((u - d) * 0)%R by field.
+replace ((d + v / 4 * (u - d) + - (d + u) / 2) * 4)%R with ((u - d) * (v - 2))%R by field.
+rewrite Rcompare_mult_l.
+2: now apply Rlt_Rminus.
+rewrite <- (Rcompare_plus_r 2).
+ring_simplify (v - 2 + 2)%R (0 + 2)%R.
+unfold v.
+case l.
+exact (Rcompare_Z2R 2 2).
+exact (Rcompare_Z2R 1 2).
+exact (Rcompare_Z2R 3 2).
+Qed.
+
+Section Fcalc_bracket_step.
+
+Variable start step : R.
+Variable nb_steps : Z.
+Variable Hstep : (0 < step)%R.
+
+Lemma ordered_steps :
+  forall k,
+  (start + Z2R k * step < start + Z2R (k + 1) * step)%R.
+Proof.
+intros k.
+apply Rplus_lt_compat_l.
+apply Rmult_lt_compat_r.
+exact Hstep.
+apply Z2R_lt.
+apply Zlt_succ.
+Qed.
+
+Lemma middle_range :
+  forall k,
+  ((start + (start + Z2R k * step)) / 2 = start + (Z2R k / 2 * step))%R.
+Proof.
+intros k.
+field.
+Qed.
+
+Hypothesis (Hnb_steps : (1 < nb_steps)%Z).
+
+Lemma inbetween_step_not_Eq :
+  forall x k l l',
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  (0 < k < nb_steps)%Z ->
+  Rcompare x (start + (Z2R nb_steps / 2 * step))%R = l' ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact l').
+Proof.
+intros x k l l' Hx Hk Hl'.
+constructor.
+(* . *)
+assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx).
+split.
+apply Rlt_le_trans with (2 := proj1 Hx').
+rewrite <- (Rplus_0_r start) at 1.
+apply Rplus_lt_compat_l.
+apply Rmult_lt_0_compat.
+now apply (Z2R_lt 0).
+exact Hstep.
+apply Rlt_le_trans with (1 := proj2 Hx').
+apply Rplus_le_compat_l.
+apply Rmult_le_compat_r.
+now apply Rlt_le.
+apply Z2R_le.
+omega.
+(* . *)
+now rewrite middle_range.
+Qed.
+
+Theorem inbetween_step_Lo :
+  forall x k l,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  (0 < k)%Z -> (2 * k + 1 < nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt).
+Proof.
+intros x k l Hx Hk1 Hk2.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+apply Rcompare_Lt.
+assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx).
+apply Rlt_le_trans with (1 := proj2 Hx').
+apply Rcompare_not_Lt_inv.
+rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l.
+apply Rcompare_not_Lt.
+change 2%R with (Z2R 2).
+rewrite <- Z2R_mult.
+apply Z2R_le.
+omega.
+exact Hstep.
+Qed.
+
+Theorem inbetween_step_Hi :
+  forall x k l,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  (nb_steps < 2 * k)%Z -> (k < nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Gt).
+Proof.
+intros x k l Hx Hk1 Hk2.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+apply Rcompare_Gt.
+assert (Hx' := inbetween_bounds _ _ (ordered_steps _) _ _ Hx).
+apply Rlt_le_trans with (2 := proj1 Hx').
+apply Rcompare_Lt_inv.
+rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l.
+apply Rcompare_Lt.
+change 2%R with (Z2R 2).
+rewrite <- Z2R_mult.
+apply Z2R_lt.
+omega.
+exact Hstep.
+Qed.
+
+Theorem inbetween_step_Lo_not_Eq :
+  forall x l,
+  inbetween start (start + step) x l ->
+  l <> loc_Exact ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt).
+Proof.
+intros x l Hx Hl.
+assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl).
+constructor.
+(* . *)
+split.
+apply Hx'.
+apply Rlt_trans with (1 := proj2 Hx').
+apply Rplus_lt_compat_l.
+rewrite <- (Rmult_1_l step) at 1.
+apply Rmult_lt_compat_r.
+exact Hstep.
+now apply (Z2R_lt 1).
+(* . *)
+apply Rcompare_Lt.
+apply Rlt_le_trans with (1 := proj2 Hx').
+rewrite middle_range.
+apply Rcompare_not_Lt_inv.
+rewrite <- (Rmult_1_l step) at 2.
+rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_l.
+rewrite Rmult_1_r.
+apply Rcompare_not_Lt.
+apply (Z2R_le 2).
+now apply (Zlt_le_succ 1).
+exact Hstep.
+Qed.
+
+Lemma middle_odd :
+  forall k,
+  (2 * k + 1 = nb_steps)%Z ->
+  (((start + Z2R k * step) + (start + Z2R (k + 1) * step))/2 = start + Z2R nb_steps /2 * step)%R.
+Proof.
+intros k Hk.
+rewrite <- Hk.
+rewrite 2!Z2R_plus, Z2R_mult.
+simpl. field.
+Qed.
+
+Theorem inbetween_step_any_Mi_odd :
+  forall x k l,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x (loc_Inexact l) ->
+  (2 * k + 1 = nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact l).
+Proof.
+intros x k l Hx Hk.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+inversion_clear Hx as [|l' _ Hl].
+now rewrite (middle_odd _ Hk) in Hl.
+Qed.
+
+Theorem inbetween_step_Lo_Mi_Eq_odd :
+  forall x k,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x loc_Exact ->
+  (2 * k + 1 = nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Lt).
+Proof.
+intros x k Hx Hk.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+inversion_clear Hx as [Hl|].
+rewrite Hl.
+rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r.
+apply Rcompare_Lt.
+change 2%R with (Z2R 2).
+rewrite <- Z2R_mult.
+apply Z2R_lt.
+rewrite <- Hk.
+apply Zlt_succ.
+exact Hstep.
+Qed.
+
+Theorem inbetween_step_Hi_Mi_even :
+  forall x k l,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  l <> loc_Exact ->
+  (2 * k = nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Gt).
+Proof.
+intros x k l Hx Hl Hk.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+apply Rcompare_Gt.
+assert (Hx' := inbetween_bounds_not_Eq _ _ _ _ Hx Hl).
+apply Rle_lt_trans with (2 := proj1 Hx').
+apply Rcompare_not_Lt_inv.
+rewrite Rcompare_plus_l, Rcompare_mult_r, Rcompare_half_r.
+change 2%R with (Z2R 2).
+rewrite <- Z2R_mult.
+apply Rcompare_not_Lt.
+apply Z2R_le.
+rewrite Hk.
+apply Zle_refl.
+exact Hstep.
+Qed.
+
+Theorem inbetween_step_Mi_Mi_even :
+  forall x k,
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x loc_Exact ->
+  (2 * k = nb_steps)%Z ->
+  inbetween start (start + Z2R nb_steps * step) x (loc_Inexact Eq).
+Proof.
+intros x k Hx Hk.
+apply inbetween_step_not_Eq with (1 := Hx).
+omega.
+apply Rcompare_Eq.
+inversion_clear Hx as [Hx'|].
+rewrite Hx', <- Hk, Z2R_mult.
+simpl (Z2R 2).
+field.
+Qed.
+
+(** Computes a new location when the interval containing a real
+    number is split into nb_steps subintervals and the real is
+    in the k-th one. (Even radix.) *)
+
+Definition new_location_even k l :=
+  if Zeq_bool k 0 then
+    match l with loc_Exact => l | _ => loc_Inexact Lt end
+  else
+    loc_Inexact
+    match Zcompare (2 * k) nb_steps with
+    | Lt => Lt
+    | Eq => match l with loc_Exact => Eq | _ => Gt end
+    | Gt => Gt
+    end.
+
+Theorem new_location_even_correct :
+  Zeven nb_steps = true ->
+  forall x k l, (0 <= k < nb_steps)%Z ->
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  inbetween start (start + Z2R nb_steps * step) x (new_location_even k l).
+Proof.
+intros He x k l Hk Hx.
+unfold new_location_even.
+destruct (Zeq_bool_spec k 0) as [Hk0|Hk0].
+(* k = 0 *)
+rewrite Hk0 in Hx.
+rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx.
+set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end).
+assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)).
+unfold l' ; case l ; try (now left) ; right ; now split.
+destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2.
+constructor.
+rewrite H1 in Hx.
+now inversion_clear Hx.
+now apply inbetween_step_Lo_not_Eq with (2 := H1).
+(* k <> 0 *)
+destruct (Zcompare_spec (2 * k) nb_steps) as [Hk1|Hk1|Hk1].
+(* . 2 * k < nb_steps *)
+apply inbetween_step_Lo with (1 := Hx).
+omega.
+destruct (Zeven_ex nb_steps).
+rewrite He in H.
+omega.
+(* . 2 * k = nb_steps *)
+set (l' := match l with loc_Exact => Eq | _ => Gt end).
+assert ((l = loc_Exact /\ l' = Eq) \/ (l <> loc_Exact /\ l' = Gt)).
+unfold l' ; case l ; try (now left) ; right ; now split.
+destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2.
+rewrite H1 in Hx.
+now apply inbetween_step_Mi_Mi_even with (1 := Hx).
+now apply inbetween_step_Hi_Mi_even with (1 := Hx).
+(* . 2 * k > nb_steps *)
+apply inbetween_step_Hi with (1 := Hx).
+exact Hk1.
+apply Hk.
+Qed.
+
+(** Computes a new location when the interval containing a real
+    number is split into nb_steps subintervals and the real is
+    in the k-th one. (Odd radix.) *)
+
+Definition new_location_odd k l :=
+  if Zeq_bool k 0 then
+    match l with loc_Exact => l | _ => loc_Inexact Lt end
+  else
+    loc_Inexact
+    match Zcompare (2 * k + 1) nb_steps with
+    | Lt => Lt
+    | Eq => match l with loc_Inexact l => l | loc_Exact => Lt end
+    | Gt => Gt
+    end.
+
+Theorem new_location_odd_correct :
+  Zeven nb_steps = false ->
+  forall x k l, (0 <= k < nb_steps)%Z ->
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  inbetween start (start + Z2R nb_steps * step) x (new_location_odd k l).
+Proof.
+intros Ho x k l Hk Hx.
+unfold new_location_odd.
+destruct (Zeq_bool_spec k 0) as [Hk0|Hk0].
+(* k = 0 *)
+rewrite Hk0 in Hx.
+rewrite Rmult_0_l, Rplus_0_r, Rmult_1_l in Hx.
+set (l' := match l with loc_Exact => l | _ => loc_Inexact Lt end).
+assert ((l = loc_Exact /\ l' = loc_Exact) \/ (l <> loc_Exact /\ l' = loc_Inexact Lt)).
+unfold l' ; case l ; try (now left) ; right ; now split.
+destruct H as [(H1,H2)|(H1,H2)] ; rewrite H2.
+constructor.
+rewrite H1 in Hx.
+now inversion_clear Hx.
+now apply inbetween_step_Lo_not_Eq with (2 := H1).
+(* k <> 0 *)
+destruct (Zcompare_spec (2 * k + 1) nb_steps) as [Hk1|Hk1|Hk1].
+(* . 2 * k + 1 < nb_steps *)
+apply inbetween_step_Lo with (1 := Hx) (3 := Hk1).
+omega.
+(* . 2 * k + 1 = nb_steps *)
+destruct l.
+apply inbetween_step_Lo_Mi_Eq_odd with (1 := Hx) (2 := Hk1).
+apply inbetween_step_any_Mi_odd with (1 := Hx) (2 := Hk1).
+(* . 2 * k + 1 > nb_steps *)
+apply inbetween_step_Hi with (1 := Hx).
+destruct (Zeven_ex nb_steps).
+rewrite Ho in H.
+omega.
+apply Hk.
+Qed.
+
+Definition new_location :=
+  if Zeven nb_steps then new_location_even else new_location_odd.
+
+Theorem new_location_correct :
+  forall x k l, (0 <= k < nb_steps)%Z ->
+  inbetween (start + Z2R k * step) (start + Z2R (k + 1) * step) x l ->
+  inbetween start (start + Z2R nb_steps * step) x (new_location k l).
+Proof.
+intros x k l Hk Hx.
+unfold new_location.
+generalize (refl_equal nb_steps) (Zle_lt_trans _ _ _ (proj1 Hk) (proj2 Hk)).
+pattern nb_steps at 2 3 5 ; case nb_steps.
+now intros _.
+(* . *)
+intros [p|p|] Hp _.
+apply new_location_odd_correct with (2 := Hk) (3 := Hx).
+now rewrite Hp.
+apply new_location_even_correct with (2 := Hk) (3 := Hx).
+now rewrite Hp.
+now rewrite Hp in Hnb_steps.
+(* . *)
+now intros p _.
+Qed.
+
+End Fcalc_bracket_step.
+
+Section Fcalc_bracket_scale.
+
+Lemma inbetween_mult_aux :
+  forall x d s,
+  ((x * s + d * s) / 2 = (x + d) / 2 * s)%R.
+Proof.
+intros x d s.
+field.
+Qed.
+
+Theorem inbetween_mult_compat :
+  forall x d u l s,
+  (0 < s)%R ->
+  inbetween x d u l ->
+  inbetween (x * s) (d * s) (u * s) l.
+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.
+now rewrite Rcompare_mult_r.
+Qed.
+
+Theorem inbetween_mult_reg :
+  forall x d u l s,
+  (0 < s)%R ->
+  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.
+Qed.
+
+End Fcalc_bracket_scale.
+
+Section Fcalc_bracket_generic.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+(** Specialization of inbetween for two consecutive floating-point numbers. *)
+
+Definition inbetween_float m e x l :=
+  inbetween (F2R (Float beta m e)) (F2R (Float beta (m + 1) e)) x l.
+
+Theorem inbetween_float_bounds :
+  forall x m e l,
+  inbetween_float m e x l ->
+  (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R.
+Proof.
+intros x m e l [Hx|l' Hx Hl].
+rewrite Hx.
+split.
+apply Rle_refl.
+apply F2R_lt_compat.
+apply Zlt_succ.
+split.
+now apply Rlt_le.
+apply Hx.
+Qed.
+
+(** Specialization of inbetween for two consecutive integers. *)
+
+Definition inbetween_int m x l :=
+  inbetween (Z2R m) (Z2R (m + 1)) x l.
+
+Theorem inbetween_float_new_location :
+  forall x m e l k,
+  (0 < k)%Z ->
+  inbetween_float m e x l ->
+  inbetween_float (Zdiv m (Zpower beta k)) (e + k) x (new_location (Zpower beta k) (Zmod m (Zpower beta k)) l).
+Proof.
+intros x m e l k Hk Hx.
+unfold inbetween_float in *.
+assert (Hr: forall m, F2R (Float beta m (e + k)) = F2R (Float beta (m * Zpower beta k) e)).
+clear -Hk. intros m.
+rewrite (F2R_change_exp beta e).
+apply (f_equal (fun r => F2R (Float beta (m * Zpower _ r) e))).
+ring.
+omega.
+assert (Hp: (Zpower beta k > 0)%Z).
+apply Zlt_gt.
+apply Zpower_gt_0.
+now apply Zlt_le_weak.
+(* . *)
+rewrite 2!Hr.
+rewrite Zmult_plus_distr_l, Zmult_1_l.
+unfold F2R at 2. simpl.
+rewrite Z2R_plus, Rmult_plus_distr_r.
+apply new_location_correct.
+apply bpow_gt_0.
+now apply Zpower_gt_1.
+now apply Z_mod_lt.
+rewrite <- 2!Rmult_plus_distr_r, <- 2!Z2R_plus.
+rewrite Zmult_comm, Zplus_assoc.
+now rewrite <- Z_div_mod_eq.
+Qed.
+
+Theorem inbetween_float_new_location_single :
+  forall x m e l,
+  inbetween_float m e x l ->
+  inbetween_float (Zdiv m beta) (e + 1) x (new_location beta (Zmod m beta) l).
+Proof.
+intros x m e l Hx.
+replace (radix_val beta) with (Zpower beta 1).
+now apply inbetween_float_new_location.
+apply Zmult_1_r.
+Qed.
+
+Theorem inbetween_float_ex :
+  forall m e l,
+  exists x,
+  inbetween_float m e x l.
+Proof.
+intros m e l.
+apply inbetween_ex.
+apply F2R_lt_compat.
+apply Zlt_succ.
+Qed.
+
+Theorem inbetween_float_unique :
+  forall x e m l m' l',
+  inbetween_float m e x l ->
+  inbetween_float m' e x l' ->
+  m = m' /\ l = l'.
+Proof.
+intros x e m l m' l' H H'.
+refine ((fun Hm => conj Hm _) _).
+rewrite <- Hm in H'. clear -H H'.
+apply inbetween_unique with (1 := H) (2 := H').
+destruct (inbetween_float_bounds x m e l H) as (H1,H2).
+destruct (inbetween_float_bounds x m' e l' H') as (H3,H4).
+cut (m < m' + 1 /\ m' < m + 1)%Z. clear ; omega.
+now split ; apply F2R_lt_reg with beta e ; apply Rle_lt_trans with x.
+Qed.
+
+End Fcalc_bracket_generic.
diff --git a/flocq/Calc/Fcalc_digits.v b/flocq/Calc/Fcalc_digits.v
new file mode 100644
index 00000000..6210bac2
--- /dev/null
+++ b/flocq/Calc/Fcalc_digits.v
@@ -0,0 +1,382 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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.
+*)
+
+(** * Functions for computing the number of digits of integers and related theorems. *)
+
+Require Import Fcore_Raux.
+Require Import Fcore_defs.
+Require Import Fcore_float_prop.
+Require Import Fcore_digits.
+
+Section Fcalc_digits.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+
+
+Theorem Zdigits_ln_beta :
+  forall n,
+  n <> Z0 ->
+  Zdigits beta n = ln_beta beta (Z2R n).
+Proof.
+intros n Hn.
+destruct (ln_beta beta (Z2R n)) as (e, He) ; simpl.
+specialize (He (Z2R_neq _ _ Hn)).
+rewrite <- Z2R_abs in He.
+assert (Hd := Zdigits_correct beta n).
+assert (Hd' := Zdigits_gt_0 beta n).
+apply Zle_antisym ; apply (bpow_lt_bpow beta).
+apply Rle_lt_trans with (2 := proj2 He).
+rewrite <- Z2R_Zpower by omega.
+now apply Z2R_le.
+apply Rle_lt_trans with (1 := proj1 He).
+rewrite <- Z2R_Zpower by omega.
+now apply Z2R_lt.
+Qed.
+
+Theorem ln_beta_F2R_Zdigits :
+  forall m e, m <> Z0 ->
+  (ln_beta beta (F2R (Float beta m e)) = Zdigits beta m + e :> Z)%Z.
+Proof.
+intros m e Hm.
+rewrite ln_beta_F2R with (1 := Hm).
+apply (f_equal (fun v => Zplus v e)).
+apply sym_eq.
+now apply Zdigits_ln_beta.
+Qed.
+
+Theorem Zdigits_mult_Zpower :
+  forall m e,
+  m <> Z0 -> (0 <= e)%Z ->
+  Zdigits beta (m * Zpower beta e) = (Zdigits beta m + e)%Z.
+Proof.
+intros m e Hm He.
+rewrite <- ln_beta_F2R_Zdigits with (1 := Hm).
+rewrite Zdigits_ln_beta.
+rewrite Z2R_mult.
+now rewrite Z2R_Zpower with (1 := He).
+contradict Hm.
+apply Zmult_integral_l with (2 := Hm).
+apply neq_Z2R.
+rewrite Z2R_Zpower with (1 := He).
+apply Rgt_not_eq.
+apply bpow_gt_0.
+Qed.
+
+Theorem Zdigits_Zpower :
+  forall e,
+  (0 <= e)%Z ->
+  Zdigits beta (Zpower beta e) = (e + 1)%Z.
+Proof.
+intros e He.
+rewrite <- (Zmult_1_l (Zpower beta e)).
+rewrite Zdigits_mult_Zpower ; try easy.
+apply Zplus_comm.
+Qed.
+
+Theorem Zdigits_le :
+  forall x y,
+  (0 <= x)%Z -> (x <= y)%Z ->
+  (Zdigits beta x <= Zdigits beta y)%Z.
+Proof.
+intros x y Hx Hxy.
+case (Z_lt_le_dec 0 x).
+clear Hx. intros Hx.
+rewrite 2!Zdigits_ln_beta.
+apply ln_beta_le.
+now apply (Z2R_lt 0).
+now apply Z2R_le.
+apply Zgt_not_eq.
+now apply Zlt_le_trans with x.
+now apply Zgt_not_eq.
+intros Hx'.
+rewrite (Zle_antisym _ _ Hx' Hx).
+apply Zdigits_ge_0.
+Qed.
+
+Theorem lt_Zdigits :
+  forall x y,
+  (0 <= y)%Z ->
+  (Zdigits beta x < Zdigits beta y)%Z ->
+  (x < y)%Z.
+Proof.
+intros x y Hy.
+cut (y <= x -> Zdigits beta y <= Zdigits beta x)%Z. omega.
+now apply Zdigits_le.
+Qed.
+
+Theorem Zpower_le_Zdigits :
+  forall e x,
+  (e < Zdigits beta x)%Z ->
+  (Zpower beta e <= Zabs x)%Z.
+Proof.
+intros e x Hex.
+destruct (Zdigits_correct beta x) as (H1,H2).
+apply Zle_trans with (2 := H1).
+apply Zpower_le.
+clear -Hex ; omega.
+Qed.
+
+Theorem Zdigits_le_Zpower :
+  forall e x,
+  (Zabs x < Zpower beta e)%Z ->
+  (Zdigits beta x <= e)%Z.
+Proof.
+intros e x.
+generalize (Zpower_le_Zdigits e x).
+omega.
+Qed.
+
+Theorem Zpower_gt_Zdigits :
+  forall e x,
+  (Zdigits beta x <= e)%Z ->
+  (Zabs x < Zpower beta e)%Z.
+Proof.
+intros e x Hex.
+destruct (Zdigits_correct beta x) as (H1,H2).
+apply Zlt_le_trans with (1 := H2).
+now apply Zpower_le.
+Qed.
+
+Theorem Zdigits_gt_Zpower :
+  forall e x,
+  (Zpower beta e <= Zabs x)%Z ->
+  (e < Zdigits beta x)%Z.
+Proof.
+intros e x Hex.
+generalize (Zpower_gt_Zdigits e x).
+omega.
+Qed.
+
+(** Characterizes the number digits of a product.
+
+This strong version is needed for proofs of division and square root
+algorithms, since they involve operation remainders.
+*)
+
+Theorem Zdigits_mult_strong :
+  forall x y,
+  (0 <= x)%Z -> (0 <= y)%Z ->
+  (Zdigits beta (x + y + x * y) <= Zdigits beta x + Zdigits beta y)%Z.
+Proof.
+intros x y Hx Hy.
+case (Z_lt_le_dec 0 x).
+clear Hx. intros Hx.
+case (Z_lt_le_dec 0 y).
+clear Hy. intros Hy.
+(* . *)
+assert (Hxy: (0 < Z2R (x + y + x * y))%R).
+apply (Z2R_lt 0).
+change Z0 with (0 + 0 + 0)%Z.
+apply Zplus_lt_compat.
+now apply Zplus_lt_compat.
+now apply Zmult_lt_0_compat.
+rewrite 3!Zdigits_ln_beta ; try now (apply sym_not_eq ; apply Zlt_not_eq).
+apply ln_beta_le_bpow with (1 := Rgt_not_eq _ _ Hxy).
+rewrite Rabs_pos_eq with (1 := Rlt_le _ _ Hxy).
+destruct (ln_beta beta (Z2R x)) as (ex, Hex). simpl.
+specialize (Hex (Rgt_not_eq _ _ (Z2R_lt _ _ Hx))).
+destruct (ln_beta beta (Z2R y)) as (ey, Hey). simpl.
+specialize (Hey (Rgt_not_eq _ _ (Z2R_lt _ _ Hy))).
+apply Rlt_le_trans with (Z2R (x + 1) * Z2R (y + 1))%R.
+rewrite <- Z2R_mult.
+apply Z2R_lt.
+apply Zplus_lt_reg_r with (- (x + y + x * y + 1))%Z.
+now ring_simplify.
+rewrite bpow_plus.
+apply Rmult_le_compat ; try (apply (Z2R_le 0) ; omega).
+rewrite <- (Rmult_1_r (Z2R (x + 1))).
+change (F2R (Float beta (x + 1) 0) <= bpow ex)%R.
+apply F2R_p1_le_bpow.
+exact Hx.
+unfold F2R. rewrite Rmult_1_r.
+apply Rle_lt_trans with (Rabs (Z2R x)).
+apply RRle_abs.
+apply Hex.
+rewrite <- (Rmult_1_r (Z2R (y + 1))).
+change (F2R (Float beta (y + 1) 0) <= bpow ey)%R.
+apply F2R_p1_le_bpow.
+exact Hy.
+unfold F2R. rewrite Rmult_1_r.
+apply Rle_lt_trans with (Rabs (Z2R y)).
+apply RRle_abs.
+apply Hey.
+apply neq_Z2R.
+now apply Rgt_not_eq.
+(* . *)
+intros Hy'.
+rewrite (Zle_antisym _ _ Hy' Hy).
+rewrite Zmult_0_r, 3!Zplus_0_r.
+apply Zle_refl.
+intros Hx'.
+rewrite (Zle_antisym _ _ Hx' Hx).
+rewrite Zmult_0_l, Zplus_0_r, 2!Zplus_0_l.
+apply Zle_refl.
+Qed.
+
+Theorem Zdigits_mult :
+  forall x y,
+  (Zdigits beta (x * y) <= Zdigits beta x + Zdigits beta y)%Z.
+Proof.
+intros x y.
+rewrite <- Zdigits_abs.
+rewrite <- (Zdigits_abs _ x).
+rewrite <- (Zdigits_abs _ y).
+apply Zle_trans with (Zdigits beta (Zabs x + Zabs y + Zabs x * Zabs y)).
+apply Zdigits_le.
+apply Zabs_pos.
+rewrite Zabs_Zmult.
+generalize (Zabs_pos x) (Zabs_pos y).
+omega.
+apply Zdigits_mult_strong ; apply Zabs_pos.
+Qed.
+
+Theorem Zdigits_mult_ge :
+  forall x y,
+  (x <> 0)%Z -> (y <> 0)%Z ->
+  (Zdigits beta x + Zdigits beta y - 1 <= Zdigits beta (x * y))%Z.
+Proof.
+intros x y Zx Zy.
+cut ((Zdigits beta x - 1) + (Zdigits beta y - 1) < Zdigits beta (x * y))%Z. omega.
+apply Zdigits_gt_Zpower.
+rewrite Zabs_Zmult.
+rewrite Zpower_exp.
+apply Zmult_le_compat.
+apply Zpower_le_Zdigits.
+apply Zlt_pred.
+apply Zpower_le_Zdigits.
+apply Zlt_pred.
+apply Zpower_ge_0.
+apply Zpower_ge_0.
+generalize (Zdigits_gt_0 beta x). omega.
+generalize (Zdigits_gt_0 beta y). omega.
+Qed.
+
+Theorem Zdigits_div_Zpower :
+  forall m e,
+  (0 <= m)%Z ->
+  (0 <= e <= Zdigits beta m)%Z ->
+  Zdigits beta (m / Zpower beta e) = (Zdigits beta m - e)%Z.
+Proof.
+intros m e Hm He.
+destruct (Zle_lt_or_eq 0 m Hm) as [Hm'|Hm'].
+(* *)
+destruct (Zle_lt_or_eq _ _ (proj2 He)) as [He'|He'].
+(* . *)
+unfold Zminus.
+rewrite <- ln_beta_F2R_Zdigits.
+2: now apply Zgt_not_eq.
+assert (Hp: (0 < Zpower beta e)%Z).
+apply lt_Z2R.
+rewrite Z2R_Zpower with (1 := proj1 He).
+apply bpow_gt_0.
+rewrite Zdigits_ln_beta.
+rewrite <- Zfloor_div with (1 := Zgt_not_eq _ _ Hp).
+rewrite Z2R_Zpower with (1 := proj1 He).
+unfold Rdiv.
+rewrite <- bpow_opp.
+change (Z2R m * bpow (-e))%R with (F2R (Float beta m (-e))).
+destruct (ln_beta beta (F2R (Float beta m (-e)))) as (e', E').
+simpl.
+specialize (E' (Rgt_not_eq _ _ (F2R_gt_0_compat beta (Float beta m (-e)) Hm'))).
+apply ln_beta_unique.
+rewrite Rabs_pos_eq in E'.
+2: now apply F2R_ge_0_compat.
+rewrite Rabs_pos_eq.
+split.
+assert (H': (0 <= e' - 1)%Z).
+(* .. *)
+cut (1 <= e')%Z. omega.
+apply bpow_lt_bpow with beta.
+apply Rle_lt_trans with (2 := proj2 E').
+simpl.
+rewrite <- (Rinv_r (bpow e)).
+rewrite <- bpow_opp.
+unfold F2R. simpl.
+apply Rmult_le_compat_r.
+apply bpow_ge_0.
+rewrite <- Z2R_Zpower with (1 := proj1 He).
+apply Z2R_le.
+rewrite <- Zabs_eq with (1 := Hm).
+now apply Zpower_le_Zdigits.
+apply Rgt_not_eq.
+apply bpow_gt_0.
+(* .. *)
+rewrite <- Z2R_Zpower with (1 := H').
+apply Z2R_le.
+apply Zfloor_lub.
+rewrite Z2R_Zpower with (1 := H').
+apply E'.
+apply Rle_lt_trans with (2 := proj2 E').
+apply Zfloor_lb.
+apply (Z2R_le 0).
+apply Zfloor_lub.
+now apply F2R_ge_0_compat.
+apply Zgt_not_eq.
+apply Zlt_le_trans with (beta^e / beta^e)%Z.
+rewrite Z_div_same_full.
+apply refl_equal.
+now apply Zgt_not_eq.
+apply Z_div_le.
+now apply Zlt_gt.
+rewrite <- Zabs_eq with (1 := Hm).
+now apply Zpower_le_Zdigits.
+(* . *)
+unfold Zminus.
+rewrite He', Zplus_opp_r.
+rewrite Zdiv_small.
+apply refl_equal.
+apply conj with (1 := Hm).
+pattern m at 1 ; rewrite <- Zabs_eq with (1 := Hm).
+apply Zpower_gt_Zdigits.
+apply Zle_refl.
+(* *)
+revert He.
+rewrite <- Hm'.
+intros (H1, H2).
+simpl.
+now rewrite (Zle_antisym _ _ H2 H1).
+Qed.
+
+End Fcalc_digits.
+
+Definition radix2 := Build_radix 2 (refl_equal _).
+
+Theorem Z_of_nat_S_digits2_Pnat :
+  forall m : positive,
+  Z_of_nat (S (digits2_Pnat m)) = Zdigits radix2 (Zpos m).
+Proof.
+intros m.
+rewrite Zdigits_ln_beta. 2: easy.
+apply sym_eq.
+apply ln_beta_unique.
+generalize (digits2_Pnat m) (digits2_Pnat_correct m).
+intros d H.
+simpl in H.
+replace (Z_of_nat (S d) - 1)%Z with (Z_of_nat d).
+rewrite <- Z2R_abs.
+rewrite <- 2!Z2R_Zpower_nat.
+split.
+now apply Z2R_le.
+now apply Z2R_lt.
+rewrite inj_S.
+apply Zpred_succ.
+Qed.
+
diff --git a/flocq/Calc/Fcalc_div.v b/flocq/Calc/Fcalc_div.v
new file mode 100644
index 00000000..6594e55b
--- /dev/null
+++ b/flocq/Calc/Fcalc_div.v
@@ -0,0 +1,165 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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 quotient of two floating-point numbers. *)
+
+Require Import Fcore_Raux.
+Require Import Fcore_defs.
+Require Import Fcore_float_prop.
+Require Import Fcore_digits.
+Require Import Fcalc_bracket.
+Require Import Fcalc_digits.
+
+Section Fcalc_div.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+(** Computes a mantissa of precision p, the corresponding exponent,
+    and the position with respect to the real quotient of the
+    input floating-point numbers.
+
+The algorithm performs the following steps:
+- Shift dividend mantissa so that it has at least p2 + p digits.
+- Perform the Euclidean division.
+- Compute the position according to the division remainder.
+
+Complexity is fine as long as p1 <= 2p and p2 <= p.
+*)
+
+Definition Fdiv_core prec m1 e1 m2 e2 :=
+  let d1 := Zdigits beta m1 in
+  let d2 := Zdigits beta m2 in
+  let e := (e1 - e2)%Z in
+  let (m, e') :=
+    match (d2 + prec - d1)%Z with
+    | Zpos p => (m1 * Zpower_pos beta p, e + Zneg p)%Z
+    | _ => (m1, e)
+    end in
+  let '(q, r) :=  Zdiv_eucl m m2 in
+  (q, e', new_location m2 r loc_Exact).
+
+Theorem Fdiv_core_correct :
+  forall prec m1 e1 m2 e2,
+  (0 < prec)%Z ->
+  (0 < m1)%Z -> (0 < m2)%Z ->
+  let '(m, e, l) := Fdiv_core prec m1 e1 m2 e2 in
+  (prec <= Zdigits beta m)%Z /\
+  inbetween_float beta m e (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) l.
+Proof.
+intros prec m1 e1 m2 e2 Hprec Hm1 Hm2.
+unfold Fdiv_core.
+set (d1 := Zdigits beta m1).
+set (d2 := Zdigits beta m2).
+case_eq
+ (match (d2 + prec - d1)%Z with
+  | Zpos p => ((m1 * Zpower_pos beta p)%Z, (e1 - e2 + Zneg p)%Z)
+  | _ => (m1, (e1 - e2)%Z)
+  end).
+intros m' e' Hme.
+(* . the shifted mantissa m' has enough digits *)
+assert (Hs: F2R (Float beta m' (e' + e2)) = F2R (Float beta m1 e1) /\ (0 < m')%Z /\ (d2 + prec <= Zdigits beta m')%Z).
+replace (d2 + prec)%Z with (d2 + prec - d1 + d1)%Z by ring.
+destruct (d2 + prec - d1)%Z as [|p|p] ;
+  unfold d1 ;
+  inversion Hme.
+ring_simplify (e1 - e2 + e2)%Z.
+repeat split.
+now rewrite <- H0.
+apply Zle_refl.
+replace (e1 - e2 + Zneg p + e2)%Z with (e1 - Zpos p)%Z by (unfold Zminus ; simpl ; ring).
+fold (Zpower beta (Zpos p)).
+split.
+pattern (Zpos p) at 1 ; replace (Zpos p) with (e1 - (e1 - Zpos p))%Z by ring.
+apply sym_eq.
+apply F2R_change_exp.
+assert (0 < Zpos p)%Z by easy.
+omega.
+split.
+apply Zmult_lt_0_compat.
+exact Hm1.
+now apply Zpower_gt_0.
+rewrite Zdigits_mult_Zpower.
+rewrite Zplus_comm.
+apply Zle_refl.
+apply sym_not_eq.
+now apply Zlt_not_eq.
+easy.
+split.
+now ring_simplify (e1 - e2 + e2)%Z.
+assert (Zneg p < 0)%Z by easy.
+omega.
+(* . *)
+destruct Hs as (Hs1, (Hs2, Hs3)).
+rewrite <- Hs1.
+generalize (Z_div_mod m' m2 (Zlt_gt _ _ Hm2)).
+destruct (Zdiv_eucl m' m2) as (q, r).
+intros (Hq, Hr).
+split.
+(* . the result mantissa q has enough digits *)
+cut (Zdigits beta m' <= d2 + Zdigits beta q)%Z. omega.
+unfold d2.
+rewrite Hq.
+assert (Hq': (0 < q)%Z).
+apply Zmult_lt_reg_r with (1 := Hm2).
+assert (m2 < m')%Z.
+apply lt_Zdigits with beta.
+now apply Zlt_le_weak.
+unfold d2 in Hs3.
+clear -Hprec Hs3 ; omega.
+cut (q * m2 = m' - r)%Z. clear -Hr H ; omega.
+rewrite Hq.
+ring.
+apply Zle_trans with (Zdigits beta (m2 + q + m2 * q)).
+apply Zdigits_le.
+rewrite <- Hq.
+now apply Zlt_le_weak.
+clear -Hr Hq'. omega.
+apply Zdigits_mult_strong ; apply Zlt_le_weak.
+now apply Zle_lt_trans with r.
+exact Hq'.
+(* . the location is correctly computed *)
+unfold inbetween_float, F2R. simpl.
+rewrite bpow_plus, Z2R_plus.
+rewrite Hq, Z2R_plus, Z2R_mult.
+replace ((Z2R m2 * Z2R q + Z2R r) * (bpow e' * bpow e2) / (Z2R m2 * bpow e2))%R
+  with ((Z2R q + Z2R r / Z2R m2) * bpow e')%R.
+apply inbetween_mult_compat.
+apply bpow_gt_0.
+destruct (Z_lt_le_dec 1 m2) as [Hm2''|Hm2''].
+replace (Z2R 1) with (Z2R m2 * /Z2R m2)%R.
+apply new_location_correct ; try easy.
+apply Rinv_0_lt_compat.
+now apply (Z2R_lt 0).
+now constructor.
+apply Rinv_r.
+apply Rgt_not_eq.
+now apply (Z2R_lt 0).
+assert (r = 0 /\ m2 = 1)%Z by (clear -Hr Hm2'' ; omega).
+rewrite (proj1 H), (proj2 H).
+unfold Rdiv.
+rewrite Rmult_0_l, Rplus_0_r.
+now constructor.
+field.
+split ; apply Rgt_not_eq.
+apply bpow_gt_0.
+now apply (Z2R_lt 0).
+Qed.
+
+End Fcalc_div.
diff --git a/flocq/Calc/Fcalc_ops.v b/flocq/Calc/Fcalc_ops.v
new file mode 100644
index 00000000..15bb2111
--- /dev/null
+++ b/flocq/Calc/Fcalc_ops.v
@@ -0,0 +1,161 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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.
+*)
+
+(** Basic operations on floats: alignment, addition, multiplication *)
+Require Import Fcore_Raux.
+Require Import Fcore_defs.
+Require Import Fcore_float_prop.
+
+Section Float_ops.
+
+Variable beta : radix.
+
+Notation bpow e := (bpow beta e).
+
+Definition Falign (f1 f2 : float beta) :=
+  let '(Float m1 e1) := f1 in
+  let '(Float m2 e2) := f2 in
+  if Zle_bool e1 e2
+  then (m1, (m2 * Zpower beta (e2 - e1))%Z, e1)
+  else ((m1 * Zpower beta (e1 - e2))%Z, m2, e2).
+
+Theorem Falign_spec :
+  forall f1 f2 : float beta,
+  let '(m1, m2, e) := Falign f1 f2 in
+  F2R f1 = F2R (Float beta m1 e) /\ F2R f2 = F2R (Float beta m2 e).
+Proof.
+unfold Falign.
+intros (m1, e1) (m2, e2).
+generalize (Zle_cases e1 e2).
+case (Zle_bool e1 e2) ; intros He ; split ; trivial.
+now rewrite <- F2R_change_exp.
+rewrite <- F2R_change_exp.
+apply refl_equal.
+omega.
+Qed.
+
+Theorem Falign_spec_exp:
+  forall f1 f2 : float beta,
+  snd (Falign f1 f2) = Zmin (Fexp f1) (Fexp f2).
+Proof.
+intros (m1,e1) (m2,e2).
+unfold Falign; simpl.
+generalize (Zle_cases e1 e2);case (Zle_bool e1 e2); intros He.
+case (Zmin_spec e1 e2); intros (H1,H2); easy.
+case (Zmin_spec e1 e2); intros (H1,H2); easy.
+Qed.
+
+Definition Fopp (f1: float beta) :=
+  let '(Float m1 e1) := f1 in
+  Float beta (-m1)%Z e1.
+
+Theorem F2R_opp :
+  forall f1 : float beta,
+  (F2R (Fopp f1) = -F2R f1)%R.
+intros (m1,e1).
+apply F2R_Zopp.
+Qed.
+
+Definition Fabs (f1: float beta) :=
+  let '(Float m1 e1) := f1 in
+  Float beta (Zabs m1)%Z e1.
+
+Theorem F2R_abs :
+  forall f1 : float beta,
+  (F2R (Fabs f1) = Rabs (F2R f1))%R.
+intros (m1,e1).
+apply F2R_Zabs.
+Qed.
+
+Definition Fplus (f1 f2 : float beta) :=
+  let '(m1, m2 ,e) := Falign f1 f2 in
+  Float beta (m1 + m2) e.
+
+Theorem F2R_plus :
+  forall f1 f2 : float beta,
+  F2R (Fplus f1 f2) = (F2R f1 + F2R f2)%R.
+Proof.
+intros f1 f2.
+unfold Fplus.
+generalize (Falign_spec f1 f2).
+destruct (Falign f1 f2) as ((m1, m2), e).
+intros (H1, H2).
+rewrite H1, H2.
+unfold F2R. simpl.
+rewrite Z2R_plus.
+apply Rmult_plus_distr_r.
+Qed.
+
+Theorem Fplus_same_exp :
+  forall m1 m2 e,
+  Fplus (Float beta m1 e) (Float beta m2 e) = Float beta (m1 + m2) e.
+Proof.
+intros m1 m2 e.
+unfold Fplus.
+simpl.
+now rewrite Zle_bool_refl, Zminus_diag, Zmult_1_r.
+Qed.
+
+Theorem Fexp_Fplus :
+  forall f1 f2 : float beta,
+  Fexp (Fplus f1 f2) = Zmin (Fexp f1) (Fexp f2).
+Proof.
+intros f1 f2.
+unfold Fplus.
+rewrite <- Falign_spec_exp.
+now destruct (Falign f1 f2) as ((p,q),e).
+Qed.
+
+Definition Fminus (f1 f2 : float beta) :=
+  Fplus f1 (Fopp f2).
+
+Theorem F2R_minus :
+  forall f1 f2 : float beta,
+  F2R (Fminus f1 f2) = (F2R f1 - F2R f2)%R.
+Proof.
+intros f1 f2; unfold Fminus.
+rewrite F2R_plus, F2R_opp.
+ring.
+Qed.
+
+Theorem Fminus_same_exp :
+  forall m1 m2 e,
+  Fminus (Float beta m1 e) (Float beta m2 e) = Float beta (m1 - m2) e.
+Proof.
+intros m1 m2 e.
+unfold Fminus.
+apply Fplus_same_exp.
+Qed.
+
+Definition Fmult (f1 f2 : float beta) :=
+  let '(Float m1 e1) := f1 in
+  let '(Float m2 e2) := f2 in
+  Float beta (m1 * m2) (e1 + e2).
+
+Theorem F2R_mult :
+  forall f1 f2 : float beta,
+  F2R (Fmult f1 f2) = (F2R f1 * F2R f2)%R.
+Proof.
+intros (m1, e1) (m2, e2).
+unfold Fmult, F2R. simpl.
+rewrite Z2R_mult, bpow_plus.
+ring.
+Qed.
+
+End Float_ops.
diff --git a/flocq/Calc/Fcalc_round.v b/flocq/Calc/Fcalc_round.v
new file mode 100644
index 00000000..3d31aea7
--- /dev/null
+++ b/flocq/Calc/Fcalc_round.v
@@ -0,0 +1,1093 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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 for computing the rounded value of a real number. *)
+
+Require Import Fcore.
+Require Import Fcore_digits.
+Require Import Fcalc_bracket.
+Require Import Fcalc_digits.
+
+Section Fcalc_round.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+Section Fcalc_round_fexp.
+
+Variable fexp : Z -> Z.
+Context { valid_exp : Valid_exp fexp }.
+Notation format := (generic_format beta fexp).
+
+(** Relates location and rounding. *)
+
+Theorem inbetween_float_round :
+  forall rnd choice,
+  ( forall x m l, inbetween_int m x l -> rnd x = choice m l ) ->
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp rnd x = F2R (Float beta (choice m l) e).
+Proof.
+intros rnd choice Hc x m l e Hl.
+unfold round, F2R. simpl.
+apply (f_equal (fun m => (Z2R m * bpow e)%R)).
+apply Hc.
+apply inbetween_mult_reg with (bpow e).
+apply bpow_gt_0.
+now rewrite scaled_mantissa_mult_bpow.
+Qed.
+
+Definition cond_incr (b : bool) m := if b then (m + 1)%Z else m.
+
+Theorem inbetween_float_round_sign :
+  forall rnd choice,
+  ( forall x m l, inbetween_int m (Rabs x) l ->
+    rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l) ) ->
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e).
+Proof.
+intros rnd choice Hc x m l e Hx.
+apply (f_equal (fun m => (Z2R m * bpow e)%R)).
+simpl.
+replace (Rlt_bool x 0) with (Rlt_bool (scaled_mantissa beta fexp x) 0).
+(* *)
+apply Hc.
+apply inbetween_mult_reg with (bpow e).
+apply bpow_gt_0.
+rewrite <- (Rabs_right (bpow e)) at 3.
+rewrite <- Rabs_mult.
+now rewrite scaled_mantissa_mult_bpow.
+apply Rle_ge.
+apply bpow_ge_0.
+(* *)
+destruct (Rlt_bool_spec x 0) as [Zx|Zx] ; simpl.
+apply Rlt_bool_true.
+rewrite <- (Rmult_0_l (bpow (-e))).
+apply Rmult_lt_compat_r with (2 := Zx).
+apply bpow_gt_0.
+apply Rlt_bool_false.
+apply Rmult_le_pos with (1 := Zx).
+apply bpow_ge_0.
+Qed.
+
+(** Relates location and rounding down. *)
+
+Theorem inbetween_int_DN :
+  forall x m l,
+  inbetween_int m x l ->
+  Zfloor x = m.
+Proof.
+intros x m l Hl.
+refine (Zfloor_imp m _ _).
+apply inbetween_bounds with (2 := Hl).
+apply Z2R_lt.
+apply Zlt_succ.
+Qed.
+
+Theorem inbetween_float_DN :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp Zfloor x = F2R (Float beta m e).
+Proof.
+apply inbetween_float_round with (choice := fun m l => m).
+exact inbetween_int_DN.
+Qed.
+
+Definition round_sign_DN s l :=
+  match l with
+  | loc_Exact => false
+  | _ => s
+  end.
+
+Theorem inbetween_int_DN_sign :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  Zfloor x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m).
+Proof.
+intros x m l Hl.
+unfold Rabs in Hl.
+destruct (Rcase_abs x) as [Zx|Zx] .
+(* *)
+rewrite Rlt_bool_true with (1 := Zx).
+inversion_clear Hl ; simpl.
+rewrite <- (Ropp_involutive x).
+rewrite H, <- Z2R_opp.
+apply Zfloor_Z2R.
+apply Zfloor_imp.
+split.
+apply Rlt_le.
+rewrite Z2R_opp.
+apply Ropp_lt_cancel.
+now rewrite Ropp_involutive.
+ring_simplify (- (m + 1) + 1)%Z.
+rewrite Z2R_opp.
+apply Ropp_lt_cancel.
+now rewrite Ropp_involutive.
+(* *)
+rewrite Rlt_bool_false.
+inversion_clear Hl ; simpl.
+rewrite H.
+apply Zfloor_Z2R.
+apply Zfloor_imp.
+split.
+now apply Rlt_le.
+apply H.
+now apply Rge_le.
+Qed.
+
+Theorem inbetween_float_DN_sign :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp Zfloor x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_DN (Rlt_bool x 0) l) m)) e).
+Proof.
+apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_DN s l) m).
+exact inbetween_int_DN_sign.
+Qed.
+
+(** Relates location and rounding up. *)
+
+Definition round_UP l :=
+  match l with
+  | loc_Exact => false
+  | _ => true
+  end.
+
+Theorem inbetween_int_UP :
+  forall x m l,
+  inbetween_int m x l ->
+  Zceil x = cond_incr (round_UP l) m.
+Proof.
+intros x m l Hl.
+assert (Hl': l = loc_Exact \/ (l <> loc_Exact /\ round_UP l = true)).
+case l ; try (now left) ; now right ; split.
+destruct Hl' as [Hl'|(Hl1, Hl2)].
+(* Exact *)
+rewrite Hl'.
+destruct Hl ; try easy.
+rewrite H.
+exact (Zceil_Z2R _).
+(* not Exact *)
+rewrite Hl2.
+simpl.
+apply Zceil_imp.
+ring_simplify (m + 1 - 1)%Z.
+refine (let H := _ in conj (proj1 H) (Rlt_le _ _ (proj2 H))).
+apply inbetween_bounds_not_Eq with (2 := Hl1) (1 := Hl).
+Qed.
+
+Theorem inbetween_float_UP :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp Zceil x = F2R (Float beta (cond_incr (round_UP l) m) e).
+Proof.
+apply inbetween_float_round with (choice := fun m l => cond_incr (round_UP l) m).
+exact inbetween_int_UP.
+Qed.
+
+Definition round_sign_UP s l :=
+  match l with
+  | loc_Exact => false
+  | _ => negb s
+  end.
+
+Theorem inbetween_int_UP_sign :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  Zceil x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m).
+Proof.
+intros x m l Hl.
+unfold Rabs in Hl.
+destruct (Rcase_abs x) as [Zx|Zx] .
+(* *)
+rewrite Rlt_bool_true with (1 := Zx).
+simpl.
+unfold Zceil.
+apply f_equal.
+inversion_clear Hl ; simpl.
+rewrite H.
+apply Zfloor_Z2R.
+apply Zfloor_imp.
+split.
+now apply Rlt_le.
+apply H.
+(* *)
+rewrite Rlt_bool_false.
+simpl.
+inversion_clear Hl ; simpl.
+rewrite H.
+apply Zceil_Z2R.
+apply Zceil_imp.
+split.
+change (m + 1 - 1)%Z with (Zpred (Zsucc m)).
+now rewrite <- Zpred_succ.
+now apply Rlt_le.
+now apply Rge_le.
+Qed.
+
+Theorem inbetween_float_UP_sign :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp Zceil x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_sign_UP (Rlt_bool x 0) l) m)) e).
+Proof.
+apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_sign_UP s l) m).
+exact inbetween_int_UP_sign.
+Qed.
+
+(** Relates location and rounding toward zero. *)
+
+Definition round_ZR (s : bool) l :=
+  match l with
+  | loc_Exact => false
+  | _ => s
+  end.
+
+Theorem inbetween_int_ZR :
+  forall x m l,
+  inbetween_int m x l ->
+  Ztrunc x = cond_incr (round_ZR (Zlt_bool m 0) l) m.
+Proof with auto with typeclass_instances.
+intros x m l Hl.
+inversion_clear Hl as [Hx|l' Hx Hl'].
+(* Exact *)
+rewrite Hx.
+rewrite Zrnd_Z2R...
+(* not Exact *)
+unfold Ztrunc.
+assert (Hm: Zfloor x = m).
+apply Zfloor_imp.
+exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)).
+rewrite Zceil_floor_neq.
+rewrite Hm.
+unfold cond_incr.
+simpl.
+case Rlt_bool_spec ; intros Hx' ;
+  case Zlt_bool_spec ; intros Hm' ; try apply refl_equal.
+elim Rlt_not_le with (1 := Hx').
+apply Rlt_le.
+apply Rle_lt_trans with (2 := proj1 Hx).
+now apply (Z2R_le 0).
+elim Rle_not_lt with (1 := Hx').
+apply Rlt_le_trans with (1 := proj2 Hx).
+apply (Z2R_le _ 0).
+now apply Zlt_le_succ.
+rewrite Hm.
+now apply Rlt_not_eq.
+Qed.
+
+Theorem inbetween_float_ZR :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp Ztrunc x = F2R (Float beta (cond_incr (round_ZR (Zlt_bool m 0) l) m) e).
+Proof.
+apply inbetween_float_round with (choice := fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m).
+exact inbetween_int_ZR.
+Qed.
+
+Theorem inbetween_int_ZR_sign :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  Ztrunc x = cond_Zopp (Rlt_bool x 0) m.
+Proof.
+intros x m l Hl.
+simpl.
+unfold Ztrunc.
+destruct (Rlt_le_dec x 0) as [Zx|Zx].
+(* *)
+rewrite Rlt_bool_true with (1 := Zx).
+simpl.
+unfold Zceil.
+apply f_equal.
+apply Zfloor_imp.
+rewrite <- Rabs_left with (1 := Zx).
+apply inbetween_bounds with (2 := Hl).
+apply Z2R_lt.
+apply Zlt_succ.
+(* *)
+rewrite Rlt_bool_false with (1 := Zx).
+simpl.
+apply Zfloor_imp.
+rewrite <- Rabs_pos_eq with (1 := Zx).
+apply inbetween_bounds with (2 := Hl).
+apply Z2R_lt.
+apply Zlt_succ.
+Qed.
+
+Theorem inbetween_float_ZR_sign :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp Ztrunc x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) m) e).
+Proof.
+apply inbetween_float_round_sign with (choice := fun s m l => m).
+exact inbetween_int_ZR_sign.
+Qed.
+
+(** Relates location and rounding to nearest. *)
+
+Definition round_N (p : bool) l :=
+  match l with
+  | loc_Exact => false
+  | loc_Inexact Lt => false
+  | loc_Inexact Eq => p
+  | loc_Inexact Gt => true
+  end.
+
+Theorem inbetween_int_N :
+  forall choice x m l,
+  inbetween_int m x l ->
+  Znearest choice x = cond_incr (round_N (choice m) l) m.
+Proof with auto with typeclass_instances.
+intros choice x m l Hl.
+inversion_clear Hl as [Hx|l' Hx Hl'].
+(* Exact *)
+rewrite Hx.
+rewrite Zrnd_Z2R...
+(* not Exact *)
+unfold Znearest.
+assert (Hm: Zfloor x = m).
+apply Zfloor_imp.
+exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)).
+rewrite Zceil_floor_neq.
+rewrite Hm.
+replace (Rcompare (x - Z2R m) (/2)) with l'.
+now case l'.
+rewrite <- Hl'.
+rewrite Z2R_plus.
+rewrite <- (Rcompare_plus_r (- Z2R m) x).
+apply f_equal.
+simpl (Z2R 1).
+field.
+rewrite Hm.
+now apply Rlt_not_eq.
+Qed.
+
+Theorem inbetween_int_N_sign :
+  forall choice x m l,
+  inbetween_int m (Rabs x) l ->
+  Znearest choice x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (if Rlt_bool x 0 then negb (choice (-(m + 1))%Z) else choice m) l) m).
+Proof with auto with typeclass_instances.
+intros choice x m l Hl.
+simpl.
+unfold Rabs in Hl.
+destruct (Rcase_abs x) as [Zx|Zx].
+(* *)
+rewrite Rlt_bool_true with (1 := Zx).
+simpl.
+rewrite <- (Ropp_involutive x).
+rewrite Znearest_opp.
+apply f_equal.
+inversion_clear Hl as [Hx|l' Hx Hl'].
+rewrite Hx.
+apply Zrnd_Z2R...
+assert (Hm: Zfloor (-x) = m).
+apply Zfloor_imp.
+exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)).
+unfold Znearest.
+rewrite Zceil_floor_neq.
+rewrite Hm.
+replace (Rcompare (- x - Z2R m) (/2)) with l'.
+now case l'.
+rewrite <- Hl'.
+rewrite Z2R_plus.
+rewrite <- (Rcompare_plus_r (- Z2R m) (-x)).
+apply f_equal.
+simpl (Z2R 1).
+field.
+rewrite Hm.
+now apply Rlt_not_eq.
+(* *)
+generalize (Rge_le _ _ Zx).
+clear Zx. intros Zx.
+rewrite Rlt_bool_false with (1 := Zx).
+simpl.
+inversion_clear Hl as [Hx|l' Hx Hl'].
+rewrite Hx.
+apply Zrnd_Z2R...
+assert (Hm: Zfloor x = m).
+apply Zfloor_imp.
+exact (conj (Rlt_le _ _ (proj1 Hx)) (proj2 Hx)).
+unfold Znearest.
+rewrite Zceil_floor_neq.
+rewrite Hm.
+replace (Rcompare (x - Z2R m) (/2)) with l'.
+now case l'.
+rewrite <- Hl'.
+rewrite Z2R_plus.
+rewrite <- (Rcompare_plus_r (- Z2R m) x).
+apply f_equal.
+simpl (Z2R 1).
+field.
+rewrite Hm.
+now apply Rlt_not_eq.
+Qed.
+
+(** Relates location and rounding to nearest even. *)
+
+Theorem inbetween_int_NE :
+  forall x m l,
+  inbetween_int m x l ->
+  ZnearestE x = cond_incr (round_N (negb (Zeven m)) l) m.
+Proof.
+intros x m l Hl.
+now apply inbetween_int_N with (choice := fun x => negb (Zeven x)).
+Qed.
+
+Theorem inbetween_float_NE :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp ZnearestE x = F2R (Float beta (cond_incr (round_N (negb (Zeven m)) l) m) e).
+Proof.
+apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (negb (Zeven m)) l) m).
+exact inbetween_int_NE.
+Qed.
+
+Theorem inbetween_int_NE_sign :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  ZnearestE x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Zeven m)) l) m).
+Proof.
+intros x m l Hl.
+erewrite inbetween_int_N_sign with (choice := fun x => negb (Zeven x)).
+2: eexact Hl.
+apply f_equal.
+case Rlt_bool.
+rewrite Zeven_opp, Zeven_plus.
+now case (Zeven m).
+apply refl_equal.
+Qed.
+
+Theorem inbetween_float_NE_sign :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e (Rabs x) l ->
+  round beta fexp ZnearestE x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (cond_incr (round_N (negb (Zeven m)) l) m)) e).
+Proof.
+apply inbetween_float_round_sign with (choice := fun s m l => cond_incr (round_N (negb (Zeven m)) l) m).
+exact inbetween_int_NE_sign.
+Qed.
+
+(** Relates location and rounding to nearest away. *)
+
+Theorem inbetween_int_NA :
+  forall x m l,
+  inbetween_int m x l ->
+  ZnearestA x = cond_incr (round_N (Zle_bool 0 m) l) m.
+Proof.
+intros x m l Hl.
+now apply inbetween_int_N with (choice := fun x => Zle_bool 0 x).
+Qed.
+
+Theorem inbetween_float_NA :
+  forall x m l,
+  let e := canonic_exp beta fexp x in
+  inbetween_float beta m e x l ->
+  round beta fexp ZnearestA x = F2R (Float beta (cond_incr (round_N (Zle_bool 0 m) l) m) e).
+Proof.
+apply inbetween_float_round with (choice := fun m l => cond_incr (round_N (Zle_bool 0 m) l) m).
+exact inbetween_int_NA.
+Qed.
+
+Theorem inbetween_int_NA_sign :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  ZnearestA x = cond_Zopp (Rlt_bool x 0) (cond_incr (round_N true l) m).
+Proof.
+intros x m l Hl.
+erewrite inbetween_int_N_sign with (choice := Zle_bool 0).
+2: eexact Hl.
+apply f_equal.
+assert (Hm: (0 <= m)%Z).
+apply Zlt_succ_le.
+apply lt_Z2R.
+apply Rle_lt_trans with (Rabs x).
+apply Rabs_pos.
+refine (proj2 (inbetween_bounds _ _ _ _ _ Hl)).
+apply Z2R_lt.
+apply Zlt_succ.
+rewrite Zle_bool_true with (1 := Hm).
+rewrite Zle_bool_false.
+now case Rlt_bool.
+omega.
+Qed.
+
+Definition truncate_aux t k :=
+  let '(m, e, l) := t in
+  let p := Zpower beta k in
+  (Zdiv m p, (e + k)%Z, new_location p (Zmod m p) l).
+
+Theorem truncate_aux_comp :
+  forall t k1 k2,
+  (0 < k1)%Z ->
+  (0 < k2)%Z ->
+  truncate_aux t (k1 + k2) = truncate_aux (truncate_aux t k1) k2.
+Proof.
+intros ((m,e),l) k1 k2 Hk1 Hk2.
+unfold truncate_aux.
+destruct (inbetween_float_ex beta m e l) as (x,Hx).
+assert (B1 := inbetween_float_new_location _ _ _ _ _ _ Hk1 Hx).
+assert (Hk3: (0 < k1 + k2)%Z).
+change Z0 with (0 + 0)%Z.
+now apply Zplus_lt_compat.
+assert (B3 := inbetween_float_new_location _ _ _ _ _ _ Hk3 Hx).
+assert (B2 := inbetween_float_new_location _ _ _ _ _ _ Hk2 B1).
+rewrite Zplus_assoc in B3.
+destruct (inbetween_float_unique _ _ _ _ _ _ _ B2 B3).
+now rewrite H, H0, Zplus_assoc.
+Qed.
+
+(** Given a triple (mantissa, exponent, position), this function
+    computes a triple with a canonic exponent, assuming the
+    original triple had enough precision. *)
+
+Definition truncate t :=
+  let '(m, e, l) := t in
+  let k := (fexp (Zdigits beta m + e) - e)%Z in
+  if Zlt_bool 0 k then truncate_aux t k
+  else t.
+
+Theorem truncate_0 :
+  forall e l,
+  let '(m', e', l') := truncate (0, e, l)%Z in
+  m' = Z0.
+Proof.
+intros e l.
+unfold truncate.
+case Zlt_bool.
+simpl.
+apply Zdiv_0_l.
+apply refl_equal.
+Qed.
+
+Theorem generic_format_truncate :
+  forall m e l,
+  (0 <= m)%Z ->
+  let '(m', e', l') := truncate (m, e, l) in
+  format (F2R (Float beta m' e')).
+Proof.
+intros m e l Hm.
+unfold truncate.
+set (k := (fexp (Zdigits beta m + e) - e)%Z).
+case Zlt_bool_spec ; intros Hk.
+(* *)
+unfold truncate_aux.
+apply generic_format_F2R.
+intros Hd.
+unfold canonic_exp.
+rewrite ln_beta_F2R_Zdigits with (1 := Hd).
+rewrite Zdigits_div_Zpower with (1 := Hm).
+replace (Zdigits beta m - k + (e + k))%Z with (Zdigits beta m + e)%Z by ring.
+unfold k.
+ring_simplify.
+apply Zle_refl.
+split.
+now apply Zlt_le_weak.
+apply Znot_gt_le.
+contradict Hd.
+apply Zdiv_small.
+apply conj with (1 := Hm).
+rewrite <- Zabs_eq with (1 := Hm).
+apply Zpower_gt_Zdigits.
+apply Zlt_le_weak.
+now apply Zgt_lt.
+(* *)
+destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|Hm'].
+apply generic_format_F2R.
+unfold canonic_exp.
+rewrite ln_beta_F2R_Zdigits.
+2: now apply Zgt_not_eq.
+unfold k in Hk. clear -Hk.
+omega.
+rewrite <- Hm', F2R_0.
+apply generic_format_0.
+Qed.
+
+Theorem truncate_correct_format :
+  forall m e,
+  m <> Z0 ->
+  let x := F2R (Float beta m e) in
+  generic_format beta fexp x ->
+  (e <= fexp (Zdigits beta m + e))%Z ->
+  let '(m', e', l') := truncate (m, e, loc_Exact) in
+  x = F2R (Float beta m' e') /\ e' = canonic_exp beta fexp x.
+Proof.
+intros m e Hm x Fx He.
+assert (Hc: canonic_exp beta fexp x = fexp (Zdigits beta m + e)).
+unfold canonic_exp, x.
+now rewrite ln_beta_F2R_Zdigits.
+unfold truncate.
+rewrite <- Hc.
+set (k := (canonic_exp beta fexp x - e)%Z).
+case Zlt_bool_spec ; intros Hk.
+(* *)
+unfold truncate_aux.
+rewrite Fx at 1.
+refine (let H := _ in conj _ H).
+unfold k. ring.
+rewrite <- H.
+apply F2R_eq_compat.
+replace (scaled_mantissa beta fexp x) with (Z2R (Zfloor (scaled_mantissa beta fexp x))).
+rewrite Ztrunc_Z2R.
+unfold scaled_mantissa.
+rewrite <- H.
+unfold x, F2R. simpl.
+rewrite Rmult_assoc, <- bpow_plus.
+ring_simplify (e + - (e + k))%Z.
+clear -Hm Hk.
+destruct k as [|k|k] ; try easy.
+simpl.
+apply Zfloor_div.
+intros H.
+generalize (Zpower_pos_gt_0 beta k) (Zle_bool_imp_le _ _ (radix_prop beta)).
+omega.
+rewrite scaled_mantissa_generic with (1 := Fx).
+now rewrite Zfloor_Z2R.
+(* *)
+split.
+apply refl_equal.
+unfold k in Hk.
+omega.
+Qed.
+
+Theorem truncate_correct_partial :
+  forall x m e l,
+  (0 < x)%R ->
+  inbetween_float beta m e x l ->
+  (e <= fexp (Zdigits beta m + e))%Z ->
+  let '(m', e', l') := truncate (m, e, l) in
+  inbetween_float beta m' e' x l' /\ e' = canonic_exp beta fexp x.
+Proof.
+intros x m e l Hx H1 H2.
+unfold truncate.
+set (k := (fexp (Zdigits beta m + e) - e)%Z).
+set (p := Zpower beta k).
+(* *)
+assert (Hx': (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R).
+apply inbetween_float_bounds with (1 := H1).
+(* *)
+assert (Hm: (0 <= m)%Z).
+cut (0 < m + 1)%Z. omega.
+apply F2R_lt_reg with beta e.
+rewrite F2R_0.
+apply Rlt_trans with  (1 := Hx).
+apply Hx'.
+assert (He: (e + k = canonic_exp beta fexp x)%Z).
+(* . *)
+unfold canonic_exp.
+destruct (Zle_lt_or_eq _ _ Hm) as [Hm'|Hm'].
+(* .. 0 < m *)
+rewrite ln_beta_F2R_bounds with (1 := Hm') (2 := Hx').
+assert (H: m <> Z0).
+apply sym_not_eq.
+now apply Zlt_not_eq.
+rewrite ln_beta_F2R with (1 := H).
+rewrite <- Zdigits_ln_beta with (1 := H).
+unfold k.
+ring.
+(* .. m = 0 *)
+rewrite <- Hm' in H2.
+destruct (ln_beta beta x) as (ex, Hex).
+simpl.
+specialize (Hex (Rgt_not_eq _ _ Hx)).
+unfold k.
+ring_simplify.
+rewrite <- Hm'.
+simpl.
+apply sym_eq.
+apply valid_exp.
+exact H2.
+apply Zle_trans with e.
+eapply bpow_lt_bpow.
+apply Rle_lt_trans with (1 := proj1 Hex).
+rewrite Rabs_pos_eq.
+rewrite <- F2R_bpow.
+rewrite <- Hm' in Hx'.
+apply Hx'.
+now apply Rlt_le.
+exact H2.
+(* . *)
+generalize (Zlt_cases 0 k).
+case (Zlt_bool 0 k) ; intros Hk ; unfold k in Hk.
+split.
+now apply inbetween_float_new_location.
+exact He.
+split.
+exact H1.
+rewrite <- He.
+unfold k.
+omega.
+Qed.
+
+Theorem truncate_correct :
+  forall x m e l,
+  (0 <= x)%R ->
+  inbetween_float beta m e x l ->
+  (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact ->
+  let '(m', e', l') := truncate (m, e, l) in
+  inbetween_float beta m' e' x l' /\
+  (e' = canonic_exp beta fexp x \/ (l' = loc_Exact /\ format x)).
+Proof.
+intros x m e l [Hx|Hx] H1 H2.
+(* 0 < x *)
+destruct (Zle_or_lt e (fexp (Zdigits beta m + e))) as [H3|H3].
+(* . enough digits *)
+generalize (truncate_correct_partial x m e l Hx H1 H3).
+destruct (truncate (m, e, l)) as ((m', e'), l').
+intros (H4, H5).
+split.
+exact H4.
+now left.
+(* . not enough digits but loc_Exact *)
+destruct H2 as [H2|H2].
+elim (Zlt_irrefl e).
+now apply Zle_lt_trans with (1 := H2).
+rewrite H2 in H1 |- *.
+unfold truncate.
+generalize (Zlt_cases 0 (fexp (Zdigits beta m + e) - e)).
+case Zlt_bool.
+intros H.
+apply False_ind.
+omega.
+intros _.
+split.
+exact H1.
+right.
+split.
+apply refl_equal.
+inversion_clear H1.
+rewrite H.
+apply generic_format_F2R.
+intros Zm.
+unfold canonic_exp.
+rewrite ln_beta_F2R_Zdigits with (1 := Zm).
+now apply Zlt_le_weak.
+(* x = 0 *)
+assert (Hm: m = Z0).
+cut (m <= 0 < m + 1)%Z. omega.
+assert (Hx': (F2R (Float beta m e) <= x < F2R (Float beta (m + 1) e))%R).
+apply inbetween_float_bounds with (1 := H1).
+rewrite <- Hx in Hx'.
+split.
+apply F2R_le_0_reg with (1 := proj1 Hx').
+apply F2R_gt_0_reg with (1 := proj2 Hx').
+rewrite Hm, <- Hx in H1 |- *.
+clear -H1.
+case H1.
+(* . *)
+intros _.
+assert (exists e', truncate (Z0, e, loc_Exact) = (Z0, e', loc_Exact)).
+unfold truncate, truncate_aux.
+case Zlt_bool.
+rewrite Zdiv_0_l, Zmod_0_l.
+eexists.
+apply f_equal.
+unfold new_location.
+now case Zeven.
+now eexists.
+destruct H as (e', H).
+rewrite H.
+split.
+constructor.
+apply sym_eq.
+apply F2R_0.
+right.
+repeat split.
+apply generic_format_0.
+(* . *)
+intros l' (H, _) _.
+rewrite F2R_0 in H.
+elim Rlt_irrefl with (1 := H).
+Qed.
+
+Section round_dir.
+
+Variable rnd : R -> Z.
+Context { valid_rnd : Valid_rnd rnd }.
+
+Variable choice : Z -> location -> Z.
+Hypothesis inbetween_int_valid :
+  forall x m l,
+  inbetween_int m x l ->
+  rnd x = choice m l.
+
+Theorem round_any_correct :
+  forall x m e l,
+  inbetween_float beta m e x l ->
+  (e = canonic_exp beta fexp x \/ (l = loc_Exact /\ format x)) ->
+  round beta fexp rnd x = F2R (Float beta (choice m l) e).
+Proof with auto with typeclass_instances.
+intros x m e l Hin [He|(Hl,Hf)].
+rewrite He in Hin |- *.
+apply inbetween_float_round with (2 := Hin).
+exact inbetween_int_valid.
+rewrite Hl in Hin.
+inversion_clear Hin.
+rewrite Hl.
+replace (choice m loc_Exact) with m.
+rewrite <- H.
+apply round_generic...
+rewrite <- (Zrnd_Z2R rnd m) at 1.
+apply inbetween_int_valid.
+now constructor.
+Qed.
+
+(** Truncating a triple is sufficient to round a real number. *)
+
+Theorem round_trunc_any_correct :
+  forall x m e l,
+  (0 <= x)%R ->
+  inbetween_float beta m e x l ->
+  (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact ->
+  round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (choice m' l') e').
+Proof.
+intros x m e l Hx Hl He.
+generalize (truncate_correct x m e l Hx Hl He).
+destruct (truncate (m, e, l)) as ((m', e'), l').
+intros (H1, H2).
+now apply round_any_correct.
+Qed.
+
+End round_dir.
+
+Section round_dir_sign.
+
+Variable rnd : R -> Z.
+Context { valid_rnd : Valid_rnd rnd }.
+
+Variable choice : bool -> Z -> location -> Z.
+Hypothesis inbetween_int_valid :
+  forall x m l,
+  inbetween_int m (Rabs x) l ->
+  rnd x = cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l).
+
+Theorem round_sign_any_correct :
+  forall x m e l,
+  inbetween_float beta m e (Rabs x) l ->
+  (e = canonic_exp beta fexp x \/ (l = loc_Exact /\ format x)) ->
+  round beta fexp rnd x = F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m l)) e).
+Proof with auto with typeclass_instances.
+intros x m e l Hin [He|(Hl,Hf)].
+rewrite He in Hin |- *.
+apply inbetween_float_round_sign with (2 := Hin).
+exact inbetween_int_valid.
+rewrite Hl in Hin.
+inversion_clear Hin.
+rewrite Hl.
+replace (choice (Rlt_bool x 0) m loc_Exact) with m.
+(* *)
+unfold Rabs in H.
+destruct (Rcase_abs x) as [Zx|Zx].
+rewrite Rlt_bool_true with (1 := Zx).
+simpl.
+rewrite F2R_Zopp.
+rewrite <- H, Ropp_involutive.
+apply round_generic...
+rewrite Rlt_bool_false.
+simpl.
+rewrite <- H.
+now apply round_generic.
+now apply Rge_le.
+(* *)
+destruct (Rlt_bool_spec x 0) as [Zx|Zx].
+(* . *)
+apply Zopp_inj.
+change (- m = cond_Zopp true (choice true m loc_Exact))%Z.
+rewrite <- (Zrnd_Z2R rnd (-m)) at 1.
+assert (Z2R (-m) < 0)%R.
+rewrite Z2R_opp.
+apply Ropp_lt_gt_0_contravar.
+apply (Z2R_lt 0).
+apply F2R_gt_0_reg with beta e.
+rewrite <- H.
+apply Rabs_pos_lt.
+now apply Rlt_not_eq.
+rewrite <- Rlt_bool_true with (1 := H0).
+apply inbetween_int_valid.
+constructor.
+rewrite Rabs_left with (1 := H0).
+rewrite Z2R_opp.
+apply Ropp_involutive.
+(* . *)
+change (m = cond_Zopp false (choice false m loc_Exact))%Z.
+rewrite <- (Zrnd_Z2R rnd m) at 1.
+assert (0 <= Z2R m)%R.
+apply (Z2R_le 0).
+apply F2R_ge_0_reg with beta e.
+rewrite <- H.
+apply Rabs_pos.
+rewrite <- Rlt_bool_false with (1 := H0).
+apply inbetween_int_valid.
+constructor.
+now apply Rabs_pos_eq.
+Qed.
+
+(** Truncating a triple is sufficient to round a real number. *)
+
+Theorem round_trunc_sign_any_correct :
+  forall x m e l,
+  inbetween_float beta m e (Rabs x) l ->
+  (e <= fexp (Zdigits beta m + e))%Z \/ l = loc_Exact ->
+  round beta fexp rnd x = let '(m', e', l') := truncate (m, e, l) in F2R (Float beta (cond_Zopp (Rlt_bool x 0) (choice (Rlt_bool x 0) m' l')) e').
+Proof.
+intros x m e l Hl He.
+generalize (truncate_correct (Rabs x) m e l (Rabs_pos _) Hl He).
+destruct (truncate (m, e, l)) as ((m', e'), l').
+intros (H1, H2).
+apply round_sign_any_correct.
+exact H1.
+destruct H2 as [H2|(H2,H3)].
+left.
+now rewrite <- canonic_exp_abs.
+right.
+split.
+exact H2.
+unfold Rabs in H3.
+destruct (Rcase_abs x) in H3.
+rewrite <- Ropp_involutive.
+now apply generic_format_opp.
+exact H3.
+Qed.
+
+End round_dir_sign.
+
+(** * Instances of the theorems above, for the usual rounding modes. *)
+
+Definition round_DN_correct :=
+  round_any_correct _ (fun m _ => m) inbetween_int_DN.
+
+Definition round_trunc_DN_correct :=
+  round_trunc_any_correct _ (fun m _ => m) inbetween_int_DN.
+
+Definition round_sign_DN_correct :=
+  round_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign.
+
+Definition round_trunc_sign_DN_correct :=
+  round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_DN s l) m) inbetween_int_DN_sign.
+
+Definition round_UP_correct :=
+  round_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP.
+
+Definition round_trunc_UP_correct :=
+  round_trunc_any_correct _ (fun m l => cond_incr (round_UP l) m) inbetween_int_UP.
+
+Definition round_sign_UP_correct :=
+  round_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign.
+
+Definition round_trunc_sign_UP_correct :=
+  round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_sign_UP s l) m) inbetween_int_UP_sign.
+
+Definition round_ZR_correct :=
+  round_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR.
+
+Definition round_trunc_ZR_correct :=
+  round_trunc_any_correct _ (fun m l => cond_incr (round_ZR (Zlt_bool m 0) l) m) inbetween_int_ZR.
+
+Definition round_sign_ZR_correct :=
+  round_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign.
+
+Definition round_trunc_sign_ZR_correct :=
+  round_trunc_sign_any_correct _ (fun s m l => m) inbetween_int_ZR_sign.
+
+Definition round_NE_correct :=
+  round_any_correct _ (fun m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE.
+
+Definition round_trunc_NE_correct :=
+  round_trunc_any_correct _ (fun m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE.
+
+Definition round_sign_NE_correct :=
+  round_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE_sign.
+
+Definition round_trunc_sign_NE_correct :=
+  round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N (negb (Zeven m)) l) m) inbetween_int_NE_sign.
+
+Definition round_NA_correct :=
+  round_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA.
+
+Definition round_trunc_NA_correct :=
+  round_trunc_any_correct _ (fun m l => cond_incr (round_N (Zle_bool 0 m) l) m) inbetween_int_NA.
+
+Definition round_sign_NA_correct :=
+  round_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign.
+
+Definition round_trunc_sign_NA_correct :=
+  round_trunc_sign_any_correct _ (fun s m l => cond_incr (round_N true l) m) inbetween_int_NA_sign.
+
+End Fcalc_round_fexp.
+
+(** Specialization of truncate for FIX formats. *)
+
+Variable emin : Z.
+
+Definition truncate_FIX t :=
+  let '(m, e, l) := t in
+  let k := (emin - e)%Z in
+  if Zlt_bool 0 k then
+    let p := Zpower beta k in
+    (Zdiv m p, (e + k)%Z, new_location p (Zmod m p) l)
+  else t.
+
+Theorem truncate_FIX_correct :
+  forall x m e l,
+  inbetween_float beta m e x l ->
+  (e <= emin)%Z \/ l = loc_Exact ->
+  let '(m', e', l') := truncate_FIX (m, e, l) in
+  inbetween_float beta m' e' x l' /\
+  (e' = canonic_exp beta (FIX_exp emin) x \/ (l' = loc_Exact /\ generic_format beta (FIX_exp emin) x)).
+Proof.
+intros x m e l H1 H2.
+unfold truncate_FIX.
+set (k := (emin - e)%Z).
+set (p := Zpower beta k).
+unfold canonic_exp, FIX_exp.
+generalize (Zlt_cases 0 k).
+case (Zlt_bool 0 k) ; intros Hk.
+(* shift *)
+split.
+now apply inbetween_float_new_location.
+clear H2.
+left.
+unfold k.
+ring.
+(* no shift *)
+split.
+exact H1.
+unfold k in Hk.
+destruct H2 as [H2|H2].
+left.
+omega.
+right.
+split.
+exact H2.
+rewrite H2 in H1.
+inversion_clear H1.
+rewrite H.
+apply generic_format_F2R.
+unfold canonic_exp.
+omega.
+Qed.
+
+End Fcalc_round.
diff --git a/flocq/Calc/Fcalc_sqrt.v b/flocq/Calc/Fcalc_sqrt.v
new file mode 100644
index 00000000..e6fe74b2
--- /dev/null
+++ b/flocq/Calc/Fcalc_sqrt.v
@@ -0,0 +1,346 @@
+(**
+This file is part of the Flocq formalization of floating-point
+arithmetic in Coq: http://flocq.gforge.inria.fr/
+
+Copyright (C) 2010-2011 Sylvie Boldo
+#<br />#
+Copyright (C) 2010-2011 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 functions and theorems for computing the rounded square root of a floating-point number. *)
+
+Require Import Fcore_Raux.
+Require Import Fcore_defs.
+Require Import Fcore_digits.
+Require Import Fcore_float_prop.
+Require Import Fcalc_bracket.
+Require Import Fcalc_digits.
+
+Section Fcalc_sqrt.
+
+Fixpoint Zsqrt_aux (p : positive) : Z * Z :=
+  match p with
+  | xH => (1, 0)%Z
+  | xO xH => (1, 1)%Z
+  | xI xH => (1, 2)%Z
+  | xO (xO p) =>
+    let (q, r) := Zsqrt_aux p in
+    let r' := (4 * r)%Z in
+    let d := (r' - (4 * q + 1))%Z in
+    if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z
+  | xO (xI p) =>
+    let (q, r) := Zsqrt_aux p in
+    let r' := (4 * r + 2)%Z in
+    let d := (r' - (4 * q + 1))%Z in
+    if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z
+  | xI (xO p) =>
+    let (q, r) := Zsqrt_aux p in
+    let r' := (4 * r + 1)%Z in
+    let d := (r' - (4 * q + 1))%Z in
+    if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z
+  | xI (xI p) =>
+    let (q, r) := Zsqrt_aux p in
+    let r' := (4 * r + 3)%Z in
+    let d := (r' - (4 * q + 1))%Z in
+    if Zlt_bool d 0 then (2 * q, r')%Z else (2 * q + 1, d)%Z
+  end.
+
+Lemma Zsqrt_ind :
+  forall P : positive -> Prop,
+  P xH -> P (xO xH) -> P (xI xH) ->
+  ( forall p, P p -> P (xO (xO p)) /\ P (xO (xI p)) /\ P (xI (xO p)) /\ P (xI (xI p)) ) ->
+  forall p, P p.
+Proof.
+intros P H1 H2 H3 Hp.
+fix 1.
+intros [[p|p|]|[p|p|]|].
+refine (proj2 (proj2 (proj2 (Hp p _)))).
+apply Zsqrt_ind.
+refine (proj1 (proj2 (proj2 (Hp p _)))).
+apply Zsqrt_ind.
+exact H3.
+refine (proj1 (proj2 (Hp p _))).
+apply Zsqrt_ind.
+refine (proj1 (Hp p _)).
+apply Zsqrt_ind.
+exact H2.
+exact H1.
+Qed.
+
+Lemma Zsqrt_aux_correct :
+  forall p,
+  let (q, r) := Zsqrt_aux p in
+  Zpos p = (q * q + r)%Z /\ (0 <= r <= 2 * q)%Z.
+Proof.
+intros p.
+elim p using Zsqrt_ind ; clear p.
+now repeat split.
+now repeat split.
+now repeat split.
+intros p.
+Opaque Zmult. simpl. Transparent Zmult.
+destruct (Zsqrt_aux p) as (q, r).
+intros (Hq, Hr).
+change (Zpos p~0~0) with (4 * Zpos p)%Z.
+change (Zpos p~0~1) with (4 * Zpos p + 1)%Z.
+change (Zpos p~1~0) with (4 * Zpos p + 2)%Z.
+change (Zpos p~1~1) with (4 * Zpos p + 3)%Z.
+rewrite Hq. clear Hq.
+repeat split.
+generalize (Zlt_cases (4 * r - (4 * q + 1)) 0).
+case Zlt_bool ; ( split ; [ ring | omega ] ).
+generalize (Zlt_cases (4 * r + 2 - (4 * q + 1)) 0).
+case Zlt_bool ; ( split ; [ ring | omega ] ).
+generalize (Zlt_cases (4 * r + 1 - (4 * q + 1)) 0).
+case Zlt_bool ; ( split ; [ ring | omega ] ).
+generalize (Zlt_cases (4 * r + 3 - (4 * q + 1)) 0).
+case Zlt_bool ; ( split ; [ ring | omega ] ).
+Qed.
+
+(** Computes the integer square root and its remainder, but
+    without carrying a proof, contrarily to the operation of
+    the standard libary. *)
+
+Definition Zsqrt p :=
+  match p with
+  | Zpos p => Zsqrt_aux p
+  | _ => (0, 0)%Z
+  end.
+
+Theorem Zsqrt_correct :
+  forall x,
+  (0 <= x)%Z ->
+  let (q, r) := Zsqrt x in
+  x = (q * q + r)%Z /\ (0 <= r <= 2 * q)%Z.
+Proof.
+unfold Zsqrt.
+intros [|p|p] Hx.
+now repeat split.
+apply Zsqrt_aux_correct.
+now elim Hx.
+Qed.
+
+Variable beta : radix.
+Notation bpow e := (bpow beta e).
+
+(** Computes a mantissa of precision p, the corresponding exponent,
+    and the position with respect to the real square root of the
+    input floating-point number.
+
+The algorithm performs the following steps:
+- Shift the mantissa so that it has at least 2p-1 digits;
+  shift it one digit more if the new exponent is not even.
+- Compute the square root s (at least p digits) of the new
+  mantissa, and its remainder r.
+- Compute the position according to the remainder:
+  -- r == 0  =>  Eq,
+  -- r <= s  =>  Lo,
+  -- r >= s  =>  Up.
+
+Complexity is fine as long as p1 <= 2p-1.
+*)
+
+Definition Fsqrt_core prec m e :=
+  let d := Zdigits beta m in
+  let s := Zmax (2 * prec - d) 0 in
+  let e' := (e - s)%Z in
+  let (s', e'') := if Zeven e' then (s, e') else (s + 1, e' - 1)%Z in
+  let m' :=
+    match s' with
+    | Zpos p => (m * Zpower_pos beta p)%Z
+    | _ => m
+    end in
+  let (q, r) := Zsqrt 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, Zdiv2 e'', l).
+
+Theorem Fsqrt_core_correct :
+  forall prec m e,
+  (0 < m)%Z ->
+  let '(m', e', l) := Fsqrt_core prec m e in
+  (prec <= Zdigits beta m')%Z /\
+  inbetween_float beta m' e' (sqrt (F2R (Float beta m e))) l.
+Proof.
+intros prec m e Hm.
+unfold Fsqrt_core.
+set (d := Zdigits beta m).
+set (s := Zmax (2 * prec - d) 0).
+(* . exponent *)
+case_eq (if Zeven (e - s) then (s, (e - s)%Z) else ((s + 1)%Z, (e - s - 1)%Z)).
+intros s' e' Hse.
+assert (He: (Zeven e' = true /\ 0 <= s' /\ 2 * prec - d <= s' /\ s' + e' = e)%Z).
+revert Hse.
+case_eq (Zeven (e - s)) ; intros He Hse ; inversion Hse.
+repeat split.
+exact He.
+unfold s.
+apply Zle_max_r.
+apply Zle_max_l.
+ring.
+assert (H: (Zmax (2 * prec - d) 0 <= s + 1)%Z).
+fold s.
+apply Zle_succ.
+repeat split.
+unfold Zminus at 1.
+now rewrite Zeven_plus, He.
+apply Zle_trans with (2 := H).
+apply Zle_max_r.
+apply Zle_trans with (2 := H).
+apply Zle_max_l.
+ring.
+clear -Hm He.
+destruct He as (He1, (He2, (He3, He4))).
+(* . shift *)
+set (m' := match s' with
+  | Z0 => m
+  | Zpos p => (m * Zpower_pos beta p)%Z
+  | Zneg _ => m
+  end).
+assert (Hs: F2R (Float beta m' e') = F2R (Float beta m e) /\ (2 * prec <= Zdigits beta m')%Z /\ (0 < m')%Z).
+rewrite <- He4.
+unfold m'.
+destruct s' as [|p|p].
+repeat split ; try easy.
+fold d.
+omega.
+fold (Zpower beta (Zpos p)).
+split.
+replace (Zpos p) with (Zpos p + e' - e')%Z by ring.
+rewrite <- F2R_change_exp.
+apply (f_equal (fun v => F2R (Float beta m v))).
+ring.
+assert (0 < Zpos p)%Z by easy.
+omega.
+split.
+rewrite Zdigits_mult_Zpower.
+fold d.
+omega.
+apply sym_not_eq.
+now apply Zlt_not_eq.
+easy.
+apply Zmult_lt_0_compat.
+exact Hm.
+now apply Zpower_gt_0.
+now elim He2.
+clearbody m'.
+destruct Hs as (Hs1, (Hs2, Hs3)).
+generalize (Zsqrt_correct m' (Zlt_le_weak _ _ Hs3)).
+destruct (Zsqrt m') as (q, r).
+intros (Hq, Hr).
+rewrite <- Hs1. clear Hs1.
+split.
+(* . mantissa width *)
+apply Zmult_le_reg_r with 2%Z.
+easy.
+rewrite Zmult_comm.
+apply Zle_trans with (1 := Hs2).
+rewrite Hq.
+apply Zle_trans with (Zdigits beta (q + q + q * q)).
+apply Zdigits_le.
+rewrite <- Hq.
+now apply Zlt_le_weak.
+omega.
+replace (Zdigits beta q * 2)%Z with (Zdigits beta q + Zdigits beta q)%Z by ring.
+apply Zdigits_mult_strong.
+omega.
+omega.
+(* . round *)
+unfold inbetween_float, F2R. simpl.
+rewrite sqrt_mult.
+2: now apply (Z2R_le 0) ; apply Zlt_le_weak.
+2: apply Rlt_le ; apply bpow_gt_0.
+destruct (Zeven_ex e') as (e2, Hev).
+rewrite He1, Zplus_0_r in Hev. clear He1.
+rewrite Hev.
+replace (Zdiv2 (2 * e2)) with e2 by now case e2.
+replace (2 * e2)%Z with (e2 + e2)%Z by ring.
+rewrite bpow_plus.
+fold (Rsqr (bpow e2)).
+rewrite sqrt_Rsqr.
+2: apply Rlt_le ; apply bpow_gt_0.
+apply inbetween_mult_compat.
+apply bpow_gt_0.
+rewrite Hq.
+case Zeq_bool_spec ; intros Hr'.
+(* .. r = 0 *)
+rewrite Hr', Zplus_0_r, Z2R_mult.
+fold (Rsqr (Z2R q)).
+rewrite sqrt_Rsqr.
+now constructor.
+apply (Z2R_le 0).
+omega.
+(* .. r <> 0 *)
+constructor.
+split.
+(* ... bounds *)
+apply Rle_lt_trans with (sqrt (Z2R (q * q))).
+rewrite Z2R_mult.
+fold (Rsqr (Z2R q)).
+rewrite sqrt_Rsqr.
+apply Rle_refl.
+apply (Z2R_le 0).
+omega.
+apply sqrt_lt_1.
+rewrite Z2R_mult.
+apply Rle_0_sqr.
+rewrite <- Hq.
+apply (Z2R_le 0).
+now apply Zlt_le_weak.
+apply Z2R_lt.
+omega.
+apply Rlt_le_trans with (sqrt (Z2R ((q + 1) * (q + 1)))).
+apply sqrt_lt_1.
+rewrite <- Hq.
+apply (Z2R_le 0).
+now apply Zlt_le_weak.
+rewrite Z2R_mult.
+apply Rle_0_sqr.
+apply Z2R_lt.
+ring_simplify.
+omega.
+rewrite Z2R_mult.
+fold (Rsqr (Z2R (q + 1))).
+rewrite sqrt_Rsqr.
+apply Rle_refl.
+apply (Z2R_le 0).
+omega.
+(* ... location *)
+rewrite Rcompare_half_r.
+rewrite <- Rcompare_sqr.
+replace ((2 * sqrt (Z2R (q * q + r))) * (2 * sqrt (Z2R (q * q + r))))%R
+  with (4 * Rsqr (sqrt (Z2R (q * q + r))))%R by (unfold Rsqr ; ring).
+rewrite Rsqr_sqrt.
+change 4%R with (Z2R 4).
+rewrite <- Z2R_plus, <- 2!Z2R_mult.
+rewrite Rcompare_Z2R.
+replace ((q + (q + 1)) * (q + (q + 1)))%Z with (4 * (q * q) + 4 * q + 1)%Z by ring.
+generalize (Zle_cases r q).
+case (Zle_bool r q) ; intros Hr''.
+change (4 * (q * q + r) < 4 * (q * q) + 4 * q + 1)%Z.
+omega.
+change (4 * (q * q + r) > 4 * (q * q) + 4 * q + 1)%Z.
+omega.
+rewrite <- Hq.
+apply (Z2R_le 0).
+now apply Zlt_le_weak.
+apply Rmult_le_pos.
+now apply (Z2R_le 0 2).
+apply sqrt_ge_0.
+rewrite <- Z2R_plus.
+apply (Z2R_le 0).
+omega.
+Qed.
+
+End Fcalc_sqrt.
-- 
cgit