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

1 files changed, 123 insertions, 0 deletions

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.