Add sqrt in Numbers

 As for power recently, we add a specification in NZ,N,Z,
 derived properties, implementations for nat, N, Z, BigN, BigZ.

  - For nat, this sqrt is brand new :-), cf NPeano.v

  - For Z, we rework what was in Zsqrt: same algorithm,
    no more refine but a pure function, based now on a sqrt
    for positive, from which we derive a Nsqrt and a Zsqrt.
    For the moment, the old Zsqrt.v file is kept as Zsqrt_compat.v.
    It is not loaded by default by Require ZArith.
    New definitions are now in Psqrt.v, Zsqrt_def.v and Nsqrt_def.v

  - For BigN, BigZ, we changed the specifications to refer to Zsqrt
    instead of using characteristic inequations.

 On the way, many extensions, in particular BinPos (lemmas about order),
 NZMulOrder (results about squares)

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/coq/trunk@13564 85f007b7-540e-0410-9357-904b9bb8a0f7

32 files changed, 919 insertions, 107 deletions

diff --git a/theories/Arith/Mult.v b/theories/Arith/Mult.v
index fed1bbca5..479138a98 100644
--- a/theories/Arith/Mult.v
+++ b/theories/Arith/Mult.v
@@ -175,19 +175,22 @@ Qed.
 Lemma mult_S_lt_compat_l : forall n m p, m < p -> S n * m < S n * p.
 Proof.
   induction n; intros; simpl in *.
-    rewrite <- 2! plus_n_O; assumption.
+    rewrite <- 2 plus_n_O; assumption.
     auto using plus_lt_compat.
 Qed.
 
 Hint Resolve mult_S_lt_compat_l: arith.
 
+Lemma mult_lt_compat_l : forall n m p, n < m -> 0 < p -> p * n < p * m.
+Proof.
+ intros m n p H Hp. destruct p. elim (lt_irrefl _ Hp).
+ now apply mult_S_lt_compat_l.
+Qed.
+
 Lemma mult_lt_compat_r : forall n m p, n < m -> 0 < p -> n * p < m * p.
 Proof.
-  intros m n p H H0.
-  induction p.
-  elim (lt_irrefl _ H0).
-  rewrite mult_comm.
-  replace (n * S p) with (S p * n); auto with arith.
+  intros m n p H Hp. destruct p. elim (lt_irrefl _ Hp).
+  rewrite (mult_comm m), (mult_comm n). now apply mult_S_lt_compat_l.
 Qed.
 
 Lemma mult_S_le_reg_l : forall n m p, S n * m <= S n * p -> m <= p.
diff --git a/theories/NArith/BinNat.v b/theories/NArith/BinNat.v
index a1a12017a..21aaabcbc 100644
--- a/theories/NArith/BinNat.v
+++ b/theories/NArith/BinNat.v
@@ -135,6 +135,11 @@ Infix "<" := Nlt : N_scope.
 Infix ">=" := Nge : N_scope.
 Infix ">" := Ngt : N_scope.
 
+Notation "x <= y <= z" := (x <= y /\ y <= z) : N_scope.
+Notation "x <= y < z" := (x <= y /\ y < z) : N_scope.
+Notation "x < y < z" := (x < y /\ y < z) : N_scope.
+Notation "x < y <= z" := (x < y /\ y <= z) : N_scope.
+
 (** Min and max *)
 
 Definition Nmin (n n' : N) := match Ncompare n n' with
@@ -463,6 +468,24 @@ apply Ncompare_Eq_eq; auto.
 apply Ngt_Nlt; auto.
 Qed.
 
+(** Order and operations *)
+
+(** NB : Many more results will come from NBinary, the Number instantiation.
+    The next lemma is useful for proving the correctness of the division.
+*)
+
+Lemma Nplus_lt_cancel_l : forall n m p, p+n < p+m -> n<m.
+Proof.
+ intros. destruct p. simpl; auto.
+ destruct n; destruct m.
+ elim (Nlt_irrefl _ H).
+ red; auto.
+ rewrite Nplus_0_r in H. simpl in H.
+ red in H. simpl in H.
+ elim (Plt_not_plus_l _ _ H).
+ now apply (Pplus_lt_mono_l p).
+Qed.
+
 (** 0 is the least natural number *)
 
 Theorem Ncompare_0 : forall n : N, Ncompare n N0 <> Lt.
diff --git a/theories/NArith/BinPos.v b/theories/NArith/BinPos.v
index 14489ebda..d334d42e9 100644
--- a/theories/NArith/BinPos.v
+++ b/theories/NArith/BinPos.v
@@ -746,6 +746,19 @@ Proof.
   intros [p|p| ] [q|q| ] H; destr_eq H; auto.
 Qed.
 
+(** Square *)
+
+Lemma Psquare_xO : forall p, p~0 * p~0 = (p*p)~0~0.
+Proof.
+ intros. simpl. now rewrite Pmult_comm.
+Qed.
+
+Lemma Psquare_xI : forall p, p~1 * p~1 = (p*p+p)~0~1.
+Proof.
+ intros. simpl. rewrite Pmult_comm. simpl. f_equal.
+ rewrite Pplus_assoc, Pplus_diag. simpl. now rewrite Pplus_comm.
+Qed.
+
 (*********************************************************************)
 (** Properties of boolean equality *)
 
@@ -959,11 +972,29 @@ Proof.
  exact Pcompare_1.
 Qed.
 
+Lemma Plt_1_succ : forall n, 1 < Psucc n.
+Proof.
+ intros. apply Pcompare_p_Sq. destruct n; auto.
+Qed.
+
 Lemma Plt_lt_succ : forall n m : positive, n < m -> n < Psucc m.
 Proof.
   unfold Plt; intros n m H; apply <- Pcompare_p_Sq; auto.
 Qed.
 
+Lemma Psucc_lt_compat : forall n m, n < m <-> Psucc n < Psucc m.
+Proof.
+ unfold Plt. induction n; destruct m; simpl; auto; split; try easy; intros.
+ now apply Pcompare_Lt_eq_Lt, Pcompare_p_Sq, IHn, Pcompare_Gt_Lt.
+ now apply Pcompare_eq_Lt, IHn, Pcompare_p_Sq, Pcompare_Lt_eq_Lt.
+ destruct (Psucc n); discriminate.
+ now apply Pcompare_eq_Lt, Pcompare_p_Sq, Pcompare_Lt_eq_Lt.
+ now apply Pcompare_Lt_eq_Lt, Pcompare_p_Sq, Pcompare_Gt_Lt.
+ destruct n; discriminate.
+ apply Plt_1_succ.
+ destruct m; auto.
+Qed.
+
 Lemma Plt_irrefl : forall p : positive, ~ p < p.
 Proof.
   unfold Plt; intro p; rewrite Pcompare_refl; discriminate.
@@ -994,6 +1025,51 @@ Proof.
   destruct ((p ?= q) Eq); intuition; discriminate.
 Qed.
 
+(** Strict order and operations *)
+
+Lemma Pplus_lt_mono_l : forall p q r, q<r <-> p+q < p+r.
+Proof.
+ induction p using Prect.
+ intros q r. rewrite <- 2 Pplus_one_succ_l. apply Psucc_lt_compat.
+ intros q r. rewrite 2 Pplus_succ_permute_l.
+ eapply iff_trans; [ eapply IHp | eapply Psucc_lt_compat ].
+Qed.
+
+Lemma Pplus_lt_mono : forall p q r s, p<q -> r<s -> p+r<q+s.
+Proof.
+ intros. apply Plt_trans with (p+s).
+ now apply Pplus_lt_mono_l.
+ rewrite (Pplus_comm p), (Pplus_comm q). now apply Pplus_lt_mono_l.
+Qed.
+
+Lemma Pmult_lt_mono_l : forall p q r, q<r -> p*q < p*r.
+Proof.
+ induction p using Prect; auto.
+ intros q r H. rewrite 2 Pmult_Sn_m.
+ apply Pplus_lt_mono; auto.
+Qed.
+
+Lemma Pmult_lt_mono : forall p q r s, p<q -> r<s -> p*r < q*s.
+Proof.
+ intros. apply Plt_trans with (p*s).
+ now apply Pmult_lt_mono_l.
+ rewrite (Pmult_comm p), (Pmult_comm q). now apply Pmult_lt_mono_l.
+Qed.
+
+Lemma Plt_plus_r : forall p q, p < p+q.
+Proof.
+ induction q using Prect.
+ rewrite <- Pplus_one_succ_r. apply Pcompare_p_Sp.
+ apply Plt_trans with (p+q); auto.
+ apply Pplus_lt_mono_l, Pcompare_p_Sp.
+Qed.
+
+Lemma Plt_not_plus_l : forall p q, ~ p+q < p.
+Proof.
+ intros p q H. elim (Plt_irrefl p).
+ apply Plt_trans with (p+q); auto using Plt_plus_r.
+Qed.
+
 
 (**********************************************************************)
 (** Properties of subtraction on binary positive numbers *)
diff --git a/theories/NArith/NArith.v b/theories/NArith/NArith.v
index 1ba3f2992..5cfd0c416 100644
--- a/theories/NArith/NArith.v
+++ b/theories/NArith/NArith.v
@@ -12,6 +12,7 @@ Require Export BinPos.
 Require Export BinNat.
 Require Export Nnat.
 Require Export Ndiv_def.
+Require Export Nsqrt_def.
 Require Export Ndigits.
 Require Export NArithRing.
 Require NBinary.
diff --git a/theories/NArith/Ndiv_def.v b/theories/NArith/Ndiv_def.v
index 03d6aceaa..2a3fd152a 100644
--- a/theories/NArith/Ndiv_def.v
+++ b/theories/NArith/Ndiv_def.v
@@ -7,7 +7,7 @@
 (************************************************************************)
 
 
-Require Import BinPos BinNat Pnat Plus Compare_dec.
+Require Import BinPos BinNat.
 
 Local Open Scope N_scope.
 
@@ -76,45 +76,11 @@ Proof.
   simpl; f_equal; apply Pplus_comm.
 Qed.
 
-Lemma Plt_add_cancel_l : forall n m p, (p+n < p+m -> n<m)%positive.
-Proof.
- intros n m p H.
- unfold Plt in *.
- rewrite nat_of_P_compare_morphism in *.
- rewrite 2 nat_of_P_plus_morphism in *.
- apply nat_compare_lt.
- apply nat_compare_lt in H.
- eapply plus_lt_reg_l; eauto.
-Qed.
-
-Lemma Plt_add_not_lt : forall n m, ~(n+m < n)%positive.
-Proof.
- intros n m H.
- unfold Plt in *.
- rewrite nat_of_P_compare_morphism in *.
- rewrite nat_of_P_plus_morphism in *.
- apply nat_compare_lt in H.
- absurd (nat_of_P m < 0)%nat.
- red; inversion 1.
- apply plus_lt_reg_l with (nat_of_P n). now rewrite plus_0_r.
-Qed.
-
-Lemma Nlt_add_cancel_l : forall n m p, p+n < p+m -> n<m.
-Proof.
- intros. destruct p. simpl; auto.
- destruct n; destruct m.
- elim (Nlt_irrefl _ H).
- red; auto.
- rewrite Nplus_0_r in H. simpl in H.
- elim (Plt_add_not_lt _ _ H).
- apply Plt_add_cancel_l with p; auto.
-Qed.
-
 Lemma NPgeb_ineq0 : forall a p, a < Npos p -> NPgeb (2*a) p = true ->
  2*a - Npos p < Npos p.
 Proof.
 intros a p LT GE.
-apply Nlt_add_cancel_l with (Npos p).
+apply Nplus_lt_cancel_l with (Npos p).
 rewrite Nplus_comm.
 generalize (NPgeb_correct (2*a) p). rewrite GE. intros <-.
 rewrite <- (Nmult_1_l (Npos p)). rewrite <- Nmult_plus_distr_r.
@@ -125,7 +91,7 @@ Lemma NPgeb_ineq1 : forall a p, a < Npos p -> NPgeb (2*a+1) p = true ->
   (2*a+1) - Npos p < Npos p.
 Proof.
 intros a p LT GE.
-apply Nlt_add_cancel_l with (Npos p).
+apply Nplus_lt_cancel_l with (Npos p).
 rewrite Nplus_comm.
 generalize (NPgeb_correct (2*a+1) p). rewrite GE. intros <-.
 rewrite <- (Nmult_1_l (Npos p)). rewrite <- Nmult_plus_distr_r.
diff --git a/theories/NArith/Nsqrt_def.v b/theories/NArith/Nsqrt_def.v
new file mode 100644
index 000000000..3e25d8316
--- /dev/null
+++ b/theories/NArith/Nsqrt_def.v
@@ -0,0 +1,50 @@
+(************************************************************************)
+(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
+(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
+(*   \VV/  **************************************************************)
+(*    //   *      This file is distributed under the terms of the       *)
+(*         *       GNU Lesser General Public License Version 2.1        *)
+(************************************************************************)
+
+(** Definition of a square root function for N. *)
+
+Require Import BinPos BinNat Psqrt.
+
+Open Scope N_scope.
+
+Definition Nsqrtrem n :=
+ match n with
+  | N0 => (N0, N0)
+  | Npos p =>
+    match Psqrtrem p with
+     | (s, IsPos r) => (Npos s, Npos r)
+     | (s, _) => (Npos s, N0)
+    end
+ end.
+
+Definition Nsqrt n :=
+ match n with
+  | N0 => N0
+  | Npos p => Npos (Psqrt p)
+ end.
+
+Lemma Nsqrtrem_spec : forall n,
+ let (s,r) := Nsqrtrem n in n = s*s + r /\ r <= 2*s.
+Proof.
+ destruct n. now split.
+ generalize (Psqrtrem_spec p). simpl.
+ destruct 1; simpl; subst; now split.
+Qed.
+
+Lemma Nsqrt_spec : forall n,
+ let s := Nsqrt n in s*s <= n < (Nsucc s)*(Nsucc s).
+Proof.
+ destruct n. now split. apply (Psqrt_spec p).
+Qed.
+
+Lemma Nsqrtrem_sqrt : forall n, fst (Nsqrtrem n) = Nsqrt n.
+Proof.
+ destruct n. reflexivity.
+ unfold Nsqrtrem, Nsqrt, Psqrt.
+ destruct (Psqrtrem p) as (s,r). now destruct r.
+Qed.
\ No newline at end of file
diff --git a/theories/NArith/Psqrt.v b/theories/NArith/Psqrt.v
new file mode 100644
index 000000000..1d85f14a2
--- /dev/null
+++ b/theories/NArith/Psqrt.v
@@ -0,0 +1,123 @@
+(************************************************************************)
+(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
+(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
+(*   \VV/  **************************************************************)
+(*    //   *      This file is distributed under the terms of the       *)
+(*         *       GNU Lesser General Public License Version 2.1        *)
+(************************************************************************)
+
+Require Import BinPos.
+
+Open Scope positive_scope.
+
+Definition Pleb x y :=
+ match Pcompare x y Eq with Gt => false | _ => true end.
+
+(** A Square Root function for positive numbers *)
+
+(** We procede by blocks of two digits : if p is written qbb'
+    then sqrt(p) will be sqrt(q)~0 or sqrt(q)~1.
+    For deciding easily in which case we are, we store the remainder
+    (as a positive_mask, since it can be null).
+    Instead of copy-pasting the following code four times, we
+    factorize as an auxiliary function, with f and g being either
+    xO or xI depending of the initial digits.
+    NB: (Pminus_mask (g (f 1)) 4) is a hack, morally it's g (f 0).
+*)
+
+Definition Psqrtrem_step (f g:positive->positive) p :=
+ match p with
+  | (s, IsPos r) =>
+    let s' := s~0~1 in
+    let r' := g (f r) in
+    if Pleb s' r' then (s~1, Pminus_mask r' s')
+    else (s~0, IsPos r')
+  | (s,_)  => (s~0, Pminus_mask (g (f 1)) 4)
+ end.
+
+Fixpoint Psqrtrem p : positive * positive_mask :=
+ match p with
+  | 1 => (1,IsNul)
+  | 2 => (1,IsPos 1)
+  | 3 => (1,IsPos 2)
+  | p~0~0 => Psqrtrem_step xO xO (Psqrtrem p)
+  | p~0~1 => Psqrtrem_step xO xI (Psqrtrem p)
+  | p~1~0 => Psqrtrem_step xI xO (Psqrtrem p)
+  | p~1~1 => Psqrtrem_step xI xI (Psqrtrem p)
+ end.
+
+Definition Psqrt p := fst (Psqrtrem p).
+
+(** An inductive relation for specifying sqrt results *)
+
+Inductive PSQRT : positive*positive_mask -> positive -> Prop :=
+ | PSQRT_exact : forall s x, x=s*s -> PSQRT (s,IsNul) x
+ | PSQRT_approx : forall s r x, x=s*s+r -> r <= s~0 -> PSQRT (s,IsPos r) x.
+
+(** Correctness proofs *)
+
+Lemma Psqrtrem_step_spec : forall f g p x,
+ (f=xO \/ f=xI) -> (g=xO \/ g=xI) ->
+ PSQRT p x -> PSQRT (Psqrtrem_step f g p) (g (f x)).
+Proof.
+intros f g _ _ Hf Hg [ s _ -> | s r _ -> Hr ].
+(* exact *)
+unfold Psqrtrem_step.
+destruct Hf,Hg; subst; simpl Pminus_mask;
+ constructor; try discriminate; now rewrite Psquare_xO.
+(* approx *)
+assert (Hfg : forall p q, g (f (p+q)) = p~0~0 + g (f q))
+ by (intros; destruct Hf, Hg; now subst).
+unfold Psqrtrem_step. unfold Pleb.
+case Pcompare_spec; [intros EQ | intros LT | intros GT].
+(* - EQ *)
+rewrite <- EQ. rewrite Pminus_mask_diag.
+destruct Hg; subst g; try discriminate.
+destruct Hf; subst f; try discriminate.
+injection EQ; clear EQ; intros <-.
+constructor. now rewrite Psquare_xI.
+(* - LT *)
+destruct (Pminus_mask_Gt (g (f r)) (s~0~1)) as (y & -> & H & _).
+change Eq with (CompOpp Eq). now rewrite <- Pcompare_antisym, LT.
+constructor.
+rewrite Hfg, <- H. now rewrite Psquare_xI, Pplus_assoc.
+apply Ple_lteq, Pcompare_p_Sq in Hr; change (r < s~1) in Hr.
+apply Ple_lteq, Pcompare_p_Sq; change (y < s~1~1).
+apply Pplus_lt_mono_l with (s~0~1).
+rewrite H. simpl. rewrite Pplus_carry_spec, Pplus_diag. simpl.
+set (s1:=s~1) in *; clearbody s1.
+destruct Hf,Hg; subst; red; simpl;
+  apply Hr || apply Pcompare_eq_Lt, Hr.
+(* - GT *)
+constructor.
+rewrite Hfg. now rewrite Psquare_xO.
+apply Ple_lteq, Pcompare_p_Sq, GT.
+Qed.
+
+Lemma Psqrtrem_spec : forall p, PSQRT (Psqrtrem p) p.
+Proof.
+fix 1.
+ destruct p; try destruct p; try (constructor; easy);
+  apply Psqrtrem_step_spec; auto.
+Qed.
+
+Lemma Psqrt_spec : forall p,
+ let s := Psqrt p in s*s <= p < (Psucc s)*(Psucc s).
+Proof.
+ intros p. simpl.
+ assert (H:=Psqrtrem_spec p).
+ unfold Psqrt in *. destruct Psqrtrem as (s,rm); simpl.
+ inversion_clear H; subst.
+ (* exact *)
+ split. red. rewrite Pcompare_refl. discriminate.
+  apply Pmult_lt_mono; apply Pcompare_p_Sp.
+ (* approx *)
+ split.
+ apply Ple_lteq; left. apply Plt_plus_r.
+ rewrite (Pplus_one_succ_l).
+ rewrite Pmult_plus_distr_r, !Pmult_plus_distr_l.
+ rewrite !Pmult_1_r. simpl (1*s).
+ rewrite <- Pplus_assoc, (Pplus_assoc s s), Pplus_diag, Pplus_assoc.
+ rewrite (Pplus_comm (_+_)). apply Pplus_lt_mono_l.
+ rewrite <- Pplus_one_succ_l. now apply Pcompare_p_Sq, Ple_lteq.
+Qed.
diff --git a/theories/NArith/vo.itarget b/theories/NArith/vo.itarget
index 7692a6368..0caf0b249 100644
--- a/theories/NArith/vo.itarget
+++ b/theories/NArith/vo.itarget
@@ -9,3 +9,5 @@ Ndiv_def.vo
 Pnat.vo
 POrderedType.vo
 Pminmax.vo
+Psqrt.vo
+Nsqrt_def.vo
diff --git a/theories/Numbers/Cyclic/ZModulo/ZModulo.v b/theories/Numbers/Cyclic/ZModulo/ZModulo.v
index 6fcb4cf91..2d9eb395a 100644
--- a/theories/Numbers/Cyclic/ZModulo/ZModulo.v
+++ b/theories/Numbers/Cyclic/ZModulo/ZModulo.v
@@ -574,27 +574,19 @@ Section ZModulo.
  generalize (Z_mod_lt [|x|] 2); omega.
  Qed.
 
- Definition sqrt x := Zsqrt_plain [|x|].
+ Definition sqrt x := Zsqrt [|x|].
  Lemma spec_sqrt : forall x,
        [|sqrt x|] ^ 2 <= [|x|] < ([|sqrt x|] + 1) ^ 2.
  Proof.
  intros.
  unfold sqrt.
  repeat rewrite Zpower_2.
- replace [|Zsqrt_plain [|x|]|] with (Zsqrt_plain [|x|]).
- apply Zsqrt_interval; auto with zarith.
+ replace [|Zsqrt [|x|]|] with (Zsqrt [|x|]).
+ apply Zsqrt_spec; auto with zarith.
  symmetry; apply Zmod_small.
- split.
- apply Zsqrt_plain_is_pos; auto with zarith.
-
- cut (Zsqrt_plain [|x|] <= (wB-1)); try omega.
- rewrite <- (Zsqrt_square_id (wB-1)).
- apply Zsqrt_le.
- split; auto.
- apply Zle_trans with (wB-1); auto with zarith.
- generalize (spec_to_Z x); auto with zarith.
- apply Zsquare_le.
- generalize wB_pos; auto with zarith.
+ split. apply Z.sqrt_nonneg; auto.
+ apply Zle_lt_trans with [|x|]; auto.
+ apply Z.sqrt_le_lin; auto.
  Qed.
 
  Definition sqrt2 x y :=
@@ -602,7 +594,7 @@ Section ZModulo.
   match z with
    | Z0 => (0, C0 0)
    | Zpos p =>
-      let (s,r,_,_) := sqrtrempos p in
+      let (s,r) := Zsqrtrem (Zpos p) in
       (s, if Z_lt_le_dec r wB then C0 r else C1 (r-wB))
    | Zneg _ => (0, C0 0)
   end.
@@ -618,11 +610,12 @@ Section ZModulo.
  remember ([|x|]*wB+[|y|]) as z.
  destruct z.
  auto with zarith.
- destruct sqrtrempos; intros.
+ generalize (Zsqrtrem_spec (Zpos p)).
+ destruct Zsqrtrem as (s,r); intros [U V]; auto with zarith.
  assert (s < wB).
   destruct (Z_lt_le_dec s wB); auto.
   assert (wB * wB <= Zpos p).
-   rewrite e.
+   rewrite U.
    apply Zle_trans with (s*s); try omega.
    apply Zmult_le_compat; generalize wB_pos; auto with zarith.
   assert (Zpos p < wB*wB).
diff --git a/theories/Numbers/Integer/Abstract/ZAxioms.v b/theories/Numbers/Integer/Abstract/ZAxioms.v
index 38855a85d..4f88008be 100644
--- a/theories/Numbers/Integer/Abstract/ZAxioms.v
+++ b/theories/Numbers/Integer/Abstract/ZAxioms.v
@@ -9,7 +9,7 @@
 (************************************************************************)
 
 Require Export NZAxioms.
-Require Import NZPow.
+Require Import NZPow NZSqrt.
 
 (** We obtain integers by postulating that successor of predecessor
     is identity. *)
@@ -76,15 +76,20 @@ Module Type ZPowSpecNeg (Import Z : ZAxiomsMiniSig')(Import P : Pow' Z).
  Axiom pow_neg : forall a b, b<0 -> a^b == 0.
 End ZPowSpecNeg.
 
+(** Same for the sqrt function. *)
+
+Module Type ZSqrtSpecNeg (Import Z : ZAxiomsMiniSig')(Import P : Sqrt' Z).
+ Axiom sqrt_neg : forall a, a<0 -> √a == 0.
+End ZSqrtSpecNeg.
 
 (** Let's group everything *)
 
 Module Type ZAxiomsSig :=
   ZAxiomsMiniSig <+ HasCompare <+ HasAbs <+ HasSgn <+ Parity
-   <+ NZPow.NZPow <+ ZPowSpecNeg.
+   <+ NZPow.NZPow <+ ZPowSpecNeg <+ NZSqrt.NZSqrt <+ ZSqrtSpecNeg.
 Module Type ZAxiomsSig' :=
   ZAxiomsMiniSig' <+ HasCompare <+ HasAbs <+ HasSgn <+ Parity
-   <+ NZPow.NZPow' <+ ZPowSpecNeg.
+   <+ NZPow.NZPow' <+ ZPowSpecNeg <+ NZSqrt.NZSqrt' <+ ZSqrtSpecNeg.
 
 
 (** Division is left apart, since many different flavours are available *)
diff --git a/theories/Numbers/Integer/Abstract/ZMulOrder.v b/theories/Numbers/Integer/Abstract/ZMulOrder.v
index 25989b2d4..1010a0f2f 100644
--- a/theories/Numbers/Integer/Abstract/ZMulOrder.v
+++ b/theories/Numbers/Integer/Abstract/ZMulOrder.v
@@ -92,14 +92,7 @@ Qed.
 
 Notation mul_nonpos := le_mul_0 (only parsing).
 
-Theorem le_0_square : forall n, 0 <= n * n.
-Proof.
-intro n; destruct (neg_nonneg_cases n).
-apply lt_le_incl; now apply mul_neg_neg.
-now apply mul_nonneg_nonneg.
-Qed.
-
-Notation square_nonneg := le_0_square (only parsing).
+Notation le_0_square := square_nonneg (only parsing).
 
 Theorem nlt_square_0 : forall n, ~ n * n < 0.
 Proof.
diff --git a/theories/Numbers/Integer/Abstract/ZProperties.v b/theories/Numbers/Integer/Abstract/ZProperties.v
index 7b9c6f452..8b34e5b2d 100644
--- a/theories/Numbers/Integer/Abstract/ZProperties.v
+++ b/theories/Numbers/Integer/Abstract/ZProperties.v
@@ -11,6 +11,5 @@ Require Export ZAxioms ZMaxMin ZSgnAbs ZParity ZPow.
 (** This functor summarizes all known facts about Z. *)
 
 Module Type ZProp (Z:ZAxiomsSig) :=
- ZMaxMinProp Z <+ ZSgnAbsProp Z <+ ZParityProp Z <+ ZPowProp Z.
-
-
+ ZMaxMinProp Z <+ ZSgnAbsProp Z <+ ZParityProp Z <+ ZPowProp Z
+ <+ NZSqrt.NZSqrtProp Z Z.
diff --git a/theories/Numbers/Integer/BigZ/ZMake.v b/theories/Numbers/Integer/BigZ/ZMake.v
index f9490cc9c..099554cd0 100644
--- a/theories/Numbers/Integer/BigZ/ZMake.v
+++ b/theories/Numbers/Integer/BigZ/ZMake.v
@@ -338,15 +338,14 @@ Module Make (N:NType) <: ZType.
   | Neg nx => Neg N.zero
   end.
 
- Theorem spec_sqrt: forall x, 0 <= to_Z x ->
-   to_Z (sqrt x) ^ 2 <= to_Z x < (to_Z (sqrt x) + 1) ^ 2.
+ Theorem spec_sqrt: forall x, to_Z (sqrt x) = Zsqrt (to_Z x).
  Proof.
- unfold to_Z, sqrt; intros [x | x] H.
-   exact (N.spec_sqrt x).
- replace (N.to_Z x) with 0.
-   rewrite N.spec_0; simpl Zpower; unfold Zpower_pos, iter_pos;
-    auto with zarith.
- generalize (N.spec_pos x); auto with zarith.
+  destruct x as [p|p]; simpl.
+  apply N.spec_sqrt.
+  rewrite N.spec_0.
+  destruct (Z_le_lt_eq_dec _ _ (N.spec_pos p)) as [LT|EQ].
+  rewrite Zsqrt_neg; auto with zarith.
+  now rewrite <- EQ.
  Qed.
 
  Definition div_eucl x y :=
diff --git a/theories/Numbers/Integer/Binary/ZBinary.v b/theories/Numbers/Integer/Binary/ZBinary.v
index 5f8728394..ee75e4aa1 100644
--- a/theories/Numbers/Integer/Binary/ZBinary.v
+++ b/theories/Numbers/Integer/Binary/ZBinary.v
@@ -10,7 +10,7 @@
 
 
 Require Import ZAxioms ZProperties.
-Require Import BinInt Zcompare Zorder ZArith_dec Zbool Zeven.
+Require Import BinInt Zcompare Zorder ZArith_dec Zbool Zeven Zsqrt_def.
 
 Local Open Scope Z_scope.
 
@@ -174,6 +174,17 @@ Definition pow_succ_r := Zpow_succ_r.
 Definition pow_neg := Zpow_neg.
 Definition pow := Zpow.
 
+(** Sqrt *)
+
+(** NB : we use a new Zsqrt defined in Zsqrt_def, the previous
+   module Zsqrt.v is now Zsqrt_compat.v *)
+
+Program Instance sqrt_wd : Proper (eq==>eq) Zsqrt.
+
+Definition sqrt_spec := Zsqrt_spec.
+Definition sqrt_neg := Zsqrt_neg.
+Definition sqrt := Zsqrt.
+
 (** We define [eq] only here to avoid refering to this [eq] above. *)
 
 Definition eq := (@eq Z).
diff --git a/theories/Numbers/Integer/SpecViaZ/ZSig.v b/theories/Numbers/Integer/SpecViaZ/ZSig.v
index be201f2d6..37f5b294e 100644
--- a/theories/Numbers/Integer/SpecViaZ/ZSig.v
+++ b/theories/Numbers/Integer/SpecViaZ/ZSig.v
@@ -78,13 +78,12 @@ Module Type ZType.
  Parameter spec_pow_pos: forall x n, [pow_pos x n] = [x] ^ Zpos n.
  Parameter spec_pow_N: forall x n, [pow_N x n] = [x] ^ Z_of_N n.
  Parameter spec_pow: forall x n, [pow x n] = [x] ^ [n].
- Parameter spec_sqrt: forall x, 0 <= [x] ->
-   [sqrt x] ^ 2 <= [x] < ([sqrt x] + 1) ^ 2.
+ Parameter spec_sqrt: forall x, [sqrt x] = Zsqrt [x].
  Parameter spec_div_eucl: forall x y,
    let (q,r) := div_eucl x y in ([q], [r]) = Zdiv_eucl [x] [y].
  Parameter spec_div: forall x y, [div x y] = [x] / [y].
  Parameter spec_modulo: forall x y, [modulo x y] = [x] mod [y].
- Parameter spec_gcd: forall a b, [gcd a b] = Zgcd (to_Z a) (to_Z b).
+ Parameter spec_gcd: forall a b, [gcd a b] = Zgcd [a] [b].
  Parameter spec_sgn : forall x, [sgn x] = Zsgn [x].
  Parameter spec_abs : forall x, [abs x] = Zabs [x].
  Parameter spec_even : forall x, even x = Zeven_bool [x].
diff --git a/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v b/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
index 3e6375543..d632d2260 100644
--- a/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
+++ b/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
@@ -18,7 +18,7 @@ Module ZTypeIsZAxioms (Import Z : ZType').
 
 Hint Rewrite
  spec_0 spec_1 spec_2 spec_add spec_sub spec_pred spec_succ
- spec_mul spec_opp spec_of_Z spec_div spec_modulo
+ spec_mul spec_opp spec_of_Z spec_div spec_modulo spec_sqrt
  spec_compare spec_eq_bool spec_max spec_min spec_abs spec_sgn
  spec_pow spec_even spec_odd
  : zsimpl.
@@ -278,6 +278,21 @@ Proof.
  intros a b. red. now rewrite spec_pow_N, spec_pow_pos.
 Qed.
 
+(** Sqrt *)
+
+Program Instance sqrt_wd : Proper (eq==>eq) sqrt.
+
+Lemma sqrt_spec : forall n, 0<=n ->
+ (sqrt n)*(sqrt n) <= n /\ n < (succ (sqrt n))*(succ (sqrt n)).
+Proof.
+ intros n. zify. apply Zsqrt_spec.
+Qed.
+
+Lemma sqrt_neg : forall n, n<0 -> sqrt n == 0.
+Proof.
+ intros n. zify. apply Zsqrt_neg.
+Qed.
+
 (** Even / Odd *)
 
 Definition Even n := exists m, n == 2*m.
diff --git a/theories/Numbers/NatInt/NZMulOrder.v b/theories/Numbers/NatInt/NZMulOrder.v
index 530044ab0..8f1b8cb6e 100644
--- a/theories/Numbers/NatInt/NZMulOrder.v
+++ b/theories/Numbers/NatInt/NZMulOrder.v
@@ -295,11 +295,60 @@ assumption. assert (F : m * m < n * n) by now apply square_lt_mono_nonneg.
 apply -> lt_nge in F. false_hyp H2 F.
 Qed.
 
-Theorem mul_2_mono_l : forall n m, n < m -> 1 + (1 + 1) * n < (1 + 1) * m.
+Theorem mul_2_mono_l : forall n m, n < m -> 1 + 2 * n < 2 * m.
 Proof.
-intros n m. rewrite <- le_succ_l, (mul_le_mono_pos_l (S n) m (1 + 1)).
-rewrite !mul_add_distr_r; nzsimpl; now rewrite le_succ_l.
-apply add_pos_pos; now apply lt_0_1.
+intros n m. rewrite <- le_succ_l, (mul_le_mono_pos_l (S n) m two).
+rewrite two_succ. nzsimpl. now rewrite le_succ_l.
+order'.
+Qed.
+
+(** A few results about squares *)
+
+Lemma square_nonneg : forall a, 0 <= a * a.
+Proof.
+ intros. rewrite <- (mul_0_r a). destruct (le_gt_cases a 0).
+ now apply mul_le_mono_nonpos_l.
+ apply mul_le_mono_nonneg_l; order.
+Qed.
+
+Lemma crossmul_le_addsquare : forall a b, 0<=a -> 0<=b -> b*a+a*b <= a*a+b*b.
+Proof.
+ assert (AUX : forall a b, 0<=a<=b -> b*a+a*b <= a*a+b*b).
+  intros a b (Ha,H).
+  destruct (le_exists_sub _ _ H) as (d & EQ & Hd).
+  rewrite EQ.
+  rewrite 2 mul_add_distr_r.
+  rewrite !add_assoc. apply add_le_mono_r.
+  rewrite add_comm. apply add_le_mono_l.
+  apply mul_le_mono_nonneg_l; trivial. order.
+ intros a b Ha Hb.
+ destruct (le_gt_cases a b).
+ apply AUX; split; order.
+ rewrite (add_comm (b*a)), (add_comm (a*a)).
+ apply AUX; split; order.
+Qed.
+
+Lemma add_square_le : forall a b, 0<=a -> 0<=b ->
+ a*a + b*b <= (a+b)*(a+b).
+Proof.
+ intros a b Ha Hb.
+ rewrite mul_add_distr_r, !mul_add_distr_l.
+ rewrite add_assoc.
+ apply add_le_mono_r.
+ rewrite <- add_assoc.
+ rewrite <- (add_0_r (a*a)) at 1.
+ apply add_le_mono_l.
+ apply add_nonneg_nonneg; now apply mul_nonneg_nonneg.
+Qed.
+
+Lemma square_add_le : forall a b, 0<=a -> 0<=b ->
+ (a+b)*(a+b) <= 2*(a*a + b*b).
+Proof.
+ intros a b Ha Hb.
+ rewrite !mul_add_distr_l, !mul_add_distr_r. nzsimpl'.
+ rewrite <- !add_assoc. apply add_le_mono_l.
+ rewrite !add_assoc. apply add_le_mono_r.
+ apply crossmul_le_addsquare; order.
 Qed.
 
 End NZMulOrderProp.
diff --git a/theories/Numbers/NatInt/NZOrder.v b/theories/Numbers/NatInt/NZOrder.v
index e2e398025..ef9057c06 100644
--- a/theories/Numbers/NatInt/NZOrder.v
+++ b/theories/Numbers/NatInt/NZOrder.v
@@ -239,9 +239,14 @@ Proof.
 apply lt_le_incl, lt_0_1.
 Qed.
 
+Theorem lt_1_2 : 1 < 2.
+Proof.
+rewrite two_succ. apply lt_succ_diag_r.
+Qed.
+
 Theorem lt_0_2 : 0 < 2.
 Proof.
-transitivity 1. apply lt_0_1. rewrite two_succ. apply lt_succ_diag_r.
+transitivity 1. apply lt_0_1. apply lt_1_2.
 Qed.
 
 Theorem le_0_2 : 0 <= 2.
@@ -256,7 +261,7 @@ Qed.
 
 (** The order tactic enriched with some knowledge of 0,1,2 *)
 
-Ltac order' := generalize lt_0_1 lt_0_2; order.
+Ltac order' := generalize lt_0_1 lt_1_2; order.
 
 
 (** More Trichotomy, decidability and double negation elimination. *)
diff --git a/theories/Numbers/NatInt/NZSqrt.v b/theories/Numbers/NatInt/NZSqrt.v
new file mode 100644
index 000000000..425e4d6b8
--- /dev/null
+++ b/theories/Numbers/NatInt/NZSqrt.v
@@ -0,0 +1,255 @@
+(************************************************************************)
+(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
+(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
+(*   \VV/  **************************************************************)
+(*    //   *      This file is distributed under the terms of the       *)
+(*         *       GNU Lesser General Public License Version 2.1        *)
+(************************************************************************)
+
+(** Square Root Function *)
+
+Require Import NZAxioms NZMulOrder.
+
+(** Interface of a sqrt function, then its specification on naturals *)
+
+Module Type Sqrt (Import T:Typ).
+ Parameters Inline sqrt : t -> t.
+End Sqrt.
+
+Module Type SqrtNotation (T:Typ)(Import NZ:Sqrt T).
+ Notation "√ x" := (sqrt x) (at level 25).
+End SqrtNotation.
+
+Module Type Sqrt' (T:Typ) := Sqrt T <+ SqrtNotation T.
+
+Module Type NZSqrtSpec (Import NZ : NZOrdAxiomsSig')(Import P : Sqrt' NZ).
+ Declare Instance sqrt_wd : Proper (eq==>eq) sqrt.
+ Axiom sqrt_spec : forall a, 0<=a -> √a * √a <= a < S (√a) * S (√a).
+End NZSqrtSpec.
+
+Module Type NZSqrt (NZ:NZOrdAxiomsSig) := Sqrt NZ <+ NZSqrtSpec NZ.
+Module Type NZSqrt' (NZ:NZOrdAxiomsSig) := Sqrt' NZ <+ NZSqrtSpec NZ.
+
+(** Derived properties of power *)
+
+Module NZSqrtProp
+ (Import NZ : NZOrdAxiomsSig')
+ (Import NZP' : NZSqrt' NZ)
+ (Import NZP : NZMulOrderProp NZ).
+
+(** First, sqrt is non-negative *)
+
+Lemma sqrt_spec_nonneg : forall a b, 0<=a ->
+ b*b <= a < S b * S b -> 0 <= b.
+Proof.
+ intros a b Ha (LE,LT).
+ destruct (le_gt_cases 0 b) as [Hb|Hb]; trivial. exfalso.
+ assert (S b * S b < b * b).
+  rewrite mul_succ_l, <- (add_0_r (b*b)).
+  apply add_lt_le_mono.
+  apply mul_lt_mono_neg_l; trivial. apply lt_succ_diag_r.
+  now apply le_succ_l.
+ order.
+Qed.
+
+Lemma sqrt_nonneg : forall a, 0<=a -> 0<=√a.
+Proof.
+ intros. now apply (sqrt_spec_nonneg a), sqrt_spec.
+Qed.
+
+(** The spec of sqrt indeed determines it *)
+
+Lemma sqrt_unique : forall a b, 0<=a -> b*b<=a<(S b)*(S b) -> √a == b.
+Proof.
+ intros a b Ha (LEb,LTb).
+ assert (0<=b) by (apply (sqrt_spec_nonneg a); try split; trivial).
+ assert (0<=√a) by now apply sqrt_nonneg.
+ destruct (sqrt_spec a Ha) as (LEa,LTa).
+ assert (b <= √a).
+  apply lt_succ_r, square_lt_simpl_nonneg; [|order].
+  now apply lt_le_incl, lt_succ_r.
+ assert (√a <= b).
+  apply lt_succ_r, square_lt_simpl_nonneg; [|order].
+  now apply lt_le_incl, lt_succ_r.
+ order.
+Qed.
+
+(** Sqrt is exact on squares *)
+
+Lemma sqrt_square : forall a, 0<=a -> √(a*a) == a.
+Proof.
+ intros a Ha.
+ apply sqrt_unique.
+ apply square_nonneg.
+ split. order.
+ apply mul_lt_mono_nonneg; trivial using lt_succ_diag_r.
+Qed.
+
+(** Sqrt is a monotone function (but not a strict one) *)
+
+Lemma sqrt_le_mono : forall a b, 0<=a<=b -> √a <= √b.
+Proof.
+ intros a b (Ha,Hab).
+ assert (Hb : 0 <= b) by order.
+ destruct (sqrt_spec a Ha) as (LE,_).
+ destruct (sqrt_spec b Hb) as (_,LT).
+ apply lt_succ_r.
+ apply square_lt_simpl_nonneg; try order.
+ now apply lt_le_incl, lt_succ_r, sqrt_nonneg.
+Qed.
+
+(** No reverse result for <=, consider for instance √2 <= √1 *)
+
+Lemma sqrt_lt_cancel : forall a b, 0<=a -> 0<=b -> √a < √b -> a < b.
+Proof.
+ intros a b Ha Hb H.
+ destruct (sqrt_spec a Ha) as (_,LT).
+ destruct (sqrt_spec b Hb) as (LE,_).
+ apply le_succ_l in H.
+ assert (S (√a) * S (√a) <= √b * √b).
+  apply mul_le_mono_nonneg; trivial.
+  now apply lt_le_incl, lt_succ_r, sqrt_nonneg.
+  now apply lt_le_incl, lt_succ_r, sqrt_nonneg.
+ order.
+Qed.
+
+(** When left side is a square, we have an equivalence for <= *)
+
+Lemma sqrt_le_square : forall a b, 0<=a -> 0<=b -> (b*b<=a <-> b <= √a).
+Proof.
+ intros a b Ha Hb. split; intros H.
+ rewrite <- (sqrt_square b); trivial.
+ apply sqrt_le_mono. split. apply square_nonneg. trivial.
+ destruct (sqrt_spec a Ha) as (LE,LT).
+ transitivity (√a * √a); trivial.
+ now apply mul_le_mono_nonneg.
+Qed.
+
+(** When right side is a square, we have an equivalence for < *)
+
+Lemma sqrt_lt_square : forall a b, 0<=a -> 0<=b -> (a<b*b <-> √a < b).
+Proof.
+ intros a b Ha Hb. split; intros H.
+ destruct (sqrt_spec a Ha) as (LE,_).
+ apply square_lt_simpl_nonneg; try order.
+ rewrite <- (sqrt_square b Hb) in H.
+ apply sqrt_lt_cancel; trivial.
+ apply square_nonneg.
+Qed.
+
+(** Sqrt and basic constants *)
+
+Lemma sqrt_0 : √0 == 0.
+Proof.
+ rewrite <- (mul_0_l 0) at 1. now apply sqrt_square.
+Qed.
+
+Lemma sqrt_1 : √1 == 1.
+Proof.
+ rewrite <- (mul_1_l 1) at 1. apply sqrt_square. order'.
+Qed.
+
+Lemma sqrt_2 : √2 == 1.
+Proof.
+ apply sqrt_unique. order'. nzsimpl. split. order'.
+ apply lt_succ_r, lt_le_incl, lt_succ_r. nzsimpl'; order.
+Qed.
+
+Lemma sqrt_lt_lin : forall a, 1<a -> √a<a.
+Proof.
+ intros a Ha. rewrite <- sqrt_lt_square; try order'.
+ rewrite <- (mul_1_r a) at 1.
+ rewrite <- mul_lt_mono_pos_l; order'.
+Qed.
+
+Lemma sqrt_le_lin : forall a, 0<=a -> √a<=a.
+Proof.
+ intros a Ha.
+ destruct (le_gt_cases a 0) as [H|H].
+ setoid_replace a with 0 by order. now rewrite sqrt_0.
+ destruct (le_gt_cases a 1) as [H'|H'].
+ rewrite <- le_succ_l, <- one_succ in H.
+ setoid_replace a with 1 by order. now rewrite sqrt_1.
+ now apply lt_le_incl, sqrt_lt_lin.
+Qed.
+
+(** Sqrt and multiplication. *)
+
+(** Due to rounding error, we don't have the usual √(a*b) = √a*√b
+    but only lower and upper bounds. *)
+
+Lemma sqrt_mul_below : forall a b, 0<=a -> 0<=b -> √a * √b <= √(a*b).
+Proof.
+ intros a b Ha Hb.
+ assert (Ha':=sqrt_nonneg _ Ha).
+ assert (Hb':=sqrt_nonneg _ Hb).
+ apply sqrt_le_square; try now apply mul_nonneg_nonneg.
+ rewrite mul_shuffle1.
+ apply mul_le_mono_nonneg; try now apply mul_nonneg_nonneg.
+  now apply sqrt_spec.
+  now apply sqrt_spec.
+Qed.
+
+Lemma sqrt_mul_above : forall a b, 0<=a -> 0<=b ->
+ √(a*b) < S (√a) * S (√b).
+Proof.
+ intros a b Ha Hb.
+ apply sqrt_lt_square.
+ now apply mul_nonneg_nonneg.
+ apply mul_nonneg_nonneg.
+ now apply lt_le_incl, lt_succ_r, sqrt_nonneg.
+ now apply lt_le_incl, lt_succ_r, sqrt_nonneg.
+ rewrite mul_shuffle1.
+ apply mul_lt_mono_nonneg; trivial; now apply sqrt_spec.
+Qed.
+
+(** And we can't find better approximations in general.
+    - The lower bound is exact for squares
+    - Concerning the upper bound, for any c>0, take a=b=c*c-1,
+      then √(a*b) = c*c -1 while S √a = S √b = c
+*)
+
+(** Sqrt and addition *)
+
+Lemma sqrt_add_le : forall a b, 0<=a -> 0<=b -> √(a+b) <= √a + √b.
+Proof.
+ intros a b Ha Hb.
+ assert (Ha':=sqrt_nonneg _ Ha).
+ assert (Hb':=sqrt_nonneg _ Hb).
+ rewrite <- lt_succ_r.
+ apply sqrt_lt_square.
+  now apply add_nonneg_nonneg.
+  now apply lt_le_incl, lt_succ_r, add_nonneg_nonneg.
+ destruct (sqrt_spec a Ha) as (_,LTa).
+ destruct (sqrt_spec b Hb) as (_,LTb).
+ revert LTa LTb. nzsimpl. rewrite 3 lt_succ_r.
+ intros LTa LTb.
+ assert (H:=add_le_mono _ _ _ _ LTa LTb).
+ etransitivity; [eexact H|]. clear LTa LTb H.
+ rewrite <- (add_assoc _ (√a) (√a)).
+ rewrite <- (add_assoc _ (√b) (√b)).
+ rewrite add_shuffle1.
+ rewrite <- (add_assoc _ (√a + √b)).
+ rewrite (add_shuffle1 (√a) (√b)).
+ apply add_le_mono_r.
+ now apply add_square_le.
+Qed.
+
+(** convexity inequality for sqrt: sqrt of middle is above middle
+    of square roots. *)
+
+Lemma add_sqrt_le : forall a b, 0<=a -> 0<=b -> √a + √b <= √(2*(a+b)).
+Proof.
+ intros a b Ha Hb.
+ assert (Ha':=sqrt_nonneg _ Ha).
+ assert (Hb':=sqrt_nonneg _ Hb).
+ apply sqrt_le_square.
+  apply mul_nonneg_nonneg. order'. now apply add_nonneg_nonneg.
+  now apply add_nonneg_nonneg.
+ transitivity (2*(√a*√a + √b*√b)).
+ now apply square_add_le.
+ apply mul_le_mono_nonneg_l. order'.
+ apply add_le_mono; now apply sqrt_spec.
+Qed.
+
+End NZSqrtProp.
diff --git a/theories/Numbers/Natural/Abstract/NAxioms.v b/theories/Numbers/Natural/Abstract/NAxioms.v
index 66ff2ded5..b360c016f 100644
--- a/theories/Numbers/Natural/Abstract/NAxioms.v
+++ b/theories/Numbers/Natural/Abstract/NAxioms.v
@@ -8,7 +8,7 @@
 (*                      Evgeny Makarov, INRIA, 2007                     *)
 (************************************************************************)
 
-Require Export NZAxioms NZPow NZDiv.
+Require Export NZAxioms NZPow NZSqrt NZDiv.
 
 (** From [NZ], we obtain natural numbers just by stating that [pred 0] == 0 *)
 
@@ -32,7 +32,7 @@ Module Type Parity (Import N : NAxiomsMiniSig').
  Axiom odd_spec : forall n, odd n = true <-> Odd n.
 End Parity.
 
-(** Power function : NZPow is enough *)
+(** For Power and Sqrt functions : NZPow and NZSqrt are enough *)
 
 (** Division Function : we reuse NZDiv.DivMod and NZDiv.NZDivCommon,
     and add to that a N-specific constraint. *)
@@ -45,10 +45,12 @@ End NDivSpecific.
 (** We now group everything together. *)
 
 Module Type NAxiomsSig := NAxiomsMiniSig <+ HasCompare <+ Parity
-  <+ NZPow.NZPow <+ DivMod <+ NZDivCommon <+ NDivSpecific.
+  <+ NZPow.NZPow <+ NZSqrt.NZSqrt
+  <+ DivMod <+ NZDivCommon <+ NDivSpecific.
 
 Module Type NAxiomsSig' := NAxiomsMiniSig' <+ HasCompare <+ Parity
-  <+ NZPow.NZPow' <+ DivMod' <+ NZDivCommon <+ NDivSpecific.
+  <+ NZPow.NZPow' <+ NZSqrt.NZSqrt'
+  <+ DivMod' <+ NZDivCommon <+ NDivSpecific.
 
 
 (** It could also be interesting to have a constructive recursor function. *)
diff --git a/theories/Numbers/Natural/Abstract/NProperties.v b/theories/Numbers/Natural/Abstract/NProperties.v
index c1977f353..fc8f27ddc 100644
--- a/theories/Numbers/Natural/Abstract/NProperties.v
+++ b/theories/Numbers/Natural/Abstract/NProperties.v
@@ -7,9 +7,9 @@
 (************************************************************************)
 
 Require Export NAxioms.
-Require Import NMaxMin NParity NPow NDiv.
+Require Import NMaxMin NParity NPow NSqrt NDiv.
 
 (** This functor summarizes all known facts about N. *)
 
 Module Type NProp (N:NAxiomsSig) :=
- NMaxMinProp N <+ NParityProp N <+ NPowProp N <+ NDivProp N.
+ NMaxMinProp N <+ NParityProp N <+ NPowProp N <+ NSqrtProp N <+ NDivProp N.
diff --git a/theories/Numbers/Natural/Abstract/NSqrt.v b/theories/Numbers/Natural/Abstract/NSqrt.v
new file mode 100644
index 000000000..d5916bdc2
--- /dev/null
+++ b/theories/Numbers/Natural/Abstract/NSqrt.v
@@ -0,0 +1,64 @@
+(************************************************************************)
+(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
+(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
+(*   \VV/  **************************************************************)
+(*    //   *      This file is distributed under the terms of the       *)
+(*         *       GNU Lesser General Public License Version 2.1        *)
+(************************************************************************)
+
+(** Properties of Square Root Function *)
+
+Require Import NAxioms NSub NZSqrt.
+
+Module NSqrtProp (Import N : NAxiomsSig')(Import NS : NSubProp N).
+
+ Module Import NZSqrtP := NZSqrtProp N N NS.
+
+ Ltac auto' := trivial; try rewrite <- neq_0_lt_0; auto using le_0_l.
+ Ltac wrap l := intros; apply l; auto'.
+
+ (** We redefine NZSqrt's results, without the non-negative hyps *)
+
+Lemma sqrt_spec' : forall a, √a*√a <= a < S (√a) * S (√a).
+Proof. wrap sqrt_spec. Qed.
+
+Lemma sqrt_unique : forall a b, b*b<=a<(S b)*(S b) -> √a == b.
+Proof. wrap sqrt_unique. Qed.
+
+Lemma sqrt_square : forall a, √(a*a) == a.
+Proof. wrap sqrt_square. Qed.
+
+Lemma sqrt_le_mono : forall a b, a<=b -> √a <= √b.
+Proof. wrap sqrt_le_mono. Qed.
+
+Lemma sqrt_lt_cancel : forall a b, √a < √b -> a < b.
+Proof. wrap sqrt_lt_cancel. Qed.
+
+Lemma sqrt_le_square : forall a b, b*b<=a <-> b <= √a.
+Proof. wrap sqrt_le_square. Qed.
+
+Lemma sqrt_lt_square : forall a b, a<b*b <-> √a < b.
+Proof. wrap sqrt_lt_square. Qed.
+
+Definition sqrt_0 := sqrt_0.
+Definition sqrt_1 := sqrt_1.
+Definition sqrt_2 := sqrt_2.
+
+Definition sqrt_lt_lin : forall a, 1<a -> √a<a := sqrt_lt_lin.
+
+Lemma sqrt_le_lin : forall a, 0<=a -> √a<=a.
+Proof. wrap sqrt_le_lin. Qed.
+
+Lemma sqrt_mul_below : forall a b, √a * √b <= √(a*b).
+Proof. wrap sqrt_mul_below. Qed.
+
+Lemma sqrt_mul_above : forall a b, √(a*b) < S (√a) * S (√b).
+Proof. wrap sqrt_mul_above. Qed.
+
+Lemma sqrt_add_le : forall a b, √(a+b) <= √a + √b.
+Proof. wrap sqrt_add_le. Qed.
+
+Lemma add_sqrt_le : forall a b, √a + √b <= √(2*(a+b)).
+Proof. wrap add_sqrt_le. Qed.
+
+End NSqrtProp.
diff --git a/theories/Numbers/Natural/BigN/NMake.v b/theories/Numbers/Natural/BigN/NMake.v
index 9e6e4b609..ec0fa89bf 100644
--- a/theories/Numbers/Natural/BigN/NMake.v
+++ b/theories/Numbers/Natural/BigN/NMake.v
@@ -738,11 +738,18 @@ Module Make (W0:CyclicType) <: NType.
  Lemma sqrt_fold : sqrt = iter_t sqrtn.
  Proof. red_t; reflexivity. Qed.
 
- Theorem spec_sqrt: forall x, [sqrt x] ^ 2 <= [x] < ([sqrt x] + 1) ^ 2.
+ Theorem spec_sqrt_aux: forall x, [sqrt x] ^ 2 <= [x] < ([sqrt x] + 1) ^ 2.
  Proof.
   intros x. rewrite sqrt_fold. destr_t x as (n,x). exact (ZnZ.spec_sqrt x).
  Qed.
 
+ Theorem spec_sqrt: forall x, [sqrt x] = Zsqrt [x].
+ Proof.
+  intros x.
+  symmetry. apply Z.sqrt_unique. apply spec_pos.
+  rewrite <- ! Zpower_2. apply spec_sqrt_aux.
+ Qed.
+
  (** * Power *)
 
  Fixpoint pow_pos (x:t)(p:positive) : t :=
diff --git a/theories/Numbers/Natural/Binary/NBinary.v b/theories/Numbers/Natural/Binary/NBinary.v
index d1217d407..348eee5ed 100644
--- a/theories/Numbers/Natural/Binary/NBinary.v
+++ b/theories/Numbers/Natural/Binary/NBinary.v
@@ -8,7 +8,7 @@
 (*                      Evgeny Makarov, INRIA, 2007                     *)
 (************************************************************************)
 
-Require Import BinPos Ndiv_def.
+Require Import BinPos Ndiv_def Nsqrt_def.
 Require Export BinNat.
 Require Import NAxioms NProperties.
 
@@ -167,6 +167,11 @@ Definition pow_0_r := Npow_0_r.
 Definition pow_succ_r n p (H:0 <= p) := Npow_succ_r n p.
 Program Instance pow_wd : Proper (eq==>eq==>eq) Npow.
 
+(** Sqrt *)
+
+Program Instance sqrt_wd : Proper (eq==>eq) Nsqrt.
+Definition sqrt_spec n (H:0<=n) := Nsqrt_spec n.
+
 (** The instantiation of operations.
     Placing them at the very end avoids having indirections in above lemmas. *)
 
@@ -192,6 +197,7 @@ Definition modulo := Nmod.
 Definition pow := Npow.
 Definition even := Neven.
 Definition odd := Nodd.
+Definition sqrt := Nsqrt.
 
 Include NProp
  <+ UsualMinMaxLogicalProperties <+ UsualMinMaxDecProperties.
diff --git a/theories/Numbers/Natural/Peano/NPeano.v b/theories/Numbers/Natural/Peano/NPeano.v
index bbf4f5cd7..b91b43e31 100644
--- a/theories/Numbers/Natural/Peano/NPeano.v
+++ b/theories/Numbers/Natural/Peano/NPeano.v
@@ -9,7 +9,7 @@
 (************************************************************************)
 
 Require Import
- Bool Peano Peano_dec Compare_dec Plus Minus Le EqNat NAxioms NProperties.
+ Bool Peano Peano_dec Compare_dec Plus Minus Le Lt EqNat NAxioms NProperties.
 
 (** Functions not already defined *)
 
@@ -134,6 +134,75 @@ Proof.
  apply le_n_S, le_minus.
 Qed.
 
+(** Square root *)
+
+(** The following square root function is linear (and tail-recursive).
+  With Peano representation, we can't do better. For faster algorithm,
+  see Psqrt/Zsqrt/Nsqrt...
+
+  We search the square root of n = k + p^2 + (q - r)
+  with q = 2p and 0<=r<=q. We starts with p=q=r=0, hence
+  looking for the square root of n = k. Then we progressively
+  decrease k and r. When k = S k' and r=0, it means we can use (S p)
+  as new sqrt candidate, since (S k')+p^2+2p = k'+(S p)^2.
+  When k reaches 0, we have found the biggest p^2 square contained
+  in n, hence the square root of n is p.
+*)
+
+Fixpoint sqrt_iter k p q r :=
+  match k with
+    | O => p
+    | S k' => match r with
+                | O => sqrt_iter k' (S p) (S (S q)) (S (S q))
+                | S r' => sqrt_iter k' p q r'
+              end
+  end.
+
+Definition sqrt n := sqrt_iter n 0 0 0.
+
+Lemma sqrt_iter_spec : forall k p q r,
+ q = p+p -> r<=q ->
+ let s := sqrt_iter k p q r in
+ s*s <= k + p*p + (q - r) < (S s)*(S s).
+Proof.
+ induction k.
+ (* k = 0 *)
+ simpl; intros p q r Hq Hr.
+ split.
+ apply le_plus_l.
+ apply le_lt_n_Sm.
+ rewrite <- mult_n_Sm.
+ rewrite plus_assoc, (plus_comm p), <- plus_assoc.
+ apply plus_le_compat; trivial.
+ rewrite <- Hq. apply le_minus.
+ (* k = S k' *)
+ destruct r.
+ (* r = 0 *)
+ intros Hq _.
+ replace (S k + p*p + (q-0)) with (k + (S p)*(S p) + (S (S q) - S (S q))).
+ apply IHk.
+ simpl. rewrite <- plus_n_Sm. congruence.
+ auto with arith.
+ rewrite minus_diag, <- minus_n_O, <- plus_n_O. simpl.
+ rewrite <- plus_n_Sm; f_equal. rewrite <- plus_assoc; f_equal.
+ rewrite <- mult_n_Sm, (plus_comm p), <- plus_assoc. congruence.
+ (* r = S r' *)
+ intros Hq Hr.
+ replace (S k + p*p + (q-S r)) with (k + p*p + (q - r)).
+ apply IHk; auto with arith.
+ simpl. rewrite plus_n_Sm. f_equal. rewrite minus_Sn_m; auto.
+Qed.
+
+Lemma sqrt_spec : forall n,
+ let s := sqrt n in s*s <= n < S s * S s.
+Proof.
+ intros.
+ replace n with (n + 0*0 + (0-0)).
+ apply sqrt_iter_spec; auto.
+ simpl. now rewrite <- 2 plus_n_O.
+Qed.
+
+
 (** * Implementation of [NAxiomsSig] by [nat] *)
 
 Module Nat
@@ -297,10 +366,14 @@ Definition odd := odd.
 Definition even_spec := even_spec.
 Definition odd_spec := odd_spec.
 
-Definition pow := pow.
 Program Instance pow_wd : Proper (eq==>eq==>eq) pow.
 Definition pow_0_r := pow_0_r.
 Definition pow_succ_r := pow_succ_r.
+Definition pow := pow.
+
+Definition sqrt_spec a (Ha:0<=a) := sqrt_spec a.
+Program Instance sqrt_wd : Proper (eq==>eq) sqrt.
+Definition sqrt := sqrt.
 
 Definition div := div.
 Definition modulo := modulo.
diff --git a/theories/Numbers/Natural/SpecViaZ/NSig.v b/theories/Numbers/Natural/SpecViaZ/NSig.v
index 7cf3e7046..703897eba 100644
--- a/theories/Numbers/Natural/SpecViaZ/NSig.v
+++ b/theories/Numbers/Natural/SpecViaZ/NSig.v
@@ -77,7 +77,7 @@ Module Type NType.
  Parameter spec_pow_pos: forall x n, [pow_pos x n] = [x] ^ Zpos n.
  Parameter spec_pow_N: forall x n, [pow_N x n] = [x] ^ Z_of_N n.
  Parameter spec_pow: forall x n, [pow x n] = [x] ^ [n].
- Parameter spec_sqrt: forall x, [sqrt x] ^ 2 <= [x] < ([sqrt x] + 1) ^ 2.
+ Parameter spec_sqrt: forall x, [sqrt x] = Zsqrt [x].
  Parameter spec_log2_0: forall x, [x] = 0 -> [log2 x] = 0.
  Parameter spec_log2: forall x, [x]<>0 -> 2^[log2 x] <= [x] < 2^([log2 x]+1).
  Parameter spec_div_eucl: forall x y,
diff --git a/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v b/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
index e1dc5349b..f072fc24a 100644
--- a/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
+++ b/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
@@ -1,4 +1,5 @@
 (************************************************************************)
+
 (*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
 (* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
 (*   \VV/  **************************************************************)
@@ -14,7 +15,7 @@ Module NTypeIsNAxioms (Import N : NType').
 
 Hint Rewrite
  spec_0 spec_1 spec_2 spec_succ spec_add spec_mul spec_pred spec_sub
- spec_div spec_modulo spec_gcd spec_compare spec_eq_bool
+ spec_div spec_modulo spec_gcd spec_compare spec_eq_bool spec_sqrt
  spec_max spec_min spec_pow_pos spec_pow_N spec_pow spec_even spec_odd
  : nsimpl.
 Ltac nsimpl := autorewrite with nsimpl.
@@ -219,6 +220,16 @@ Proof.
  intros. now zify.
 Qed.
 
+(** Sqrt *)
+
+Program Instance sqrt_wd : Proper (eq==>eq) sqrt.
+
+Lemma sqrt_spec : forall n, 0<=n ->
+ (sqrt n)*(sqrt n) <= n /\ n < (succ (sqrt n))*(succ (sqrt n)).
+Proof.
+ intros n. zify. apply Zsqrt_spec.
+Qed.
+
 (** Even / Odd *)
 
 Definition Even n := exists m, n == 2*m.
diff --git a/theories/Numbers/vo.itarget b/theories/Numbers/vo.itarget
index a54e85872..7867b8caa 100644
--- a/theories/Numbers/vo.itarget
+++ b/theories/Numbers/vo.itarget
@@ -48,6 +48,7 @@ NatInt/NZProperties.vo
 NatInt/NZDomain.vo
 NatInt/NZDiv.vo
 NatInt/NZPow.vo
+NatInt/NZSqrt.vo
 Natural/Abstract/NAddOrder.vo
 Natural/Abstract/NAdd.vo
 Natural/Abstract/NAxioms.vo
@@ -63,6 +64,7 @@ Natural/Abstract/NDiv.vo
 Natural/Abstract/NMaxMin.vo
 Natural/Abstract/NParity.vo
 Natural/Abstract/NPow.vo
+Natural/Abstract/NSqrt.vo
 Natural/BigN/BigN.vo
 Natural/BigN/Nbasic.vo
 Natural/BigN/NMake_gen.vo
diff --git a/theories/ZArith/ZArith.v b/theories/ZArith/ZArith.v
index ece3475dc..93121d48f 100644
--- a/theories/ZArith/ZArith.v
+++ b/theories/ZArith/ZArith.v
@@ -13,7 +13,7 @@ Require Export ZArith_base.
 (** Extra modules using [Omega] or [Ring]. *)
 
 Require Export Zcomplements.
-Require Export Zsqrt.
+Require Export Zsqrt_def.
 Require Export Zpower.
 Require Export Zdiv.
 Require Export Zlogarithm.
diff --git a/theories/ZArith/Zsqrt.v b/theories/ZArith/Zsqrt_compat.v
index 20fa6bb5a..ce46aa6d4 100644
--- a/theories/ZArith/Zsqrt.v
+++ b/theories/ZArith/Zsqrt_compat.v
@@ -11,6 +11,16 @@ Require Import Omega.
 Require Export ZArith_base.
 Open Local Scope Z_scope.
 
+(**    !! This file is deprecated !!
+
+    Please use rather Zsqrt_def.Zsqrt (or Zsqrtrem).
+    Unlike here, proofs there are fully separated from functions.
+    Some equivalence proofs between the old and the new versions
+    can be found below. A Require Import ZArith provides by default
+    the new versions.
+
+*)
+
 (**********************************************************************)
 (** Definition and properties of square root on Z *)
 
@@ -211,3 +221,12 @@ Proof.
 Qed.
 
 
+(** Equivalence between Zsqrt_plain and [Zsqrt_def.Zsqrt] *)
+
+Lemma Zsqrt_equiv : forall n, Zsqrt_plain n = Zsqrt_def.Zsqrt n.
+Proof.
+ intros. destruct (Z_le_gt_dec 0 n).
+ symmetry. apply Z.sqrt_unique; trivial.
+ now apply Zsqrt_interval.
+ now destruct n.
+Qed.
\ No newline at end of file
diff --git a/theories/ZArith/Zsqrt_def.v b/theories/ZArith/Zsqrt_def.v
new file mode 100644
index 000000000..15f4e5275
--- /dev/null
+++ b/theories/ZArith/Zsqrt_def.v
@@ -0,0 +1,60 @@
+(************************************************************************)
+(*  v      *   The Coq Proof Assistant  /  The Coq Development Team     *)
+(* <O___,, *   INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010     *)
+(*   \VV/  **************************************************************)
+(*    //   *      This file is distributed under the terms of the       *)
+(*         *       GNU Lesser General Public License Version 2.1        *)
+(************************************************************************)
+
+(** Definition of a square root function for Z. *)
+
+Require Import BinPos BinInt Psqrt.
+
+Open Scope Z_scope.
+
+Definition Zsqrtrem n :=
+ match n with
+  | 0 => (0, 0)
+  | Zpos p =>
+    match Psqrtrem p with
+     | (s, IsPos r) => (Zpos s, Zpos r)
+     | (s, _) => (Zpos s, 0)
+    end
+  | Zneg _ => (0,0)
+ end.
+
+Definition Zsqrt n :=
+ match n with
+  | 0 => 0
+  | Zpos p => Zpos (Psqrt p)
+  | Zneg _ => 0
+ end.
+
+Lemma Zsqrtrem_spec : forall n, 0<=n ->
+ let (s,r) := Zsqrtrem n in n = s*s + r /\ 0 <= r <= 2*s.
+Proof.
+ destruct n. now repeat split.
+ generalize (Psqrtrem_spec p). simpl.
+ destruct 1; simpl; subst; now repeat split.
+ now destruct 1.
+Qed.
+
+Lemma Zsqrt_spec : forall n, 0<=n ->
+ let s := Zsqrt n in s*s <= n < (Zsucc s)*(Zsucc s).
+Proof.
+ destruct n. now repeat split. unfold Zsqrt.
+ rewrite <- Zpos_succ_morphism. intros _. apply (Psqrt_spec p).
+ now destruct 1.
+Qed.
+
+Lemma Zsqrt_neg : forall n, n<0 -> Zsqrt n = 0.
+Proof.
+ intros. now destruct n.
+Qed.
+
+Lemma Zsqrtrem_sqrt : forall n, fst (Zsqrtrem n) = Zsqrt n.
+Proof.
+ destruct n; try reflexivity.
+ unfold Zsqrtrem, Zsqrt, Psqrt.
+ destruct (Psqrtrem p) as (s,r). now destruct r.
+Qed.
\ No newline at end of file
diff --git a/theories/ZArith/vo.itarget b/theories/ZArith/vo.itarget
index 59994e32f..28db6848d 100644
--- a/theories/ZArith/vo.itarget
+++ b/theories/ZArith/vo.itarget
@@ -27,5 +27,6 @@ Zorder.vo
 Zpow_def.vo
 Zpower.vo
 Zpow_facts.vo
-Zsqrt.vo
+Zsqrt_compat.vo
 Zwf.vo
+Zsqrt_def.vo