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

8 files changed, 175 insertions, 12 deletions

diff --git a/theories/Numbers/Natural/Abstract/NAxioms.v b/theories/Numbers/Natural/Abstract/NAxioms.v
index 66ff2ded54..b360c016f3 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 c1977f3533..fc8f27ddc9 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 0000000000..d5916bdc2d
--- /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 9e6e4b6092..ec0fa89bff 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 d1217d407e..348eee5ed2 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 bbf4f5cd7b..b91b43e310 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 7cf3e7046b..703897eba5 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 e1dc5349bc..f072fc24aa 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.