Numbers: axiomatization, properties and implementations of gcd

 - For nat, we create a brand-new gcd function, structural in
   the sense of Coq, even if it's Euclid algorithm. Cool...
 - We re-organize the Zgcd that was in Znumtheory, create out of it
   files Pgcd, Ngcd_def, Zgcd_def. Proofs of correctness are revised
   in order to be much simpler (no omega, no advanced lemmas of
   Znumtheory, etc).
 - Abstract Properties NZGcd / ZGcd / NGcd could still be completed,
   for the moment they contain up to Gauss thm. We could add stuff
   about (relative) primality, relationship between gcd and div,mod,
   or stuff about parity, etc etc.
 - Znumtheory remains as it was, apart for Zgcd and correctness proofs
   gone elsewhere. We could later take advantage of ZGcd in it.
   Someday, we'll have to switch from the current Zdivide inductive,
   to Zdivide' via exists. To be continued...

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

27 files changed, 1547 insertions, 391 deletions

diff --git a/doc/stdlib/index-list.html.template b/doc/stdlib/index-list.html.template
index e0f445d03..43761b0ab 100644
--- a/doc/stdlib/index-list.html.template
+++ b/doc/stdlib/index-list.html.template
@@ -140,6 +140,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/PArith/POrderedType.v
     theories/PArith/Pminmax.v
     theories/PArith/Psqrt.v
+    theories/PArith/Pgcd.v
     (theories/PArith/PArith.v)
   </dd>
 
@@ -153,6 +154,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/NArith/Ndist.v
     theories/NArith/Ndec.v
     theories/NArith/Ndiv_def.v
+    theories/NArith/Ngcd_def.v
     theories/NArith/Nsqrt_def.v
     (theories/NArith/NArith.v)
   </dd>
@@ -186,6 +188,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/ZArith/Zlog_def.v
     theories/ZArith/Zlogarithm.v
     (theories/ZArith/ZArith.v)
+    theories/ZArith/Zgcd_def.v
     theories/ZArith/Zgcd_alt.v
     theories/ZArith/Zwf.v
     theories/ZArith/Znumtheory.v
@@ -240,6 +243,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/Numbers/NatInt/NZPow.v
     theories/Numbers/NatInt/NZSqrt.v
     theories/Numbers/NatInt/NZLog.v
+    theories/Numbers/NatInt/NZGcd.v
     </dd>
     </dt>
 
@@ -284,6 +288,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/Numbers/Natural/Abstract/NPow.v
     theories/Numbers/Natural/Abstract/NSqrt.v
     theories/Numbers/Natural/Abstract/NLog.v
+    theories/Numbers/Natural/Abstract/NGcd.v
     theories/Numbers/Natural/Abstract/NProperties.v
     theories/Numbers/Natural/Binary/NBinary.v
     theories/Numbers/Natural/Peano/NPeano.v
@@ -310,6 +315,7 @@ through the <tt>Require Import</tt> command.</p>
     theories/Numbers/Integer/Abstract/ZMaxMin.v
     theories/Numbers/Integer/Abstract/ZParity.v
     theories/Numbers/Integer/Abstract/ZPow.v
+    theories/Numbers/Integer/Abstract/ZGcd.v
     theories/Numbers/Integer/Abstract/ZProperties.v
     theories/Numbers/Integer/Abstract/ZDivEucl.v
     theories/Numbers/Integer/Abstract/ZDivFloor.v
diff --git a/theories/NArith/NArith.v b/theories/NArith/NArith.v
index 5cfd0c416..22c9012a7 100644
--- a/theories/NArith/NArith.v
+++ b/theories/NArith/NArith.v
@@ -13,6 +13,7 @@ Require Export BinNat.
 Require Export Nnat.
 Require Export Ndiv_def.
 Require Export Nsqrt_def.
+Require Export Ngcd_def.
 Require Export Ndigits.
 Require Export NArithRing.
 Require NBinary.
diff --git a/theories/NArith/Ngcd_def.v b/theories/NArith/Ngcd_def.v
new file mode 100644
index 000000000..fe5185c6b
--- /dev/null
+++ b/theories/NArith/Ngcd_def.v
@@ -0,0 +1,85 @@
+(************************************************************************)
+(*  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 BinNat Pgcd.
+
+Local Open Scope N_scope.
+
+(** * Divisibility *)
+
+Definition Ndivide p q := exists r, p*r = q.
+Notation "( p | q )" := (Ndivide p q) (at level 0) : N_scope.
+
+(** * Definition of a Pgcd function for binary natural numbers *)
+
+Definition Ngcd (a b : N) :=
+ match a, b with
+  | 0, _ => b
+  | _, 0 => a
+  | Npos p, Npos q => Npos (Pgcd p q)
+ end.
+
+(** * Generalized Gcd, also computing rests of a and b after
+    division by gcd. *)
+
+Definition Nggcd (a b : N) :=
+ match a, b with
+  | 0, _ => (b,(0,1))
+  | _, 0 => (a,(1,0))
+  | Npos p, Npos q =>
+     let '(g,(aa,bb)) := Pggcd p q in
+     (Npos g, (Npos aa, Npos bb))
+ end.
+
+(** The first component of Nggcd is Ngcd *)
+
+Lemma Nggcd_gcd : forall  a b, fst (Nggcd a b) = Ngcd a b.
+Proof.
+ intros [ |p] [ |q]; simpl; auto.
+ generalize (Pggcd_gcd p q). destruct Pggcd as (g,(aa,bb)); simpl; congruence.
+Qed.
+
+(** The other components of Nggcd are indeed the correct factors. *)
+
+Lemma Nggcd_correct_divisors : forall a b,
+  let '(g,(aa,bb)) := Nggcd a b in
+  a=g*aa /\ b=g*bb.
+Proof.
+ intros [ |p] [ |q]; simpl; auto.
+ now rewrite Pmult_1_r.
+ now rewrite Pmult_1_r.
+ generalize (Pggcd_correct_divisors p q).
+ destruct Pggcd as (g,(aa,bb)); simpl. destruct 1; split; congruence.
+Qed.
+
+(** We can use this fact to prove a part of the gcd correctness *)
+
+Lemma Ngcd_divide_l : forall a b, (Ngcd a b | a).
+Proof.
+ intros a b. rewrite <- Nggcd_gcd. generalize (Nggcd_correct_divisors a b).
+ destruct Nggcd as (g,(aa,bb)); simpl. intros (H,_). exists aa; auto.
+Qed.
+
+Lemma Ngcd_divide_r : forall a b, (Ngcd a b | b).
+Proof.
+ intros a b. rewrite <- Nggcd_gcd. generalize (Nggcd_correct_divisors a b).
+ destruct Nggcd as (g,(aa,bb)); simpl. intros (_,H). exists bb; auto.
+Qed.
+
+(** We now prove directly that gcd is the greatest amongst common divisors *)
+
+Lemma Ngcd_greatest : forall a b c, (c|a) -> (c|b) -> (c|Ngcd a b).
+Proof.
+ intros [ |p] [ |q]; simpl; auto.
+ intros [ |r]. intros (s,H). discriminate.
+ intros ([ |s],Hs) ([ |t],Ht); try discriminate; simpl in *.
+ destruct (Pgcd_greatest p q r) as (u,H).
+ exists s. now inversion Hs.
+ exists t. now inversion Ht.
+ exists (Npos u). simpl; now f_equal.
+Qed.
diff --git a/theories/NArith/vo.itarget b/theories/NArith/vo.itarget
index 1797b25e4..3e285772d 100644
--- a/theories/NArith/vo.itarget
+++ b/theories/NArith/vo.itarget
@@ -6,3 +6,4 @@ Ndist.vo
 Nnat.vo
 Ndiv_def.vo
 Nsqrt_def.vo
+Ngcd_def.vo
\ No newline at end of file
diff --git a/theories/Numbers/Integer/Abstract/ZAxioms.v b/theories/Numbers/Integer/Abstract/ZAxioms.v
index 47286c729..754c0b3d1 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 NZSqrt NZLog.
+Require Import NZPow NZSqrt NZLog NZGcd.
 
 (** We obtain integers by postulating that successor of predecessor
     is identity. *)
@@ -70,16 +70,16 @@ Module Type Parity (Import Z : ZAxiomsMiniSig').
  Axiom odd_spec : forall n, odd n = true <-> Odd n.
 End Parity.
 
-(** For pow sqrt log2, the NZ axiomatizations are enough. *)
+(** For pow sqrt log2 gcd, the NZ axiomatizations are enough. *)
 
 (** Let's group everything *)
 
 Module Type ZAxiomsSig :=
   ZAxiomsMiniSig <+ HasCompare <+ HasAbs <+ HasSgn <+ Parity
-   <+ NZPow.NZPow <+ NZSqrt.NZSqrt <+ NZLog.NZLog2.
+   <+ NZPow.NZPow <+ NZSqrt.NZSqrt <+ NZLog.NZLog2 <+ NZGcd.NZGcd.
 Module Type ZAxiomsSig' :=
   ZAxiomsMiniSig' <+ HasCompare <+ HasAbs <+ HasSgn <+ Parity
-   <+ NZPow.NZPow' <+ NZSqrt.NZSqrt' <+ NZLog.NZLog2.
+   <+ NZPow.NZPow' <+ NZSqrt.NZSqrt' <+ NZLog.NZLog2 <+ NZGcd.NZGcd'.
 
 
 (** Division is left apart, since many different flavours are available *)
diff --git a/theories/Numbers/Integer/Abstract/ZGcd.v b/theories/Numbers/Integer/Abstract/ZGcd.v
new file mode 100644
index 000000000..d58d1f1e2
--- /dev/null
+++ b/theories/Numbers/Integer/Abstract/ZGcd.v
@@ -0,0 +1,276 @@
+(************************************************************************)
+(*  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 the greatest common divisor *)
+
+Require Import ZAxioms ZMulOrder ZSgnAbs NZGcd.
+
+Module ZGcdProp
+ (Import A : ZAxiomsSig')
+ (Import B : ZMulOrderProp A)
+ (Import D : ZSgnAbsProp A B).
+
+ Include NZGcdProp A A B.
+
+(** Results concerning divisibility*)
+
+Lemma divide_opp_l : forall n m, (-n | m) <-> (n | m).
+Proof.
+ intros n m. split; intros (p,Hp); exists (-p); rewrite <- Hp.
+  now rewrite mul_opp_l, mul_opp_r.
+  now rewrite mul_opp_opp.
+Qed.
+
+Lemma divide_opp_r : forall n m, (n | -m) <-> (n | m).
+Proof.
+ intros n m. split; intros (p,Hp); exists (-p).
+  now rewrite mul_opp_r, Hp, opp_involutive.
+  now rewrite <- Hp, mul_opp_r.
+Qed.
+
+Lemma divide_abs_l : forall n m, (abs n | m) <-> (n | m).
+Proof.
+ intros n m. destruct (abs_eq_or_opp n) as [H|H]; rewrite H.
+ easy. apply divide_opp_l.
+Qed.
+
+Lemma divide_abs_r : forall n m, (n | abs m) <-> (n | m).
+Proof.
+ intros n m. destruct (abs_eq_or_opp m) as [H|H]; rewrite H.
+ easy. apply divide_opp_r.
+Qed.
+
+Lemma divide_1_r_abs : forall n, (n | 1) -> abs n == 1.
+Proof.
+ intros n Hn. apply divide_1_r_nonneg. apply abs_nonneg.
+ now apply divide_abs_l.
+Qed.
+
+Lemma divide_1_r : forall n, (n | 1) -> n==1 \/ n==-(1).
+Proof.
+ intros n (m,Hm). now apply eq_mul_1 with m.
+Qed.
+
+Lemma divide_antisym_abs : forall n m,
+ (n | m) -> (m | n) -> abs n == abs m.
+Proof.
+ intros. apply divide_antisym_nonneg; try apply abs_nonneg.
+ now apply divide_abs_l, divide_abs_r.
+ now apply divide_abs_l, divide_abs_r.
+Qed.
+
+Lemma divide_antisym : forall n m,
+ (n | m) -> (m | n) -> n == m \/ n == -m.
+Proof.
+ intros. now apply abs_eq_cases, divide_antisym_abs.
+Qed.
+
+Lemma divide_sub_r : forall n m p, (n | m) -> (n | p) -> (n | m - p).
+Proof.
+ intros n m p H H'. rewrite <- add_opp_r.
+ apply divide_add_r; trivial. now apply divide_opp_r.
+Qed.
+
+Lemma divide_add_cancel_r : forall n m p, (n | m) -> (n | m + p) -> (n | p).
+Proof.
+ intros n m p H H'. rewrite <- (add_simpl_l m p). now apply divide_sub_r.
+Qed.
+
+(** Properties of gcd *)
+
+Lemma gcd_opp_l : forall n m, gcd (-n) m == gcd n m.
+Proof.
+ intros. apply gcd_unique_alt; try apply gcd_nonneg.
+ intros. rewrite divide_opp_r. apply gcd_divide_iff.
+Qed.
+
+Lemma gcd_opp_r : forall n m, gcd n (-m) == gcd n m.
+Proof.
+ intros. now rewrite gcd_comm, gcd_opp_l, gcd_comm.
+Qed.
+
+Lemma gcd_abs_l : forall n m, gcd (abs n) m == gcd n m.
+Proof.
+ intros. destruct (abs_eq_or_opp n) as [H|H]; rewrite H.
+ easy. apply gcd_opp_l.
+Qed.
+
+Lemma gcd_abs_r : forall n m, gcd n (abs m) == gcd n m.
+Proof.
+ intros. now rewrite gcd_comm, gcd_abs_l, gcd_comm.
+Qed.
+
+Lemma gcd_0_l : forall n, gcd 0 n == abs n.
+Proof.
+ intros. rewrite <- gcd_abs_r. apply gcd_0_l_nonneg, abs_nonneg.
+Qed.
+
+Lemma gcd_0_r : forall n, gcd n 0 == abs n.
+Proof.
+ intros. now rewrite gcd_comm, gcd_0_l.
+Qed.
+
+Lemma gcd_diag : forall n, gcd n n == abs n.
+Proof.
+ intros. rewrite <- gcd_abs_l, <- gcd_abs_r.
+ apply gcd_diag_nonneg, abs_nonneg.
+Qed.
+
+Lemma gcd_add_mult_diag_r : forall n m p, gcd n (m+p*n) == gcd n m.
+Proof.
+ intros. apply gcd_unique_alt; try apply gcd_nonneg.
+ intros. rewrite gcd_divide_iff. split; intros (U,V); split; trivial.
+ apply divide_add_r; trivial. now apply divide_mul_r.
+ apply divide_add_cancel_r with (p*n); trivial.
+ now apply divide_mul_r. now rewrite add_comm.
+Qed.
+
+Lemma gcd_add_diag_r : forall n m, gcd n (m+n) == gcd n m.
+Proof.
+ intros n m. rewrite <- (mul_1_l n) at 2. apply gcd_add_mult_diag_r.
+Qed.
+
+Lemma gcd_sub_diag_r : forall n m, gcd n (m-n) == gcd n m.
+Proof.
+ intros n m. rewrite <- (mul_1_l n) at 2.
+ rewrite <- add_opp_r, <- mul_opp_l. apply gcd_add_mult_diag_r.
+Qed.
+
+Definition Bezout n m p := exists a, exists b, a*n + b*m == p.
+
+Instance Bezout_wd : Proper (eq==>eq==>eq==>iff) Bezout.
+Proof.
+ unfold Bezout. intros x x' Hx y y' Hy z z' Hz.
+ setoid_rewrite Hx. setoid_rewrite Hy. now setoid_rewrite Hz.
+Qed.
+
+Lemma bezout_1_gcd : forall n m, Bezout n m 1 -> gcd n m == 1.
+Proof.
+ intros n m (q & r & H).
+ apply gcd_unique; trivial using divide_1_l, le_0_1.
+ intros p Hn Hm.
+ rewrite <- H. apply divide_add_r; now apply divide_mul_r.
+Qed.
+
+Lemma gcd_bezout : forall n m p, gcd n m == p -> Bezout n m p.
+Proof.
+ (* First, a version restricted to natural numbers *)
+ assert (aux : forall n, 0<=n -> forall m, 0<=m -> Bezout n m (gcd n m)).
+  intros n Hn; pattern n.
+  apply strong_right_induction with (z:=0); trivial.
+  unfold Bezout. solve_predicate_wd.
+  clear n Hn. intros n Hn IHn.
+  apply le_lteq in Hn; destruct Hn as [Hn|Hn].
+  intros m Hm; pattern m.
+  apply strong_right_induction with (z:=0); trivial.
+  unfold Bezout. solve_predicate_wd.
+  clear m Hm. intros m Hm IHm.
+  destruct (lt_trichotomy n m) as [LT|[EQ|LT]].
+  (* n < m *)
+  destruct (IHm (m-n)) as (a & b & EQ).
+  apply sub_nonneg; order.
+  now apply lt_sub_pos.
+  exists (a-b). exists b.
+  rewrite gcd_sub_diag_r in EQ. rewrite <- EQ.
+  rewrite mul_sub_distr_r, mul_sub_distr_l.
+  now rewrite add_sub_assoc, add_sub_swap.
+  (* n = m *)
+  rewrite EQ. rewrite gcd_diag_nonneg; trivial.
+  exists 1. exists 0. now nzsimpl.
+  (* m < n *)
+  destruct (IHn m Hm LT n) as (a & b & EQ). order.
+  exists b. exists a. now rewrite gcd_comm, <- EQ, add_comm.
+  (* n = 0 *)
+  intros m Hm. rewrite <- Hn, gcd_0_l_nonneg; trivial.
+  exists 0. exists 1. now nzsimpl.
+ (* Then we relax the positivity condition on n *)
+ assert (aux' : forall n m, 0<=m -> Bezout n m (gcd n m)).
+  intros n m Hm.
+  destruct (le_ge_cases 0 n). now apply aux.
+  assert (Hn' : 0 <= -n) by now apply opp_nonneg_nonpos.
+  destruct (aux (-n) Hn' m Hm) as (a & b & EQ).
+  exists (-a). exists b. now rewrite <- gcd_opp_l, <- EQ, mul_opp_r, mul_opp_l.
+ (* And finally we do the same for m *)
+ intros n m p Hp. rewrite <- Hp; clear Hp.
+ destruct (le_ge_cases 0 m). now apply aux'.
+ assert (Hm' : 0 <= -m) by now apply opp_nonneg_nonpos.
+ destruct (aux' n (-m) Hm') as (a & b & EQ).
+ exists a. exists (-b). now rewrite <- gcd_opp_r, <- EQ, mul_opp_r, mul_opp_l.
+Qed.
+
+Lemma gcd_mul_mono_l :
+  forall n m p, gcd (p * n) (p * m) == abs p * gcd n m.
+Proof.
+ intros n m p.
+ apply gcd_unique.
+ apply mul_nonneg_nonneg; trivial using gcd_nonneg, abs_nonneg.
+ destruct (gcd_divide_l n m) as (q,Hq).
+ rewrite <- Hq at 2. rewrite <- (abs_sgn p) at 2.
+ rewrite mul_shuffle1. apply divide_factor_l.
+ destruct (gcd_divide_r n m) as (q,Hq).
+ rewrite <- Hq at 2. rewrite <- (abs_sgn p) at 2.
+ rewrite mul_shuffle1. apply divide_factor_l.
+ intros q H H'.
+ destruct (gcd_bezout n m (gcd n m) (eq_refl _)) as (a & b & EQ).
+ rewrite <- EQ, <- sgn_abs, mul_add_distr_l. apply divide_add_r.
+ rewrite mul_shuffle2. now apply divide_mul_l.
+ rewrite mul_shuffle2. now apply divide_mul_l.
+Qed.
+
+Lemma gcd_mul_mono_l_nonneg :
+ forall n m p, 0<=p -> gcd (p*n) (p*m) == p * gcd n m.
+Proof.
+ intros. rewrite <- (abs_eq p) at 3; trivial. apply gcd_mul_mono_l.
+Qed.
+
+Lemma gcd_mul_mono_r :
+ forall n m p, gcd (n * p) (m * p) == gcd n m * abs p.
+Proof.
+ intros n m p. now rewrite !(mul_comm _ p), gcd_mul_mono_l, mul_comm.
+Qed.
+
+Lemma gcd_mul_mono_r_nonneg :
+ forall n m p, 0<=p -> gcd (n*p) (m*p) == gcd n m * p.
+Proof.
+ intros. rewrite <- (abs_eq p) at 3; trivial. apply gcd_mul_mono_r.
+Qed.
+
+Lemma gauss : forall n m p, (n | m * p) -> gcd n m == 1 -> (n | p).
+Proof.
+ intros n m p H G.
+ destruct (gcd_bezout n m 1 G) as (a & b & EQ).
+ rewrite <- (mul_1_l p), <- EQ, mul_add_distr_r.
+ apply divide_add_r. rewrite mul_shuffle0. apply divide_factor_r.
+ rewrite <- mul_assoc. now apply divide_mul_r.
+Qed.
+
+Lemma divide_mul_split : forall n m p, n ~= 0 -> (n | m * p) ->
+ exists q, exists r, n == q*r /\ (q | m) /\ (r | p).
+Proof.
+ intros n m p Hn H.
+ assert (G := gcd_nonneg n m).
+ apply le_lteq in G; destruct G as [G|G].
+ destruct (gcd_divide_l n m) as (q,Hq).
+ exists (gcd n m). exists q.
+ split. easy.
+ split. apply gcd_divide_r.
+ destruct (gcd_divide_r n m) as (r,Hr).
+ rewrite <- Hr in H. rewrite <- Hq in H at 1.
+ rewrite <- mul_assoc in H. apply mul_divide_cancel_l in H; [|order].
+ apply gauss with r; trivial.
+ apply mul_cancel_l with (gcd n m); [order|].
+ rewrite mul_1_r.
+ rewrite <- gcd_mul_mono_l_nonneg, Hq, Hr; order.
+ symmetry in G. apply gcd_eq_0 in G. destruct G as (Hn',_); order.
+Qed.
+
+(** TODO : relation between gcd and division and modulo *)
+
+(** TODO : more about rel_prime (i.e. gcd == 1), about prime ... *)
+
+End ZGcdProp.
diff --git a/theories/Numbers/Integer/Abstract/ZProperties.v b/theories/Numbers/Integer/Abstract/ZProperties.v
index d2e962673..6fbf0f23d 100644
--- a/theories/Numbers/Integer/Abstract/ZProperties.v
+++ b/theories/Numbers/Integer/Abstract/ZProperties.v
@@ -6,10 +6,10 @@
 (*         *       GNU Lesser General Public License Version 2.1        *)
 (************************************************************************)
 
-Require Export ZAxioms ZMaxMin ZSgnAbs ZParity ZPow.
+Require Export ZAxioms ZMaxMin ZSgnAbs ZParity ZPow ZGcd.
 
 (** This functor summarizes all known facts about Z. *)
 
 Module Type ZProp (Z:ZAxiomsSig) :=
  ZMaxMinProp Z <+ ZSgnAbsProp Z <+ ZParityProp Z <+ ZPowProp Z
- <+ NZSqrt.NZSqrtProp Z Z <+ NZLog.NZLog2Prop Z Z Z.
+ <+ NZSqrt.NZSqrtProp Z Z <+ NZLog.NZLog2Prop Z Z Z <+ ZGcdProp Z.
diff --git a/theories/Numbers/Integer/Binary/ZBinary.v b/theories/Numbers/Integer/Binary/ZBinary.v
index b92b303f0..01d36854b 100644
--- a/theories/Numbers/Integer/Binary/ZBinary.v
+++ b/theories/Numbers/Integer/Binary/ZBinary.v
@@ -10,7 +10,7 @@
 
 
 Require Import ZAxioms ZProperties BinInt Zcompare Zorder ZArith_dec
- Zbool Zeven Zsqrt_def Zpow_def Zlog_def.
+ Zbool Zeven Zsqrt_def Zpow_def Zlog_def Zgcd_def.
 
 Local Open Scope Z_scope.
 
@@ -161,6 +161,15 @@ Definition log2_spec := Zlog2_spec.
 Definition log2_nonpos := Zlog2_nonpos.
 Definition log2 := Zlog2.
 
+(** Gcd *)
+
+Definition gcd_divide_l := Zgcd_divide_l.
+Definition gcd_divide_r := Zgcd_divide_r.
+Definition gcd_greatest := Zgcd_greatest.
+Definition gcd_nonneg := Zgcd_nonneg.
+Definition divide := Zdivide'.
+Definition gcd := Zgcd.
+
 (** We define [eq] only here to avoid refering to this [eq] above. *)
 
 Definition eq := (@eq Z).
diff --git a/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v b/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
index 6689875cb..999e7cd03 100644
--- a/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
+++ b/theories/Numbers/Integer/SpecViaZ/ZSigZAxioms.v
@@ -20,7 +20,7 @@ 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_sqrt
  spec_compare spec_eq_bool spec_max spec_min spec_abs spec_sgn
- spec_pow spec_log2 spec_even spec_odd
+ spec_pow spec_log2 spec_even spec_odd spec_gcd
  : zsimpl.
 
 Ltac zsimpl := autorewrite with zsimpl.
@@ -349,6 +349,38 @@ Proof.
 intros a b. zify. intros. apply Z_mod_neg; auto with zarith.
 Qed.
 
+(** Gcd *)
+
+Definition divide n m := exists p, n*p == m.
+Local Notation "( x | y )" := (divide x y) (at level 0).
+
+Lemma spec_divide : forall n m, (n|m) <-> Zdivide' [n] [m].
+Proof.
+ intros n m. split.
+ intros (p,H). exists [p]. revert H; now zify.
+ intros (z,H). exists (of_Z z). now zify.
+Qed.
+
+Lemma gcd_divide_l : forall n m, (gcd n m | n).
+Proof.
+ intros n m. apply spec_divide. zify. apply Zgcd_divide_l.
+Qed.
+
+Lemma gcd_divide_r : forall n m, (gcd n m | m).
+Proof.
+ intros n m. apply spec_divide. zify. apply Zgcd_divide_r.
+Qed.
+
+Lemma gcd_greatest : forall n m p, (p|n) -> (p|m) -> (p|gcd n m).
+Proof.
+ intros n m p. rewrite !spec_divide. zify. apply Zgcd_greatest.
+Qed.
+
+Lemma gcd_nonneg : forall n m, 0 <= gcd n m.
+Proof.
+ intros. zify. apply Zgcd_nonneg.
+Qed.
+
 End ZTypeIsZAxioms.
 
 Module ZType_ZAxioms (Z : ZType)
diff --git a/theories/Numbers/NatInt/NZGcd.v b/theories/Numbers/NatInt/NZGcd.v
new file mode 100644
index 000000000..9c022646b
--- /dev/null
+++ b/theories/Numbers/NatInt/NZGcd.v
@@ -0,0 +1,299 @@
+(************************************************************************)
+(*  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        *)
+(************************************************************************)
+
+(** Greatest Common Divisor *)
+
+Require Import NZAxioms NZMulOrder.
+
+(** Interface of a gcd function, then its specification on naturals *)
+
+Module Type Gcd (Import A : Typ).
+ Parameters Inline gcd : t -> t -> t.
+End Gcd.
+
+Module Type NZGcdSpec (A : NZOrdAxiomsSig')(B : Gcd A).
+ Import A B.
+ Definition divide n m := exists p, n*p == m.
+ Local Notation "( n | m )" := (divide n m) (at level 0).
+ Axiom gcd_divide_l : forall n m, (gcd n m | n).
+ Axiom gcd_divide_r : forall n m, (gcd n m | m).
+ Axiom gcd_greatest : forall n m p, (p | n) -> (p | m) -> (p | gcd n m).
+ Axiom gcd_nonneg : forall n m, 0 <= gcd n m.
+End NZGcdSpec.
+
+Module Type DivideNotation (A:NZOrdAxiomsSig')(B:Gcd A)(C:NZGcdSpec A B).
+ Import A B C.
+ Notation "( n | m )" := (divide n m) (at level 0).
+End DivideNotation.
+
+Module Type NZGcd (A : NZOrdAxiomsSig) := Gcd A <+ NZGcdSpec A.
+Module Type NZGcd' (A : NZOrdAxiomsSig) :=
+ Gcd A <+ NZGcdSpec A <+ DivideNotation A.
+
+(** Derived properties of gcd *)
+
+Module NZGcdProp
+ (Import A : NZOrdAxiomsSig')
+ (Import B : NZGcd' A)
+ (Import C : NZMulOrderProp A).
+
+(** Results concerning divisibility*)
+
+Instance divide_wd : Proper (eq==>eq==>iff) divide.
+Proof.
+ unfold divide. intros x x' Hx y y' Hy.
+ setoid_rewrite Hx. setoid_rewrite Hy. easy.
+Qed.
+
+Lemma divide_1_l : forall n, (1 | n).
+Proof.
+ intros n. exists n. now nzsimpl.
+Qed.
+
+Lemma divide_0_r : forall n, (n | 0).
+Proof.
+ intros n. exists 0. now nzsimpl.
+Qed.
+
+Hint Rewrite divide_1_l divide_0_r : nz.
+
+Lemma divide_0_l : forall n, (0 | n) -> n==0.
+Proof.
+ intros n (m,Hm). revert Hm. now nzsimpl.
+Qed.
+
+Lemma eq_mul_1_nonneg : forall n m,
+ 0<=n -> n*m == 1 -> n==1 /\ m==1.
+Proof.
+ intros n m Hn H.
+ apply le_lteq in Hn. destruct Hn as [Hn|Hn].
+ destruct (lt_ge_cases m 0) as [Hm|Hm].
+ generalize (mul_pos_neg n m Hn Hm). order'.
+ apply le_lteq in Hm. destruct Hm as [Hm|Hm].
+ apply le_succ_l in Hn. rewrite <- one_succ in Hn.
+ apply le_lteq in Hn. destruct Hn as [Hn|Hn].
+ generalize (lt_1_mul_pos n m Hn Hm). order.
+ rewrite <- Hn, mul_1_l in H. now split.
+ rewrite <- Hm, mul_0_r in H. order'.
+ rewrite <- Hn, mul_0_l in H. order'.
+Qed.
+
+Lemma divide_1_r_nonneg : forall n, 0<=n -> (n | 1) -> n==1.
+Proof.
+ intros n Hn (m,Hm). now apply (eq_mul_1_nonneg n m).
+Qed.
+
+Lemma divide_refl : forall n, (n | n).
+Proof.
+ intros n. exists 1. now nzsimpl.
+Qed.
+
+Lemma divide_trans : forall n m p, (n | m) -> (m | p) -> (n | p).
+Proof.
+ intros n m p (q,Hq) (r,Hr). exists (q*r).
+ now rewrite mul_assoc, Hq.
+Qed.
+
+Instance divide_reflexive : Reflexive divide := divide_refl.
+Instance divide_transitive : Transitive divide := divide_trans.
+
+(** Due to sign, no general antisymmetry result *)
+
+Lemma divide_antisym_nonneg : forall n m,
+ 0<=n -> 0<=m -> (n | m) -> (m | n) -> n == m.
+Proof.
+ intros n m Hn Hm (q,Hq) (r,Hr).
+ apply le_lteq in Hn. destruct Hn as [Hn|Hn].
+ destruct (lt_ge_cases q 0) as [Hq'|Hq'].
+ generalize (mul_pos_neg n q Hn Hq'). order.
+ rewrite <- Hq, <- mul_assoc in Hr.
+ apply mul_id_r in Hr; [|order].
+ destruct (eq_mul_1_nonneg q r) as [H _]; trivial.
+ now rewrite H, mul_1_r in Hq.
+ rewrite <- Hn, mul_0_l in Hq. now rewrite <- Hn.
+Qed.
+
+Lemma mul_divide_mono_l : forall n m p, (n | m) -> (p * n | p * m).
+Proof.
+ intros n m p (q,Hq). exists q. now rewrite <- mul_assoc, Hq.
+Qed.
+
+Lemma mul_divide_mono_r : forall n m p, (n | m) -> (n * p | m * p).
+Proof.
+ intros n m p (q,Hq). exists q. now rewrite mul_shuffle0, Hq.
+Qed.
+
+Lemma mul_divide_cancel_l : forall n m p, p ~= 0 ->
+ ((p * n | p * m) <-> (n | m)).
+Proof.
+ intros n m p Hp. split.
+ intros (q,Hq). exists q. now rewrite <- mul_assoc, mul_cancel_l in Hq.
+ apply mul_divide_mono_l.
+Qed.
+
+Lemma mul_divide_cancel_r : forall n m p, p ~= 0 ->
+ ((n * p | m * p) <-> (n | m)).
+Proof.
+ intros. rewrite 2 (mul_comm _ p). now apply mul_divide_cancel_l.
+Qed.
+
+Lemma divide_add_r : forall n m p, (n | m) -> (n | p) -> (n | m + p).
+Proof.
+ intros n m p (q,Hq) (r,Hr). exists (q+r).
+ now rewrite mul_add_distr_l, Hq, Hr.
+Qed.
+
+Lemma divide_mul_l : forall n m p, (n | m) -> (n | m * p).
+Proof.
+ intros n m p (q,Hq). exists (q*p). now rewrite mul_assoc, Hq.
+Qed.
+
+Lemma divide_mul_r : forall n m p, (n | p) -> (n | m * p).
+Proof.
+ intros n m p. rewrite mul_comm. apply divide_mul_l.
+Qed.
+
+Lemma divide_factor_l : forall n m, (n | n * m).
+Proof.
+ intros. apply divide_mul_l, divide_refl.
+Qed.
+
+Lemma divide_factor_r : forall n m, (n | m * n).
+Proof.
+ intros. apply divide_mul_r, divide_refl.
+Qed.
+
+Lemma divide_pos_le : forall n m, 0 < m -> (n | m) -> n <= m.
+Proof.
+ intros n m Hm (q,Hq).
+ destruct (le_gt_cases n 0) as [Hn|Hn]. order.
+ rewrite <- Hq.
+ destruct (lt_ge_cases q 0) as [Hq'|Hq'].
+ generalize (mul_pos_neg n q Hn Hq'). order.
+ apply le_lteq in Hq'. destruct Hq' as [Hq'|Hq'].
+ rewrite <- (mul_1_r n) at 1. apply mul_le_mono_pos_l; trivial.
+ now rewrite one_succ, le_succ_l.
+ rewrite <- Hq', mul_0_r in Hq. order.
+Qed.
+
+(** Basic properties of gcd *)
+
+Lemma gcd_unique : forall n m p,
+ 0<=p -> (p|n) -> (p|m) ->
+ (forall q, (q|n) -> (q|m) -> (q|p)) ->
+ gcd n m == p.
+Proof.
+ intros n m p Hp Hn Hm H.
+ apply divide_antisym_nonneg; trivial. apply gcd_nonneg.
+ apply H. apply gcd_divide_l. apply gcd_divide_r.
+ now apply gcd_greatest.
+Qed.
+
+Instance gcd_wd : Proper (eq==>eq==>eq) gcd.
+Proof.
+ intros x x' Hx y y' Hy.
+ apply gcd_unique.
+ apply gcd_nonneg.
+ rewrite Hx. apply gcd_divide_l.
+ rewrite Hy. apply gcd_divide_r.
+ intro. rewrite Hx, Hy. apply gcd_greatest.
+Qed.
+
+Lemma gcd_divide_iff : forall n m p,
+  (p | gcd n m) <-> (p | n) /\ (p | m).
+Proof.
+ intros. split. split.
+ transitivity (gcd n m); trivial using gcd_divide_l.
+ transitivity (gcd n m); trivial using gcd_divide_r.
+ intros (H,H'). now apply gcd_greatest.
+Qed.
+
+Lemma gcd_unique_alt : forall n m p, 0<=p ->
+ (forall q, (q|p) <-> (q|n) /\ (q|m)) ->
+ gcd n m == p.
+Proof.
+ intros n m p Hp H.
+ apply gcd_unique; trivial.
+ apply -> H. apply divide_refl.
+ apply -> H. apply divide_refl.
+ intros. apply H. now split.
+Qed.
+
+Lemma gcd_comm : forall n m, gcd n m == gcd m n.
+Proof.
+ intros. apply gcd_unique_alt; try apply gcd_nonneg.
+ intros. rewrite and_comm. apply gcd_divide_iff.
+Qed.
+
+Lemma gcd_assoc : forall n m p, gcd n (gcd m p) == gcd (gcd n m) p.
+Proof.
+ intros. apply gcd_unique_alt; try apply gcd_nonneg.
+ intros. now rewrite !gcd_divide_iff, and_assoc.
+Qed.
+
+Lemma gcd_0_l_nonneg : forall n, 0<=n -> gcd 0 n == n.
+Proof.
+ intros. apply gcd_unique; trivial.
+ apply divide_0_r.
+ apply divide_refl.
+Qed.
+
+Lemma gcd_0_r_nonneg : forall n, 0<=n -> gcd n 0 == n.
+Proof.
+ intros. now rewrite gcd_comm, gcd_0_l_nonneg.
+Qed.
+
+Lemma gcd_1_l : forall n, gcd 1 n == 1.
+Proof.
+ intros. apply gcd_unique; trivial using divide_1_l, le_0_1.
+Qed.
+
+Lemma gcd_1_r : forall n, gcd n 1 == 1.
+Proof.
+ intros. now rewrite gcd_comm, gcd_1_l.
+Qed.
+
+Lemma gcd_diag_nonneg : forall n, 0<=n -> gcd n n == n.
+Proof.
+ intros. apply gcd_unique; trivial using divide_refl.
+Qed.
+
+Lemma gcd_eq_0_l : forall n m, gcd n m == 0 -> n == 0.
+Proof.
+ intros.
+ generalize (gcd_divide_l n m). rewrite H. apply divide_0_l.
+Qed.
+
+Lemma gcd_eq_0_r : forall n m, gcd n m == 0 -> m == 0.
+Proof.
+ intros. apply gcd_eq_0_l with n. now rewrite gcd_comm.
+Qed.
+
+Lemma gcd_eq_0 : forall n m, gcd n m == 0 <-> n == 0 /\ m == 0.
+Proof.
+ intros. split. split.
+ now apply gcd_eq_0_l with m.
+ now apply gcd_eq_0_r with n.
+ intros (EQ,EQ'). rewrite EQ, EQ'. now apply gcd_0_r_nonneg.
+Qed.
+
+Lemma gcd_mul_diag_l : forall n m, 0<=n -> gcd n (n*m) == n.
+Proof.
+ intros n m Hn. apply gcd_unique_alt; trivial.
+ intros q. split. split; trivial. now apply divide_mul_l.
+ now destruct 1.
+Qed.
+
+Lemma divide_gcd_iff : forall n m, 0<=n -> ((n|m) <-> gcd n m == n).
+Proof.
+ intros n m Hn. split. intros (q,Hq). rewrite <- Hq.
+ now apply gcd_mul_diag_l.
+ intros EQ. rewrite <- EQ. apply gcd_divide_r.
+Qed.
+
+End NZGcdProp.
diff --git a/theories/Numbers/Natural/Abstract/NAxioms.v b/theories/Numbers/Natural/Abstract/NAxioms.v
index 2fbfb04c2..ee2a92e84 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 NZSqrt NZLog NZDiv.
+Require Export NZAxioms NZPow NZSqrt NZLog NZDiv NZGcd.
 
 (** From [NZ], we obtain natural numbers just by stating that [pred 0] == 0 *)
 
@@ -39,16 +39,16 @@ Module Type NDivSpecific (Import N : NAxiomsMiniSig')(Import DM : DivMod' N).
  Axiom mod_upper_bound : forall a b, b ~= 0 -> a mod b < b.
 End NDivSpecific.
 
-(** For pow sqrt log2, the NZ axiomatizations are enough. *)
+(** For gcd pow sqrt log2, the NZ axiomatizations are enough. *)
 
 (** We now group everything together. *)
 
 Module Type NAxiomsSig := NAxiomsMiniSig <+ HasCompare <+ Parity
-  <+ NZPow.NZPow <+ NZSqrt.NZSqrt <+ NZLog.NZLog2
+  <+ NZPow.NZPow <+ NZSqrt.NZSqrt <+ NZLog.NZLog2 <+ NZGcd.NZGcd
   <+ DivMod <+ NZDivCommon <+ NDivSpecific.
 
 Module Type NAxiomsSig' := NAxiomsMiniSig' <+ HasCompare <+ Parity
-  <+ NZPow.NZPow' <+ NZSqrt.NZSqrt' <+ NZLog.NZLog2
+  <+ NZPow.NZPow' <+ NZSqrt.NZSqrt' <+ NZLog.NZLog2 <+ NZGcd.NZGcd'
   <+ DivMod' <+ NZDivCommon <+ NDivSpecific.
 
 
diff --git a/theories/Numbers/Natural/Abstract/NGcd.v b/theories/Numbers/Natural/Abstract/NGcd.v
new file mode 100644
index 000000000..5bf33d4d2
--- /dev/null
+++ b/theories/Numbers/Natural/Abstract/NGcd.v
@@ -0,0 +1,214 @@
+(************************************************************************)
+(*  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 the greatest common divisor *)
+
+Require Import NAxioms NSub NZGcd.
+
+Module NGcdProp
+ (Import A : NAxiomsSig')
+ (Import B : NSubProp A).
+
+ Include NZGcdProp A A B.
+
+(** Results concerning divisibility*)
+
+Definition divide_1_r n : (n | 1) -> n == 1
+ := divide_1_r_nonneg n (le_0_l n).
+
+Definition divide_antisym n m : (n | m) -> (m | n) -> n == m
+ := divide_antisym_nonneg n m (le_0_l n) (le_0_l m).
+
+Lemma divide_add_cancel_r : forall n m p, (n | m) -> (n | m + p) -> (n | p).
+Proof.
+ intros n m p (q,Hq) (r,Hr).
+ exists (r-q). rewrite mul_sub_distr_l, Hq, Hr.
+ now rewrite add_comm, add_sub.
+Qed.
+
+Lemma divide_sub_r : forall n m p, (n | m) -> (n | p) -> (n | m - p).
+Proof.
+ intros n m p H H'.
+ destruct (le_ge_cases m p) as [LE|LE].
+ apply sub_0_le in LE. rewrite LE. apply divide_0_r.
+ apply divide_add_cancel_r with p; trivial.
+ now rewrite add_comm, sub_add.
+Qed.
+
+(** Properties of gcd *)
+
+Definition gcd_0_l n : gcd 0 n == n := gcd_0_l_nonneg n (le_0_l n).
+Definition gcd_0_r n : gcd n 0 == n := gcd_0_r_nonneg n (le_0_l n).
+Definition gcd_diag n : gcd n n == n := gcd_diag_nonneg n (le_0_l n).
+Definition gcd_unique' n m p := gcd_unique n m p (le_0_l p).
+Definition gcd_unique_alt' n m p := gcd_unique_alt n m p (le_0_l p).
+Definition divide_gcd_iff' n m := divide_gcd_iff n m (le_0_l n).
+
+Lemma gcd_add_mult_diag_r : forall n m p, gcd n (m+p*n) == gcd n m.
+Proof.
+ intros. apply gcd_unique_alt'.
+ intros. rewrite gcd_divide_iff. split; intros (U,V); split; trivial.
+ apply divide_add_r; trivial. now apply divide_mul_r.
+ apply divide_add_cancel_r with (p*n); trivial.
+ now apply divide_mul_r. now rewrite add_comm.
+Qed.
+
+Lemma gcd_add_diag_r : forall n m, gcd n (m+n) == gcd n m.
+Proof.
+ intros n m. rewrite <- (mul_1_l n) at 2. apply gcd_add_mult_diag_r.
+Qed.
+
+Lemma gcd_sub_diag_r : forall n m, n<=m -> gcd n (m-n) == gcd n m.
+Proof.
+ intros n m H. symmetry.
+ rewrite <- (sub_add n m H) at 1. apply gcd_add_diag_r.
+Qed.
+
+(** On natural numbers, we should use a particular form
+  for the Bezout identity, since we don't have full subtraction. *)
+
+Definition Bezout n m p := exists a, exists b, a*n == p + b*m.
+
+Instance Bezout_wd : Proper (eq==>eq==>eq==>iff) Bezout.
+Proof.
+ unfold Bezout. intros x x' Hx y y' Hy z z' Hz.
+ setoid_rewrite Hx. setoid_rewrite Hy. now setoid_rewrite Hz.
+Qed.
+
+Lemma bezout_1_gcd : forall n m, Bezout n m 1 -> gcd n m == 1.
+Proof.
+ intros n m (q & r & H).
+ apply gcd_unique; trivial using divide_1_l, le_0_1.
+ intros p Hn Hm.
+ apply divide_add_cancel_r with (r*m).
+ now apply divide_mul_r.
+ rewrite add_comm, <- H. now apply divide_mul_r.
+Qed.
+
+(** For strictly positive numbers, we have Bezout in the two directions. *)
+
+Lemma gcd_bezout_pos_pos : forall n, 0<n -> forall m, 0<m ->
+ Bezout n m (gcd n m) /\ Bezout m n (gcd n m).
+Proof.
+ intros n Hn. rewrite <- le_succ_l, <- one_succ in Hn.
+ pattern n. apply strong_right_induction with (z:=1); trivial.
+ unfold Bezout. solve_predicate_wd.
+ clear n Hn. intros n Hn IHn.
+ intros m Hm. rewrite <- le_succ_l, <- one_succ in Hm.
+ pattern m. apply strong_right_induction with (z:=1); trivial.
+ unfold Bezout. solve_predicate_wd.
+ clear m Hm. intros m Hm IHm.
+ destruct (lt_trichotomy n m) as [LT|[EQ|LT]].
+ (* n < m *)
+ destruct (IHm (m-n)) as ((a & b & EQ), (a' & b' & EQ')).
+ rewrite one_succ, le_succ_l.
+ apply lt_add_lt_sub_l; now nzsimpl.
+ apply sub_lt; order'.
+ split.
+ exists (a+b). exists b.
+ rewrite mul_add_distr_r, EQ, mul_sub_distr_l, <- add_assoc.
+ rewrite gcd_sub_diag_r by order.
+ rewrite sub_add. reflexivity. apply mul_le_mono_l; order.
+ exists a'. exists (a'+b').
+ rewrite gcd_sub_diag_r in EQ' by order.
+ rewrite (add_comm a'), mul_add_distr_r, add_assoc, <- EQ'.
+ rewrite mul_sub_distr_l, sub_add. reflexivity. apply mul_le_mono_l; order.
+ (* n = m *)
+ rewrite EQ. rewrite gcd_diag.
+ split.
+ exists 1. exists 0. now nzsimpl.
+ exists 1. exists 0. now nzsimpl.
+ (* m < n *)
+ rewrite gcd_comm, and_comm.
+ apply IHn; trivial.
+ now rewrite <- le_succ_l, <- one_succ.
+Qed.
+
+Lemma gcd_bezout_pos : forall n m, 0<n -> Bezout n m (gcd n m).
+Proof.
+ intros n m Hn.
+ destruct (eq_0_gt_0_cases m) as [EQ|LT].
+ rewrite EQ, gcd_0_r. exists 1. exists 0. now nzsimpl.
+ now apply gcd_bezout_pos_pos.
+Qed.
+
+(** For arbitrary natural numbers, we could only say that at least
+  one of the Bezout identities holds. *)
+
+Lemma gcd_bezout : forall n m,
+ Bezout n m (gcd n m) \/ Bezout m n (gcd n m).
+Proof.
+ intros n m.
+ destruct (eq_0_gt_0_cases n) as [EQ|LT].
+ right. rewrite EQ, gcd_0_l. exists 1. exists 0. now nzsimpl.
+ left. now apply gcd_bezout_pos.
+Qed.
+
+Lemma gcd_mul_mono_l :
+  forall n m p, gcd (p * n) (p * m) == p * gcd n m.
+Proof.
+ intros n m p.
+ apply gcd_unique'.
+ apply mul_divide_mono_l, gcd_divide_l.
+ apply mul_divide_mono_l, gcd_divide_r.
+ intros q H H'.
+ destruct (eq_0_gt_0_cases n) as [EQ|LT].
+ rewrite EQ in *. now rewrite gcd_0_l.
+ destruct (gcd_bezout_pos n m) as (a & b & EQ); trivial.
+ apply divide_add_cancel_r with (p*m*b).
+ now apply divide_mul_l.
+ rewrite <- mul_assoc, <- mul_add_distr_l, add_comm, (mul_comm m), <- EQ.
+ rewrite (mul_comm a), mul_assoc.
+ now apply divide_mul_l.
+Qed.
+
+Lemma gcd_mul_mono_r :
+ forall n m p, gcd (n*p) (m*p) == gcd n m * p.
+Proof.
+ intros. rewrite !(mul_comm _ p). apply gcd_mul_mono_l.
+Qed.
+
+Lemma gauss : forall n m p, (n | m * p) -> gcd n m == 1 -> (n | p).
+Proof.
+ intros n m p H G.
+ destruct (eq_0_gt_0_cases n) as [EQ|LT].
+ rewrite EQ in *. rewrite gcd_0_l in G. now rewrite <- (mul_1_l p), <- G.
+ destruct (gcd_bezout_pos n m) as (a & b & EQ); trivial.
+ rewrite G in EQ.
+ apply divide_add_cancel_r with (m*p*b).
+ now apply divide_mul_l.
+ rewrite (mul_comm _ b), mul_assoc. rewrite <- (mul_1_l p) at 2.
+ rewrite <- mul_add_distr_r, add_comm, <- EQ.
+ now apply divide_mul_l, divide_factor_r.
+Qed.
+
+Lemma divide_mul_split : forall n m p, n ~= 0 -> (n | m * p) ->
+ exists q, exists r, n == q*r /\ (q | m) /\ (r | p).
+Proof.
+ intros n m p Hn H.
+ assert (G := gcd_nonneg n m).
+ apply le_lteq in G; destruct G as [G|G].
+ destruct (gcd_divide_l n m) as (q,Hq).
+ exists (gcd n m). exists q.
+ split. easy.
+ split. apply gcd_divide_r.
+ destruct (gcd_divide_r n m) as (r,Hr).
+ rewrite <- Hr in H. rewrite <- Hq in H at 1.
+ rewrite <- mul_assoc in H. apply mul_divide_cancel_l in H; [|order].
+ apply gauss with r; trivial.
+ apply mul_cancel_l with (gcd n m); [order|].
+ rewrite mul_1_r.
+ rewrite <- gcd_mul_mono_l, Hq, Hr; order.
+ symmetry in G. apply gcd_eq_0 in G. destruct G as (Hn',_); order.
+Qed.
+
+(** TODO : relation between gcd and division and modulo *)
+
+(** TODO : more about rel_prime (i.e. gcd == 1), about prime ... *)
+
+End NGcdProp.
diff --git a/theories/Numbers/Natural/Abstract/NProperties.v b/theories/Numbers/Natural/Abstract/NProperties.v
index 35b6af834..3fc44124f 100644
--- a/theories/Numbers/Natural/Abstract/NProperties.v
+++ b/theories/Numbers/Natural/Abstract/NProperties.v
@@ -7,10 +7,10 @@
 (************************************************************************)
 
 Require Export NAxioms.
-Require Import NMaxMin NParity NPow NSqrt NLog NDiv.
+Require Import NMaxMin NParity NPow NSqrt NLog NDiv NGcd.
 
 (** This functor summarizes all known facts about N. *)
 
 Module Type NProp (N:NAxiomsSig) :=
  NMaxMinProp N <+ NParityProp N <+ NPowProp N <+ NSqrtProp N
-  <+ NLog2Prop N <+ NDivProp N.
+  <+ NLog2Prop N <+ NDivProp N <+ NGcdProp N.
diff --git a/theories/Numbers/Natural/Binary/NBinary.v b/theories/Numbers/Natural/Binary/NBinary.v
index cc57171b3..73ed4e814 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 Nsqrt_def.
+Require Import BinPos Ndiv_def Nsqrt_def Ngcd_def.
 Require Export BinNat.
 Require Import NAxioms NProperties.
 
@@ -169,6 +169,14 @@ Definition sqrt_spec n (H:0<=n) := Nsqrt_spec n.
 Lemma sqrt_neg : forall a, a<0 -> Nsqrt a = 0.
 Proof. destruct a; discriminate. Qed.
 
+(** Gcd *)
+
+Definition gcd_divide_l := Ngcd_divide_l.
+Definition gcd_divide_r := Ngcd_divide_r.
+Definition gcd_greatest := Ngcd_greatest.
+Lemma gcd_nonneg : forall a b, 0 <= Ngcd a b.
+Proof. intros. now destruct (Ngcd a b). Qed.
+
 (** The instantiation of operations.
     Placing them at the very end avoids having indirections in above lemmas. *)
 
@@ -196,6 +204,8 @@ Definition even := Neven.
 Definition odd := Nodd.
 Definition sqrt := Nsqrt.
 Definition log2 := Nlog2.
+Definition divide := Ndivide.
+Definition gcd := Ngcd.
 
 Include NProp
  <+ UsualMinMaxLogicalProperties <+ UsualMinMaxDecProperties.
diff --git a/theories/Numbers/Natural/Peano/NPeano.v b/theories/Numbers/Natural/Peano/NPeano.v
index a828ec821..ce8e03a5a 100644
--- a/theories/Numbers/Natural/Peano/NPeano.v
+++ b/theories/Numbers/Natural/Peano/NPeano.v
@@ -9,7 +9,8 @@
 (************************************************************************)
 
 Require Import
- Bool Peano Peano_dec Compare_dec Plus Minus Le Lt EqNat NAxioms NProperties.
+ Bool Peano Peano_dec Compare_dec Plus Mult Minus Le Lt EqNat
+ NAxioms NProperties.
 
 (** Functions not already defined *)
 
@@ -274,6 +275,63 @@ Proof.
  inversion 1; now subst.
 Qed.
 
+(** * Gcd *)
+
+(** We use Euclid algorithm, which is normally not structural,
+    but Coq is now clever enough to accept this (behind modulo
+    there is a subtraction, which now preserve being a subterm)
+*)
+
+Fixpoint gcd a b :=
+  match a with
+   | O => b
+   | S a' => gcd (b mod (S a')) (S a')
+  end.
+
+Definition divide x y := exists z, x*z=y.
+Notation "( x | y )" := (divide x y) (at level 0) : nat_scope.
+
+Lemma gcd_divide : forall a b, (gcd a b | a) /\ (gcd a b | b).
+Proof.
+ fix 1.
+ intros [|a] b; simpl.
+ split.
+  exists 0; now rewrite <- mult_n_O.
+  exists 1; now rewrite <- mult_n_Sm, <- mult_n_O.
+ fold (b mod (S a)).
+ destruct (gcd_divide (b mod (S a)) (S a)) as (H,H').
+ set (a':=S a) in *.
+ split; auto.
+ rewrite (div_mod b a') at 2 by discriminate.
+ destruct H as (u,Hu), H' as (v,Hv).
+ exists ((b/a')*v + u).
+ rewrite mult_plus_distr_l.
+ now rewrite (mult_comm _ v), mult_assoc, Hv, Hu.
+Qed.
+
+Lemma gcd_divide_l : forall a b, (gcd a b | a).
+Proof.
+ intros. apply gcd_divide.
+Qed.
+
+Lemma gcd_divide_r : forall a b, (gcd a b | b).
+Proof.
+ intros. apply gcd_divide.
+Qed.
+
+Lemma gcd_greatest : forall a b c, (c|a) -> (c|b) -> (c|gcd a b).
+Proof.
+ fix 1.
+ intros [|a] b; simpl; auto.
+ fold (b mod (S a)).
+ intros c H H'. apply gcd_greatest; auto.
+ set (a':=S a) in *.
+ rewrite (div_mod b a') in H' by discriminate.
+ destruct H as (u,Hu), H' as (v,Hv).
+ exists (v - u * (b/a')).
+ now rewrite mult_minus_distr_l, mult_assoc, Hu, Hv, minus_plus.
+Qed.
+
 
 (** * Implementation of [NAxiomsSig] by [nat] *)
 
@@ -459,6 +517,14 @@ Program Instance mod_wd : Proper (eq==>eq==>eq) modulo.
 Definition div_mod := div_mod.
 Definition mod_upper_bound := mod_upper_bound.
 
+Definition divide := divide.
+Definition gcd := gcd.
+Definition gcd_divide_l := gcd_divide_l.
+Definition gcd_divide_r := gcd_divide_r.
+Definition gcd_greatest := gcd_greatest.
+Lemma gcd_nonneg : forall a b, 0<=gcd a b.
+Proof. intros. apply le_O_n. Qed.
+
 (** Generic Properties *)
 
 Include NProp
diff --git a/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v b/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
index 547f51b3d..7b90aa09e 100644
--- a/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
+++ b/theories/Numbers/Natural/SpecViaZ/NSigNAxioms.v
@@ -299,6 +299,42 @@ destruct (Z_mod_lt [a] [b]); auto.
 generalize (spec_pos b); auto with zarith.
 Qed.
 
+(** Gcd *)
+
+Definition divide n m := exists p, n*p == m.
+Local Notation "( x | y )" := (divide x y) (at level 0).
+
+Lemma spec_divide : forall n m, (n|m) <-> Zdivide' [n] [m].
+Proof.
+ intros n m. split.
+ intros (p,H). exists [p]. revert H; now zify.
+ intros (z,H). exists (of_N (Zabs_N z)). zify.
+ rewrite spec_of_N, Z_of_N_abs.
+ rewrite <- (Zabs_eq [n]) by apply spec_pos.
+ rewrite <- Zabs_Zmult, H.
+ apply Zabs_eq, spec_pos.
+Qed.
+
+Lemma gcd_divide_l : forall n m, (gcd n m | n).
+Proof.
+ intros n m. apply spec_divide. zify. apply Zgcd_divide_l.
+Qed.
+
+Lemma gcd_divide_r : forall n m, (gcd n m | m).
+Proof.
+ intros n m. apply spec_divide. zify. apply Zgcd_divide_r.
+Qed.
+
+Lemma gcd_greatest : forall n m p, (p|n) -> (p|m) -> (p|gcd n m).
+Proof.
+ intros n m p. rewrite !spec_divide. zify. apply Zgcd_greatest.
+Qed.
+
+Lemma gcd_nonneg : forall n m, 0 <= gcd n m.
+Proof.
+ intros. zify. apply Zgcd_nonneg.
+Qed.
+
 (** Recursion *)
 
 Definition recursion (A : Type) (a : A) (f : N.t -> A -> A) (n : N.t) :=
diff --git a/theories/Numbers/vo.itarget b/theories/Numbers/vo.itarget
index fbfaa8007..e380d835e 100644
--- a/theories/Numbers/vo.itarget
+++ b/theories/Numbers/vo.itarget
@@ -30,6 +30,7 @@ Integer/Abstract/ZDivEucl.vo
 Integer/Abstract/ZMaxMin.vo
 Integer/Abstract/ZParity.vo
 Integer/Abstract/ZPow.vo
+Integer/Abstract/ZGcd.vo
 Integer/BigZ/BigZ.vo
 Integer/BigZ/ZMake.vo
 Integer/Binary/ZBinary.vo
@@ -50,6 +51,7 @@ NatInt/NZDiv.vo
 NatInt/NZPow.vo
 NatInt/NZSqrt.vo
 NatInt/NZLog.vo
+NatInt/NZGcd.vo
 Natural/Abstract/NAddOrder.vo
 Natural/Abstract/NAdd.vo
 Natural/Abstract/NAxioms.vo
@@ -67,6 +69,7 @@ Natural/Abstract/NParity.vo
 Natural/Abstract/NPow.vo
 Natural/Abstract/NSqrt.vo
 Natural/Abstract/NLog.vo
+Natural/Abstract/NGcd.vo
 Natural/BigN/BigN.vo
 Natural/BigN/Nbasic.vo
 Natural/BigN/NMake_gen.vo
diff --git a/theories/PArith/BinPos.v b/theories/PArith/BinPos.v
index 9c8775f82..c42f9cceb 100644
--- a/theories/PArith/BinPos.v
+++ b/theories/PArith/BinPos.v
@@ -1047,11 +1047,18 @@ Proof.
  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.
+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.
+ assert (IMPL : forall p q r, q<r -> p*q < p*r).
+  induction p using Prect; auto.
+  intros q r H. rewrite 2 Pmult_Sn_m.
+  apply Pplus_lt_mono; auto.
+ split; auto.
+ intros H.
+ destruct (Pcompare_spec q r) as [EQ|LT|GT]; trivial.
+ rewrite EQ in H. elim (Plt_irrefl _ H).
+ apply (IMPL p) in GT.
+ elim (Plt_irrefl (p*q)). eapply Plt_trans; eauto.
 Qed.
 
 Lemma Pmult_lt_mono : forall p q r s, p<q -> r<s -> p*r < q*s.
@@ -1216,6 +1223,43 @@ Proof.
   unfold Pminus; rewrite U; simpl; auto.
 Qed.
 
+Lemma Pplus_minus_eq : forall p q, p+q-q = p.
+Proof.
+ intros. apply Pplus_reg_l with q.
+ rewrite Pplus_minus.
+ apply Pplus_comm.
+ now rewrite ZC4, Pplus_comm, Plt_plus_r.
+Qed.
+
+Lemma Pmult_minus_distr_l : forall p q r, r<q -> p*(q-r) = p*q-p*r.
+Proof.
+ intros p q r H.
+ apply Pplus_reg_l with (p*r).
+ rewrite <- Pmult_plus_distr_l.
+ rewrite Pplus_minus.
+ symmetry. apply Pplus_minus.
+ apply (Pmult_lt_mono_l p) in H.
+ now rewrite ZC4, H.
+ now rewrite ZC4, H.
+Qed.
+
+Lemma Pminus_decr : forall n m, m<n -> n-m < n.
+Proof.
+ intros n m LT.
+ apply Pplus_lt_mono_l with m.
+ rewrite Pplus_minus.
+ rewrite Pplus_comm.
+ apply Plt_plus_r.
+ now rewrite ZC4, LT.
+Qed.
+
+Lemma Pminus_xI_xI : forall n m, m<n -> n~1 - m~1 = (n-m)~0.
+Proof.
+ intros. unfold Pminus. simpl.
+ destruct (Pminus_mask_Gt n m) as (p & -> & _); trivial.
+ now rewrite ZC4, H.
+Qed.
+
 (** When x<y, the substraction of x by y returns 1 *)
 
 Lemma Pminus_mask_Lt : forall p q:positive, p<q -> Pminus_mask p q = IsNeg.
@@ -1249,7 +1293,7 @@ Fixpoint Psize (p:positive) : nat :=
     | p~0 => S (Psize p)
   end.
 
-Lemma Psize_monotone : forall p q, (p?=q) Eq = Lt -> (Psize p <= Psize q)%nat.
+Lemma Psize_monotone : forall p q, p<q -> (Psize p <= Psize q)%nat.
 Proof.
   assert (le0 : forall n, (0<=n)%nat) by (induction n; auto).
   assert (leS : forall n m, (n<=m -> S n <= S m)%nat) by (induction 1; auto).
diff --git a/theories/PArith/PArith.v b/theories/PArith/PArith.v
index f453af386..8688c5013 100644
--- a/theories/PArith/PArith.v
+++ b/theories/PArith/PArith.v
@@ -8,8 +8,4 @@
 
 (** Library for positive natural numbers *)
 
-Require Export BinPos.
-Require Export Pnat.
-Require Export Pminmax.
-Require Export Psqrt.
-Require Export POrderedType.
+Require Export BinPos Pnat Pminmax Psqrt Pgcd POrderedType.
diff --git a/theories/PArith/Pgcd.v b/theories/PArith/Pgcd.v
new file mode 100644
index 000000000..22d25dd8c
--- /dev/null
+++ b/theories/PArith/Pgcd.v
@@ -0,0 +1,265 @@
+(************************************************************************)
+(*  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 Le Plus.
+
+Local Open Scope positive_scope.
+
+(** * Divisibility *)
+
+Definition Pdivide p q := exists r, p*r = q.
+Notation "( p | q )" := (Pdivide p q) (at level 0) : positive_scope.
+
+Lemma Pdivide_add_cancel_r : forall p q r, (p | q) -> (p | q + r) -> (p | r).
+Proof.
+ intros p q r (s,Hs) (t,Ht).
+ exists (t-s).
+ rewrite Pmult_minus_distr_l.
+ rewrite Hs, Ht.
+ rewrite Pplus_comm. apply Pplus_minus_eq.
+ apply Pmult_lt_mono_l with p.
+ rewrite Hs, Ht.
+ apply Plt_plus_r.
+Qed.
+
+Lemma Pdivide_xO_xI : forall p q r, (p | q~0) -> (p | r~1) -> (p | q).
+Proof.
+ intros p q r (s,Hs) (t,Ht).
+ destruct p.
+  destruct s; try discriminate.
+  rewrite Pmult_xO_permute_r in Hs. exists s; congruence.
+ simpl in Ht. discriminate.
+ exists q; auto.
+Qed.
+
+Lemma Pdivide_xO_xO : forall p q, (p~0|q~0) <-> (p|q).
+Proof.
+ intros p q. split; intros (r,H); simpl in *.
+ injection H. exists r; auto.
+ exists r; simpl; f_equal; auto.
+Qed.
+
+Lemma Pdivide_mult_l : forall p q r, (p|q) -> (p|q*r).
+Proof.
+ intros p q r (s,H). exists (s*r). rewrite Pmult_assoc; now f_equal.
+Qed.
+
+Lemma Pdivide_mult_r : forall p q r, (p|r) -> (p|q*r).
+Proof.
+ intros p q r. rewrite Pmult_comm. apply Pdivide_mult_l.
+Qed.
+
+
+(** * Definition of a Pgcd function for positive numbers *)
+
+(** Instead of the Euclid algorithm, we use here the Stein binary
+   algorithm, which is faster for this representation. This algorithm
+   is almost structural, but in the last cases we do some recursive
+   calls on subtraction, hence the need for a counter.
+*)
+
+Fixpoint Pgcdn (n: nat) (a b : positive) : positive :=
+  match n with
+    | O => 1
+    | S n =>
+      match a,b with
+	| 1, _ => 1
+	| _, 1 => 1
+	| a~0, b~0 => (Pgcdn n a b)~0
+	| _  , b~0 => Pgcdn n a b
+	| a~0, _   => Pgcdn n a b
+	| a'~1, b'~1 =>
+          match (a' ?= b') Eq with
+	    | Eq => a
+	    | Lt => Pgcdn n (b'-a') a
+	    | Gt => Pgcdn n (a'-b') b
+          end
+      end
+  end.
+
+(** We'll show later that we need at most (log2(a.b)) loops *)
+
+Definition Pgcd (a b: positive) := Pgcdn (Psize a + Psize b)%nat a b.
+
+
+(** * Generalized Gcd, also computing the division of a and b by the gcd *)
+
+Fixpoint Pggcdn (n: nat) (a b : positive) : (positive*(positive*positive)) :=
+  match n with
+    | O => (1,(a,b))
+    | S n =>
+      match a,b with
+	| 1, _ => (1,(1,b))
+	| _, 1 => (1,(a,1))
+	| a~0, b~0 =>
+           let (g,p) := Pggcdn n a b in
+           (g~0,p)
+	| _, b~0 =>
+           let '(g,(aa,bb)) := Pggcdn n a b in
+           (g,(aa, bb~0))
+	| a~0, _ =>
+           let '(g,(aa,bb)) := Pggcdn n a b in
+           (g,(aa~0, bb))
+	| a'~1, b'~1 =>
+           match Pcompare a' b' Eq with
+	     | Eq => (a,(1,1))
+	     | Lt =>
+	        let '(g,(ba,aa)) := Pggcdn n (b'-a') a in
+	        (g,(aa, aa + ba~0))
+	     | Gt =>
+		let '(g,(ab,bb)) := Pggcdn n (a'-b') b in
+		(g,(bb + ab~0, bb))
+	   end
+      end
+  end.
+
+Definition Pggcd (a b: positive) := Pggcdn (Psize a + Psize b)%nat a b.
+
+(** The first component of Pggcd is Pgcd *)
+
+Lemma Pggcdn_gcdn : forall n a b,
+  fst (Pggcdn n a b) = Pgcdn n a b.
+Proof.
+ induction n.
+ simpl; auto.
+ destruct a, b; simpl; auto; try case Pcompare_spec; simpl; trivial;
+  rewrite <- IHn; destruct Pggcdn as (g,(u,v)); simpl; auto.
+Qed.
+
+Lemma Pggcd_gcd : forall a b, fst (Pggcd a b) = Pgcd a b.
+Proof.
+ unfold Pggcd, Pgcd. intros. apply Pggcdn_gcdn.
+Qed.
+
+(** The other components of Pggcd are indeed the correct factors. *)
+
+Ltac destr_pggcdn IHn :=
+ match goal with |- context [ Pggcdn _ ?x ?y ] =>
+  generalize (IHn x y); destruct Pggcdn as (g,(u,v)); simpl
+ end.
+
+Lemma Pggcdn_correct_divisors : forall n a b,
+  let '(g,(aa,bb)) := Pggcdn n a b in
+  a = g*aa /\ b = g*bb.
+Proof.
+ induction n.
+ simpl; auto.
+ destruct a, b; simpl; auto; try case Pcompare_spec; try destr_pggcdn IHn.
+ (* Eq *)
+ intros ->. now rewrite Pmult_comm.
+ (* Lt *)
+ intros (H',H) LT; split; auto.
+ rewrite Pmult_plus_distr_l, Pmult_xO_permute_r, <- H, <- H'.
+ simpl. f_equal. symmetry.
+ apply Pplus_minus; auto. rewrite ZC4; rewrite LT; auto.
+ (* Gt *)
+ intros (H',H) LT; split; auto.
+ rewrite Pmult_plus_distr_l, Pmult_xO_permute_r, <- H, <- H'.
+ simpl. f_equal. symmetry.
+ apply Pplus_minus; auto. rewrite ZC4; rewrite LT; auto.
+ (* Then... *)
+ intros (H,H'); split; auto. rewrite Pmult_xO_permute_r, H'; auto.
+ intros (H,H'); split; auto. rewrite Pmult_xO_permute_r, H; auto.
+ intros (H,H'); split; subst; auto.
+Qed.
+
+Lemma Pggcd_correct_divisors : forall a b,
+  let '(g,(aa,bb)) := Pggcd a b in
+  a=g*aa /\ b=g*bb.
+Proof.
+ unfold Pggcd. intros. apply Pggcdn_correct_divisors.
+Qed.
+
+(** We can use this fact to prove a part of the Pgcd correctness *)
+
+Lemma Pgcd_divide_l : forall a b, (Pgcd a b | a).
+Proof.
+ intros a b. rewrite <- Pggcd_gcd. generalize (Pggcd_correct_divisors a b).
+ destruct Pggcd as (g,(aa,bb)); simpl. intros (H,_). exists aa; auto.
+Qed.
+
+Lemma Pgcd_divide_r : forall a b, (Pgcd a b | b).
+Proof.
+ intros a b. rewrite <- Pggcd_gcd. generalize (Pggcd_correct_divisors a b).
+ destruct Pggcd as (g,(aa,bb)); simpl. intros (_,H). exists bb; auto.
+Qed.
+
+(** We now prove directly that gcd is the greatest amongst common divisors *)
+
+Lemma Pgcdn_greatest : forall n a b, (Psize a + Psize b <= n)%nat ->
+ forall p, (p|a) -> (p|b) -> (p|Pgcdn n a b).
+Proof.
+ induction n.
+ destruct a, b; simpl; inversion 1.
+ destruct a, b; simpl; try case Pcompare_spec; simpl; auto.
+ (* Lt *)
+ intros LT LE p Hp1 Hp2. apply IHn; clear IHn; trivial.
+ apply le_S_n in LE. eapply le_trans; [|eapply LE].
+ rewrite plus_comm, <- plus_n_Sm, <- plus_Sn_m.
+ apply plus_le_compat; trivial.
+ apply Psize_monotone, Pminus_decr, LT.
+ apply Pdivide_xO_xI with a; trivial.
+ apply (Pdivide_add_cancel_r p a~1); trivial.
+ rewrite <- Pminus_xI_xI, Pplus_minus; trivial.
+ simpl. now rewrite ZC4, LT.
+ (* Gt *)
+ intros LT LE p Hp1 Hp2. apply IHn; clear IHn; trivial.
+ apply le_S_n in LE. eapply le_trans; [|eapply LE].
+ apply plus_le_compat; trivial.
+ apply Psize_monotone, Pminus_decr, LT.
+ apply Pdivide_xO_xI with b; trivial.
+ apply (Pdivide_add_cancel_r p b~1); trivial.
+ rewrite <- Pminus_xI_xI, Pplus_minus; trivial.
+ simpl. now rewrite ZC4, LT.
+ (* a~1 b~0 *)
+ intros LE p Hp1 Hp2. apply IHn; clear IHn; trivial.
+ apply le_S_n in LE. simpl. now rewrite plus_n_Sm.
+ apply Pdivide_xO_xI with a; trivial.
+ (* a~0 b~1 *)
+ intros LE p Hp1 Hp2. apply IHn; clear IHn; trivial.
+ simpl. now apply le_S_n.
+ apply Pdivide_xO_xI with b; trivial.
+ (* a~0 b~0 *)
+ intros LE p Hp1 Hp2.
+ destruct p.
+ change (Pgcdn n a b)~0 with (2*(Pgcdn n a b)).
+ apply Pdivide_mult_r.
+ apply IHn; clear IHn.
+ apply le_S_n in LE. apply le_Sn_le. now rewrite plus_n_Sm.
+ apply Pdivide_xO_xI with p; trivial. exists 1; now rewrite Pmult_1_r.
+ apply Pdivide_xO_xI with p; trivial. exists 1; now rewrite Pmult_1_r.
+ apply Pdivide_xO_xO.
+ apply IHn; clear IHn.
+ apply le_S_n in LE. apply le_Sn_le. now rewrite plus_n_Sm.
+ now apply Pdivide_xO_xO.
+ now apply Pdivide_xO_xO.
+ exists (Pgcdn n a b)~0. auto.
+Qed.
+
+Lemma Pgcd_greatest : forall a b p, (p|a) -> (p|b) -> (p|Pgcd a b).
+Proof.
+ intros. apply Pgcdn_greatest; auto.
+Qed.
+
+(** As a consequence, the rests after division by gcd are relatively prime *)
+
+Lemma Pggcd_greatest : forall a b,
+ let (aa,bb) := snd (Pggcd a b) in
+ forall p, (p|aa) -> (p|bb) -> p=1.
+Proof.
+ intros. generalize (Pgcd_greatest a b) (Pggcd_correct_divisors a b).
+ rewrite <- Pggcd_gcd. destruct Pggcd as (g,(aa,bb)); simpl.
+ intros H (EQa,EQb) p Hp1 Hp2; subst.
+ assert (H' : (g*p | g)).
+  apply H.
+  destruct Hp1 as (r,Hr). exists r. now rewrite <- Hr, Pmult_assoc.
+  destruct Hp2 as (r,Hr). exists r. now rewrite <- Hr, Pmult_assoc.
+ destruct H' as (q,H').
+ apply Pmult_1_inversion_l with q.
+ apply Pmult_reg_l with g. now rewrite Pmult_assoc, Pmult_1_r.
+Qed.
diff --git a/theories/PArith/vo.itarget b/theories/PArith/vo.itarget
index e9139da8a..fee542494 100644
--- a/theories/PArith/vo.itarget
+++ b/theories/PArith/vo.itarget
@@ -3,4 +3,5 @@ Pnat.vo
 POrderedType.vo
 Pminmax.vo
 Psqrt.vo
+Pgcd.vo
 PArith.vo
\ No newline at end of file
diff --git a/theories/QArith/Qreduction.v b/theories/QArith/Qreduction.v
index ad16f8a25..4888ce4b1 100644
--- a/theories/QArith/Qreduction.v
+++ b/theories/QArith/Qreduction.v
@@ -9,7 +9,7 @@
 (** Normalisation functions for rational numbers. *)
 
 Require Export QArith_base.
-Require Import Znumtheory.
+Require Import Zgcd_def Znumtheory.
 
 (** First, a function that (tries to) build a positive back from a Z. *)
 
@@ -41,23 +41,16 @@ Definition Qred (q:Q) :=
 Lemma Qred_correct : forall q, (Qred q) == q.
 Proof.
   unfold Qred, Qeq; intros (n,d); simpl.
-  generalize (Zggcd_gcd n ('d)) (Zgcd_is_pos n ('d))
-    (Zgcd_is_gcd n ('d)) (Zggcd_correct_divisors n ('d)).
+  generalize (Zggcd_gcd n ('d)) (Zgcd_nonneg n ('d))
+    (Zggcd_correct_divisors n ('d)).
   destruct (Zggcd n (Zpos d)) as (g,(nn,dd)); simpl.
   Open Scope Z_scope.
-  intuition.
-  rewrite <- H in H0,H1; clear H.
-  rewrite H3; rewrite H4.
-  assert (0 <> g).
-  intro; subst g; discriminate.
-
-  assert (0 < dd).
-  apply Zmult_gt_0_lt_0_reg_r with g.
-  omega.
-  rewrite Zmult_comm.
-  rewrite <- H4; compute; auto.
-  rewrite Z2P_correct; auto.
-  ring.
+  intros Hg LE (Hn,Hd). rewrite Hd, Hn.
+  rewrite <- Hg in LE; clear Hg.
+  assert (0 <> g) by (intro; subst; discriminate).
+  rewrite Z2P_correct. ring.
+  apply Zmult_gt_0_lt_0_reg_r with g; auto with zarith.
+  now rewrite Zmult_comm, <- Hd.
   Close Scope Z_scope.
 Qed.
 
@@ -67,10 +60,10 @@ Proof.
   unfold Qred, Qeq in *; simpl in *.
   Open Scope Z_scope.
   generalize (Zggcd_gcd a ('b)) (Zgcd_is_gcd a ('b))
-    (Zgcd_is_pos a ('b)) (Zggcd_correct_divisors a ('b)).
+    (Zgcd_nonneg a ('b)) (Zggcd_correct_divisors a ('b)).
   destruct (Zggcd a (Zpos b)) as (g,(aa,bb)).
   generalize (Zggcd_gcd c ('d)) (Zgcd_is_gcd c ('d))
-    (Zgcd_is_pos c ('d)) (Zggcd_correct_divisors c ('d)).
+    (Zgcd_nonneg c ('d)) (Zggcd_correct_divisors c ('d)).
   destruct (Zggcd c (Zpos d)) as (g',(cc,dd)).
   simpl.
   intro H; rewrite <- H; clear H.
@@ -78,10 +71,9 @@ Proof.
   intro H; rewrite <- H; clear H.
   intros Hg1 Hg2 (Hg3,Hg4).
   intros.
-  assert (g <> 0).
-  intro; subst g; discriminate.
-  assert (g' <> 0).
-  intro; subst g'; discriminate.
+  assert (g <> 0) by (intro; subst g; discriminate).
+  assert (g' <> 0) by (intro; subst g'; discriminate).
+
   elim (rel_prime_cross_prod aa bb cc dd).
   congruence.
   unfold rel_prime in |- *.
diff --git a/theories/ZArith/ZArith.v b/theories/ZArith/ZArith.v
index 8199ed369..3956e9160 100644
--- a/theories/ZArith/ZArith.v
+++ b/theories/ZArith/ZArith.v
@@ -12,7 +12,7 @@ Require Export ZArith_base.
 
 (** Extra definitions *)
 
-Require Export Zpow_def Zsqrt_def Zlog_def.
+Require Export Zpow_def Zsqrt_def Zlog_def Zgcd_def.
 
 (** Extra modules using [Omega] or [Ring]. *)
 
diff --git a/theories/ZArith/Zgcd_alt.v b/theories/ZArith/Zgcd_alt.v
index 1e0258fcd..bcb9e42ae 100644
--- a/theories/ZArith/Zgcd_alt.v
+++ b/theories/ZArith/Zgcd_alt.v
@@ -6,7 +6,7 @@
 (*         *       GNU Lesser General Public License Version 2.1        *)
 (************************************************************************)
 
-(** * Zgcd_alt : an alternate version of Zgcd, based on Euler's algorithm *)
+(** * Zgcd_alt : an alternate version of Zgcd, based on Euclid's algorithm *)
 
 (**
 Author: Pierre Letouzey
@@ -15,13 +15,14 @@ Author: Pierre Letouzey
 (** The alternate [Zgcd_alt] given here used to be the main [Zgcd]
     function (see file [Znumtheory]), but this main [Zgcd] is now
     based on a modern binary-efficient algorithm. This earlier
-    version, based on Euler's algorithm of iterated modulo, is kept
+    version, based on Euclid's algorithm of iterated modulo, is kept
     here due to both its intrinsic interest and its use as reference
     point when proving gcd on Int31 numbers *)
 
 Require Import ZArith_base.
 Require Import ZArithRing.
 Require Import Zdiv.
+Require Import Zgcd_def.
 Require Import Znumtheory.
 
 Open Scope Z_scope.
diff --git a/theories/ZArith/Zgcd_def.v b/theories/ZArith/Zgcd_def.v
new file mode 100644
index 000000000..1785d4e66
--- /dev/null
+++ b/theories/ZArith/Zgcd_def.v
@@ -0,0 +1,139 @@
+(************************************************************************)
+(*  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 BinInt Pgcd.
+
+Open Local Scope Z_scope.
+
+(** * Definition of a Pgcd function for relative numbers *)
+
+Definition Zgcd (a b : Z) : Z :=
+  match a,b with
+    | Z0, _ => Zabs b
+    | _, Z0 => Zabs a
+    | Zpos a, Zpos b => Zpos (Pgcd a b)
+    | Zpos a, Zneg b => Zpos (Pgcd a b)
+    | Zneg a, Zpos b => Zpos (Pgcd a b)
+    | Zneg a, Zneg b => Zpos (Pgcd a b)
+  end.
+
+(** * Generalized Gcd, also computing division of a and b by gcd. *)
+
+Definition Zggcd (a b : Z) : Z*(Z*Z) :=
+  match a,b with
+    | Z0, _ => (Zabs b,(Z0, Zsgn b))
+    | _, Z0 => (Zabs a,(Zsgn a, Z0))
+    | Zpos a, Zpos b =>
+       let '(g,(aa,bb)) := Pggcd a b in (Zpos g, (Zpos aa, Zpos bb))
+    | Zpos a, Zneg b =>
+       let '(g,(aa,bb)) := Pggcd a b in (Zpos g, (Zpos aa, Zneg bb))
+    | Zneg a, Zpos b =>
+       let '(g,(aa,bb)) := Pggcd a b in (Zpos g, (Zneg aa, Zpos bb))
+    | Zneg a, Zneg b =>
+       let '(g,(aa,bb)) := Pggcd a b in (Zpos g, (Zneg aa, Zneg bb))
+  end.
+
+(** The first component of Zggcd is Zgcd *)
+
+Lemma Zggcd_gcd : forall a b, fst (Zggcd a b) = Zgcd a b.
+Proof.
+ intros [ |p|p] [ |q|q]; simpl; auto;
+  generalize (Pggcd_gcd p q); destruct Pggcd as (g,(aa,bb)); simpl; congruence.
+Qed.
+
+(** The other components of Zggcd are indeed the correct factors. *)
+
+Lemma Zggcd_correct_divisors : forall a b,
+  let '(g,(aa,bb)) := Zggcd a b in
+  a = g*aa /\ b = g*bb.
+Proof.
+ intros [ |p|p] [ |q|q]; simpl; rewrite ?Pmult_1_r; auto;
+  generalize (Pggcd_correct_divisors p q);
+  destruct Pggcd as (g,(aa,bb)); simpl; destruct 1; subst; auto.
+Qed.
+
+(** Divisibility. We use here a simple "exists", while the historical
+    Znumtheory.Zdivide is a specialized inductive type. *)
+
+Definition Zdivide' x y := exists z, x*z = y.
+
+Local Notation "( x | y )" := (Zdivide' x y) (at level 0).
+
+Lemma Zgcd_divide_l : forall a b, (Zgcd a b | a).
+Proof.
+ intros a b. rewrite <- Zggcd_gcd. generalize (Zggcd_correct_divisors a b).
+ destruct Zggcd as (g,(aa,bb)); simpl. intros (H,_). exists aa; auto.
+Qed.
+
+Lemma Zgcd_divide_r : forall a b, (Zgcd a b | b).
+Proof.
+ intros a b. rewrite <- Zggcd_gcd. generalize (Zggcd_correct_divisors a b).
+ destruct Zggcd as (g,(aa,bb)); simpl. intros (_,H). exists bb; auto.
+Qed.
+
+Lemma Zdivide'_0_l : forall x, (0|x) -> x = 0.
+Proof.
+ intros x (y,H). auto.
+Qed.
+
+Lemma Zdivide'_Zpos_Zneg_r : forall x y, (x|Zpos y) <-> (x|Zneg y).
+Proof.
+ intros x y; split; intros (z,H); exists (-z);
+  now rewrite <- Zopp_mult_distr_r, H.
+Qed.
+
+Lemma Zdivide'_Zpos_Zneg_l : forall x y, (Zpos x|y) <-> (Zneg x|y).
+Proof.
+ intros x y; split; intros (z,H); exists (-z);
+  now rewrite <- Zmult_opp_comm.
+Qed.
+
+Lemma Zdivide'_Pdivide : forall p q, (Zpos p|Zpos q) <-> (p|q)%positive.
+Proof.
+ intros p q. split.
+ intros ([ |r|r],H); simpl in *; discriminate H || injection H. exists r; auto.
+ intros (r,H). exists (Zpos r); simpl; f_equal; auto.
+Qed.
+
+Lemma Zgcd_greatest : forall a b c, (c|a) -> (c|b) -> (c|Zgcd a b).
+Proof.
+ assert (D : forall p q c, (c|Zpos p) -> (c|Zpos q) -> (c|Zpos (Pgcd p q))).
+  intros p q [|r|r] H H'.
+  apply Zdivide'_0_l in H. discriminate.
+  apply Zdivide'_Pdivide, Pgcd_greatest; now apply Zdivide'_Pdivide.
+  apply Zdivide'_Zpos_Zneg_l, Zdivide'_Pdivide, Pgcd_greatest;
+   now apply Zdivide'_Pdivide, Zdivide'_Zpos_Zneg_l.
+ intros [ |p|p] [ |q|q]; simpl; auto; intros; try apply D; trivial;
+  now apply Zdivide'_Zpos_Zneg_r.
+Qed.
+
+Lemma Zgcd_nonneg : forall a b, 0 <= Zgcd a b.
+Proof.
+ intros [ |p|p] [ |q|q]; discriminate.
+Qed.
+
+(*
+(** Zggcd produces coefficients that are relatively prime *)
+
+Lemma Zggcd_greatest : forall a b,
+ let (aa,bb) := snd (Zggcd a b) in
+ forall z, (z|aa) -> (z|bb) -> z=1 \/ z=-1.
+Proof.
+ intros [ |p|p] [ |q|q]; simpl.
+*)
+
+
+(** Zggcd and opp : an auxiliary result used in QArith *)
+
+Theorem Zggcd_opp: forall x y,
+  Zggcd (-x) y = let '(g,(aa,bb)) := Zggcd x y in
+                 (g,(-aa,bb)).
+Proof.
+intros [|x|x] [|y|y]; unfold Zggcd, Zopp; auto;
+ destruct (Pggcd x y) as (g,(aa,bb)); auto.
+Qed.
diff --git a/theories/ZArith/Znumtheory.v b/theories/ZArith/Znumtheory.v
index 7f66bd0a6..24fdb3143 100644
--- a/theories/ZArith/Znumtheory.v
+++ b/theories/ZArith/Znumtheory.v
@@ -10,8 +10,12 @@ Require Import ZArith_base.
 Require Import ZArithRing.
 Require Import Zcomplements.
 Require Import Zdiv.
+Require Import Zgcd_def.
 Require Import Wf_nat.
-Open Local Scope Z_scope.
+
+(** For compatibility reasons, this Open Scope isn't local as it should *)
+
+Open Scope Z_scope.
 
 (** This file contains some notions of number theory upon Z numbers:
      - a divisibility predicate [Zdivide]
@@ -19,7 +23,7 @@ Open Local Scope Z_scope.
      - Euclid algorithm [euclid]
      - a relatively prime predicate [rel_prime]
      - a prime predicate [prime]
-     - an efficient [Zgcd] function
+     - properties of the efficient [Zgcd] function defined in [Zgcd_def]
 *)
 
 (** * Divisibility *)
@@ -33,6 +37,11 @@ Notation "( a | b )" := (Zdivide a b) (at level 0) : Z_scope.
 
 (** Results concerning divisibility*)
 
+Lemma Zdivide_equiv : forall a b, Zdivide' a b <-> Zdivide a b.
+Proof.
+ intros a b; split; intros (c,H); exists c; rewrite Zmult_comm; auto.
+Qed.
+
 Lemma Zdivide_refl : forall a:Z, (a | a).
 Proof.
   intros; apply Zdivide_intro with 1; ring.
@@ -427,7 +436,7 @@ Section extended_euclid_algorithm.
   (** The recursive part of Euclid's algorithm uses well-founded
       recursion of non-negative integers. It maintains 6 integers
       [u1,u2,u3,v1,v2,v3] such that the following invariant holds:
-      [u1*a+u2*b=u3] and [v1*a+v2*b=v3] and [gcd(u2,v3)=gcd(a,b)].
+      [u1*a+u2*b=u3] and [v1*a+v2*b=v3] and [gcd(u3,v3)=gcd(a,b)].
       *)
 
   Lemma euclid_rec :
@@ -873,183 +882,17 @@ Proof.
   contradict Hp; auto with zarith.
 Qed.
 
+(** we now prove that Zgcd (defined in Zgcd_def) is indeed a gcd in
+   the sense of Zis_gcd. *)
 
-(** We could obtain a [Zgcd] function via Euclid algorithm. But we propose
-  here a binary version of [Zgcd], faster and executable within Coq.
-
-   Algorithm:
-
-   gcd 0 b = b
-   gcd a 0 = a
-   gcd (2a) (2b) = 2(gcd a b)
-   gcd (2a+1) (2b) = gcd (2a+1) b
-   gcd (2a) (2b+1) = gcd a (2b+1)
-   gcd (2a+1) (2b+1) = gcd (b-a) (2*a+1)
-                    or gcd (a-b) (2*b+1), depending on whether a<b
-*)
-
-Open Scope positive_scope.
-
-Fixpoint Pgcdn (n: nat) (a b : positive) : positive :=
-  match n with
-    | O => 1
-    | S n =>
-      match a,b with
-	| xH, _ => 1
-	| _, xH => 1
-	| xO a, xO b => xO (Pgcdn n a b)
-	| a, xO b => Pgcdn n a b
-	| xO a, b => Pgcdn n a b
-	| xI a', xI b' =>
-          match Pcompare a' b' Eq with
-	    | Eq => a
-	    | Lt => Pgcdn  n (b'-a') a
-	    | Gt => Pgcdn n (a'-b') b
-          end
-      end
-  end.
-
-Definition Pgcd (a b: positive) := Pgcdn (Psize a + Psize b)%nat a b.
-
-Close Scope positive_scope.
-
-Definition Zgcd (a b : Z) : Z :=
-  match a,b with
-    | Z0, _ => Zabs b
-    | _, Z0 => Zabs a
-    | Zpos a, Zpos b => Zpos (Pgcd a b)
-    | Zpos a, Zneg b => Zpos (Pgcd a b)
-    | Zneg a, Zpos b => Zpos (Pgcd a b)
-    | Zneg a, Zneg b => Zpos (Pgcd a b)
-  end.
-
-Lemma Zgcd_is_pos : forall a b, 0 <= Zgcd a b.
-Proof.
-  unfold Zgcd; destruct a; destruct b; auto with zarith.
-Qed.
-
-Lemma Zis_gcd_even_odd : forall a b g, Zis_gcd (Zpos a) (Zpos (xI b)) g ->
-  Zis_gcd (Zpos (xO a)) (Zpos (xI b)) g.
-Proof.
-  intros.
-  destruct H.
-  constructor; auto.
-  destruct H as (e,H2); exists (2*e); auto with zarith.
-  rewrite Zpos_xO; rewrite H2; ring.
-  intros.
-  apply H1; auto.
-  rewrite Zpos_xO in H2.
-  rewrite Zpos_xI in H3.
-  apply Gauss with 2; auto.
-  apply bezout_rel_prime.
-  destruct H3 as (bb, H3).
-  apply Bezout_intro with bb (-Zpos b).
-  omega.
-Qed.
-
-Lemma Pgcdn_correct : forall n a b, (Psize a + Psize b<=n)%nat ->
-  Zis_gcd (Zpos a) (Zpos b) (Zpos (Pgcdn n a b)).
-Proof.
-  intro n; pattern n; apply lt_wf_ind; clear n; intros.
-  destruct n.
-  simpl.
-  destruct a; simpl in *; try inversion H0.
-  destruct a.
-  destruct b; simpl.
-  case_eq (Pcompare a b Eq); intros.
-  (* a = xI, b = xI, compare = Eq *)
-  rewrite (Pcompare_Eq_eq _ _ H1); apply Zis_gcd_refl.
-  (* a = xI, b = xI, compare = Lt *)
-  apply Zis_gcd_sym.
-  apply Zis_gcd_for_euclid with 1.
-  apply Zis_gcd_sym.
-  replace (Zpos (xI b) - 1 * Zpos (xI a)) with (Zpos(xO (b - a))).
-  apply Zis_gcd_even_odd.
-  apply H; auto.
-  simpl in *.
-  assert (Psize (b-a) <= Psize b)%nat.
-  apply Psize_monotone.
-  change (Zpos (b-a) < Zpos b).
-  rewrite (Zpos_minus_morphism _ _ H1).
-  assert (0 < Zpos a) by (compute; auto).
-  omega.
-  omega.
-  rewrite Zpos_xO; do 2 rewrite Zpos_xI.
-  rewrite Zpos_minus_morphism; auto.
-  omega.
-  (* a = xI, b = xI, compare = Gt *)
-  apply Zis_gcd_for_euclid with 1.
-  replace (Zpos (xI a) - 1 * Zpos (xI b)) with (Zpos(xO (a - b))).
-  apply Zis_gcd_sym.
-  apply Zis_gcd_even_odd.
-  apply H; auto.
-  simpl in *.
-  assert (Psize (a-b) <= Psize a)%nat.
-  apply Psize_monotone.
-  change (Zpos (a-b) < Zpos a).
-  rewrite (Zpos_minus_morphism b a).
-  assert (0 < Zpos b) by (compute; auto).
-  omega.
-  rewrite ZC4; rewrite H1; auto.
-  omega.
-  rewrite Zpos_xO; do 2 rewrite Zpos_xI.
-  rewrite Zpos_minus_morphism; auto.
-  omega.
-  rewrite ZC4; rewrite H1; auto.
-  (* a = xI, b = xO *)
-  apply Zis_gcd_sym.
-  apply Zis_gcd_even_odd.
-  apply Zis_gcd_sym.
-  apply H; auto.
-  simpl in *; omega.
-  (* a = xI, b = xH *)
-  apply Zis_gcd_1.
-  destruct b; simpl.
-  (* a = xO, b = xI *)
-  apply Zis_gcd_even_odd.
-  apply H; auto.
-  simpl in *; omega.
-  (* a = xO, b = xO *)
-  rewrite (Zpos_xO a); rewrite (Zpos_xO b); rewrite (Zpos_xO (Pgcdn n a b)).
-  apply Zis_gcd_mult.
-  apply H; auto.
-  simpl in *; omega.
-  (* a = xO, b = xH *)
-  apply Zis_gcd_1.
-  (* a = xH *)
-  simpl; apply Zis_gcd_sym; apply Zis_gcd_1.
-Qed.
-
-Lemma Pgcd_correct : forall a b, Zis_gcd (Zpos a) (Zpos b) (Zpos (Pgcd a b)).
-Proof.
-  unfold Pgcd; intros.
-  apply Pgcdn_correct; auto.
-Qed.
+Notation Zgcd_is_pos := Zgcd_nonneg (only parsing).
 
 Lemma Zgcd_is_gcd : forall a b, Zis_gcd a b (Zgcd a b).
 Proof.
-  destruct a.
-  intros.
-  simpl.
-  apply Zis_gcd_0_abs.
-  destruct b; simpl.
-  apply Zis_gcd_0.
-  apply Pgcd_correct.
-  apply Zis_gcd_sym.
-  apply Zis_gcd_minus; simpl.
-  apply Pgcd_correct.
-  destruct b; simpl.
-  apply Zis_gcd_minus; simpl.
-  apply Zis_gcd_sym.
-  apply Zis_gcd_0.
-  apply Zis_gcd_minus; simpl.
-  apply Zis_gcd_sym.
-  apply Pgcd_correct.
-  apply Zis_gcd_sym.
-  apply Zis_gcd_minus; simpl.
-  apply Zis_gcd_minus; simpl.
-  apply Zis_gcd_sym.
-  apply Pgcd_correct.
+ constructor; intros; apply Zdivide_equiv.
+ apply Zgcd_divide_l.
+ apply Zgcd_divide_r.
+ apply Zgcd_greatest; now apply Zdivide_equiv.
 Qed.
 
 Theorem Zgcd_spec : forall x y : Z, {z : Z | Zis_gcd x y z /\ 0 <= z}.
@@ -1247,167 +1090,3 @@ Proof.
   absurd (n = 0); auto with zarith.
   case Hr1; auto with zarith.
 Qed.
-
-(** A Generalized Gcd that also computes Bezout coefficients.
-   The algorithm is the same as for Zgcd. *)
-
-Open Scope positive_scope.
-
-Fixpoint Pggcdn (n: nat) (a b : positive) : (positive*(positive*positive)) :=
-  match n with
-    | O => (1,(a,b))
-    | S n =>
-      match a,b with
-	| xH, b => (1,(1,b))
-	| a, xH => (1,(a,1))
-	| xO a, xO b =>
-           let (g,p) := Pggcdn n a b in
-           (xO g,p)
-	| a, xO b =>
-           let (g,p) := Pggcdn n a b in
-           let (aa,bb) := p in
-           (g,(aa, xO bb))
-	| xO a, b =>
-           let (g,p) := Pggcdn n a b in
-           let (aa,bb) := p in
-           (g,(xO aa, bb))
-	| xI a', xI b' =>
-           match Pcompare a' b' Eq with
-	     | Eq => (a,(1,1))
-	     | Lt =>
-	        let (g,p) := Pggcdn n (b'-a') a in
-	        let (ba,aa) := p in
-	        (g,(aa, aa + xO ba))
-	     | Gt =>
-		let (g,p) := Pggcdn n (a'-b') b in
-		let (ab,bb) := p in
-		(g,(bb+xO ab, bb))
-	   end
-      end
-  end.
-
-Definition Pggcd (a b: positive) := Pggcdn (Psize a + Psize b)%nat a b.
-
-Open Scope Z_scope.
-
-Definition Zggcd (a b : Z) : Z*(Z*Z) :=
-  match a,b with
-    | Z0, _ => (Zabs b,(0, Zsgn b))
-    | _, Z0 => (Zabs a,(Zsgn a, 0))
-    | Zpos a, Zpos b =>
-       let (g,p) := Pggcd a b in
-       let (aa,bb) := p in
-       (Zpos g, (Zpos aa, Zpos bb))
-    | Zpos a, Zneg b =>
-       let (g,p) := Pggcd a b in
-       let (aa,bb) := p in
-       (Zpos g, (Zpos aa, Zneg bb))
-    | Zneg a, Zpos b =>
-       let (g,p) := Pggcd a b in
-       let (aa,bb) := p in
-       (Zpos g, (Zneg aa, Zpos bb))
-    | Zneg a, Zneg b =>
-       let (g,p) := Pggcd a b in
-       let (aa,bb) := p in
-       (Zpos g, (Zneg aa, Zneg bb))
-  end.
-
-
-Lemma Pggcdn_gcdn : forall n a b,
-  fst (Pggcdn n a b) = Pgcdn n a b.
-Proof.
-  induction n.
-  simpl; auto.
-  destruct a; destruct b; simpl; auto.
-  destruct (Pcompare a b Eq); simpl; auto.
-  rewrite <- IHn; destruct (Pggcdn n (b-a) (xI a)) as (g,(aa,bb)); simpl; auto.
-  rewrite <- IHn; destruct (Pggcdn n (a-b) (xI b)) as (g,(aa,bb)); simpl; auto.
-  rewrite <- IHn; destruct (Pggcdn n (xI a) b) as (g,(aa,bb)); simpl; auto.
-  rewrite <- IHn; destruct (Pggcdn n a (xI b)) as (g,(aa,bb)); simpl; auto.
-  rewrite <- IHn; destruct (Pggcdn n a b) as (g,(aa,bb)); simpl; auto.
-Qed.
-
-Lemma Pggcd_gcd : forall a b, fst (Pggcd a b) = Pgcd a b.
-Proof.
-  intros; exact (Pggcdn_gcdn (Psize a+Psize b)%nat a b).
-Qed.
-
-Lemma Zggcd_gcd : forall a b, fst (Zggcd a b) = Zgcd a b.
-Proof.
-  destruct a; destruct b; simpl; auto; rewrite <- Pggcd_gcd;
-    destruct (Pggcd p p0) as (g,(aa,bb)); simpl; auto.
-Qed.
-
-Open Scope positive_scope.
-
-Lemma Pggcdn_correct_divisors : forall n a b,
-  let (g,p) := Pggcdn n a b in
-  let (aa,bb):=p in
-  (a=g*aa) /\ (b=g*bb).
-Proof.
-  induction n.
-  simpl; auto.
-  destruct a; destruct b; simpl; auto.
-  case_eq (Pcompare a b Eq); intros.
-  (* Eq *)
-  rewrite Pmult_comm; simpl; auto.
-  rewrite (Pcompare_Eq_eq _ _ H); auto.
-  (* Lt *)
-  generalize (IHn (b-a) (xI a)); destruct (Pggcdn n (b-a) (xI a)) as (g,(ba,aa)); simpl.
-  intros (H0,H1); split; auto.
-  rewrite Pmult_plus_distr_l.
-  rewrite Pmult_xO_permute_r.
-  rewrite <- H1; rewrite <- H0.
-  simpl; f_equal; symmetry.
-  apply Pplus_minus; auto.
-  rewrite ZC4; rewrite H; auto.
-  (* Gt *)
-  generalize (IHn (a-b) (xI b)); destruct (Pggcdn n (a-b) (xI b)) as (g,(ab,bb)); simpl.
-  intros (H0,H1); split; auto.
-  rewrite Pmult_plus_distr_l.
-  rewrite Pmult_xO_permute_r.
-  rewrite <- H1; rewrite <- H0.
-  simpl; f_equal; symmetry.
-  apply Pplus_minus; auto.
-  (* Then... *)
-  generalize (IHn (xI a) b); destruct (Pggcdn n (xI a) b) as (g,(ab,bb)); simpl.
-  intros (H0,H1); split; auto.
-  rewrite Pmult_xO_permute_r; rewrite H1; auto.
-  generalize (IHn a (xI b)); destruct (Pggcdn n a (xI b)) as (g,(ab,bb)); simpl.
-  intros (H0,H1); split; auto.
-  rewrite Pmult_xO_permute_r; rewrite H0; auto.
-  generalize (IHn a b); destruct (Pggcdn n a b) as (g,(ab,bb)); simpl.
-  intros (H0,H1); split; subst; auto.
-Qed.
-
-Lemma Pggcd_correct_divisors : forall a b,
-  let (g,p) := Pggcd a b in
-  let (aa,bb):=p in
-  (a=g*aa) /\ (b=g*bb).
-Proof.
-  intros a b; exact (Pggcdn_correct_divisors (Psize a + Psize b)%nat a b).
-Qed.
-
-Close Scope positive_scope.
-
-Lemma Zggcd_correct_divisors : forall a b,
-  let (g,p) := Zggcd a b in
-  let (aa,bb):=p in
-  (a=g*aa) /\ (b=g*bb).
-Proof.
-  destruct a; destruct b; simpl; auto; try solve [rewrite Pmult_comm; simpl; auto];
-    generalize (Pggcd_correct_divisors p p0); destruct (Pggcd p p0) as (g,(aa,bb));
-      destruct 1; subst; auto.
-Qed.
-
-Theorem Zggcd_opp: forall x y,
-  Zggcd (-x) y = let (p1,p) := Zggcd x y in
-                 let (p2,p3) := p in
-                 (p1,(-p2,p3)).
-Proof.
-intros [|x|x] [|y|y]; unfold Zggcd, Zopp; auto.
-case Pggcd; intros p1 (p2, p3); auto.
-case Pggcd; intros p1 (p2, p3); auto.
-case Pggcd; intros p1 (p2, p3); auto.
-case Pggcd; intros p1 (p2, p3); auto.
-Qed.
diff --git a/theories/ZArith/vo.itarget b/theories/ZArith/vo.itarget
index 7802a8655..bfc52d3e4 100644
--- a/theories/ZArith/vo.itarget
+++ b/theories/ZArith/vo.itarget
@@ -31,3 +31,4 @@ Zsqrt_compat.vo
 Zwf.vo
 Zsqrt_def.vo
 Zlog_def.vo
+Zgcd_def.vo
\ No newline at end of file