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

7 files changed, 278 insertions, 37 deletions

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