nouvel algorithme pour Zgcd (plus rapide) + un Qcompare

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

1 files changed, 382 insertions, 497 deletions

diff --git a/theories/ZArith/Znumtheory.v b/theories/ZArith/Znumtheory.v
index 3a610226e..d61cc84bc 100644
--- a/theories/ZArith/Znumtheory.v
+++ b/theories/ZArith/Znumtheory.v
@@ -12,15 +12,17 @@ Require Import ZArith_base.
 Require Import ZArithRing.
 Require Import Zcomplements.
 Require Import Zdiv.
+Require Import Ndigits.
+Require Import Wf_nat.
 Open Local Scope Z_scope.
 
 (** This file contains some notions of number theory upon Z numbers: 
      - a divisibility predicate [Zdivide]
      - a gcd predicate [gcd]
      - Euclid algorithm [euclid]
-     - an efficient [Zgcd] function 
      - a relatively prime predicate [rel_prime]
      - a prime predicate [prime]
+     - an efficient [Zgcd] function 
 *)
 
 (** * Divisibility *)
@@ -215,6 +217,16 @@ Proof.
 constructor; auto with zarith.
 Qed.
 
+Lemma Zis_gcd_1 : forall a, Zis_gcd a 1 1.
+Proof.
+constructor; auto with zarith.
+Qed.
+
+Lemma Zis_gcd_refl : forall a, Zis_gcd a a a.
+Proof.
+constructor; auto with zarith.
+Qed.
+
 Lemma Zis_gcd_minus : forall a b d:Z, Zis_gcd a (- b) d -> Zis_gcd b a d.
 Proof.
 simple induction 1; constructor; intuition.
@@ -225,6 +237,14 @@ Proof.
 simple induction 1; constructor; intuition.
 Qed.
 
+Lemma Zis_gcd_0_abs : forall a:Z, Zis_gcd 0 a (Zabs a).
+Proof.
+intros a.
+apply Zabs_ind.
+intros; apply Zis_gcd_sym; apply Zis_gcd_0; auto.
+intros; apply Zis_gcd_opp; apply Zis_gcd_0; auto.
+Qed.
+
 Hint Resolve Zis_gcd_sym Zis_gcd_0 Zis_gcd_minus Zis_gcd_opp: zarith.
 
 (** * Extended Euclid algorithm. *)
@@ -366,478 +386,8 @@ replace (c * (u * a + v * b)) with (u * (c * a) + v * (c * b)).
 rewrite H6; rewrite H7; ring.
 ring.
 Qed.
-
-Lemma Zis_gcd_0_abs : forall b, 
- Zis_gcd 0 b (Zabs b) /\ Zabs b >= 0 /\ 0 = Zabs b * 0 /\ b = Zabs b * Zsgn b.
-Proof.
-intro b.
-elim (Z_le_lt_eq_dec _ _ (Zabs_pos b)).
-intros H0; split.
-apply Zabs_ind.
-intros; apply Zis_gcd_sym; apply Zis_gcd_0; auto.
-intros; apply Zis_gcd_opp; apply Zis_gcd_0; auto.
-repeat split; auto with zarith.
-symmetry; apply Zabs_Zsgn.
-
-intros H0; rewrite <- H0.
-rewrite <- (Zabs_Zsgn b); rewrite <- H0; simpl in |- *.
-split; [ apply Zis_gcd_0 | idtac ]; auto with zarith.
-Qed.
   
 
-(** We could obtain a [Zgcd] function via [euclid]. But we propose 
-  here a more direct version of a [Zgcd], that can compute within Coq.
-  For that, we use an explicit measure in [nat], and we proved later
-  that using [2(d+1)] is enough, where  [d] is the number of binary digits 
-  of the first argument. *)
-
-Fixpoint Zgcdn (n:nat) : Z -> Z -> Z := fun a b => 
- match n with 
-  | O => 1 (* arbitrary, since n should be big enough *)
-  | S n => match a with 
-     | Z0 => Zabs b
-     | Zpos _ => Zgcdn n (Zmod b a) a
-     | Zneg a => Zgcdn n (Zmod b (Zpos a)) (Zpos a)
-     end
-   end.
-
-(* For technical reason, we don't use [Ndigit.Psize] but this 
-   ad-hoc version: [Psize p = S (Psiz p)]. *)
-
-Fixpoint Psiz (p:positive) : nat := 
-  match p with 
-    | xH => O
-    | xI p => S (Psiz p) 
-    | xO p => S (Psiz p)
-  end.
-
-Definition Zgcd_bound (a:Z) := match a with 
- | Z0 => S O
- | Zpos p => let n := Psiz p in S (S (n+n))
- | Zneg p => let n := Psiz p in S (S (n+n))
-end.
-
-Definition Zgcd a b := Zgcdn (Zgcd_bound a) a b.
-
-(** A first obvious fact : [Zgcd a b] is positive. *)
-
-Lemma Zgcdn_is_pos : forall n a b, 
-  0 <= Zgcdn n a b.
-Proof.
-induction n.
-simpl; auto with zarith.
-destruct a; simpl; intros; auto with zarith; auto.
-Qed.
-
-Lemma Zgcd_is_pos : forall a b, 0 <= Zgcd a b.
-Proof.
-intros; unfold Zgcd; apply Zgcdn_is_pos; auto.
-Qed.
-
-(** We now prove that Zgcd is indeed a gcd. *)
-
-(** 1) We prove a weaker & easier bound. *)
-
-Lemma Zgcdn_linear_bound : forall n a b, 
- Zabs a < Z_of_nat n -> Zis_gcd a b (Zgcdn n a b).
-Proof.
-induction n.
-simpl; intros.
-elimtype False; generalize (Zabs_pos a); omega.
-destruct a; intros; simpl; 
- [ generalize (Zis_gcd_0_abs b); intuition | | ]; 
- unfold Zmod; 
- generalize (Z_div_mod b (Zpos p) (refl_equal Gt));
- destruct (Zdiv_eucl b (Zpos p)) as (q,r);
- intros (H0,H1); 
- rewrite inj_S in H; simpl Zabs in H; 
- assert (H2: Zabs r < Z_of_nat n) by (rewrite Zabs_eq; auto with zarith);  
- assert (IH:=IHn r (Zpos p) H2); clear IHn; 
- simpl in IH |- *; 
- rewrite H0.
- apply Zis_gcd_for_euclid2; auto.
- apply Zis_gcd_minus; apply Zis_gcd_sym.
- apply Zis_gcd_for_euclid2; auto.
-Qed.
-
-(** 2) For Euclid's algorithm, the worst-case situation corresponds 
-     to Fibonacci numbers. Let's define them: *)
-
-Fixpoint fibonacci (n:nat) : Z := 
- match n with 
-   | O => 1
-   | S O => 1
-   | S (S n as p) => fibonacci p + fibonacci n 
- end.
-
-Lemma fibonacci_pos : forall n, 0 <= fibonacci n.
-Proof.
-cut (forall N n, (n<N)%nat -> 0<=fibonacci n).
-eauto.
-induction N.
-inversion 1.
-intros.
-destruct n.
-simpl; auto with zarith.
-destruct n.
-simpl; auto with zarith.
-change (0 <= fibonacci (S n) + fibonacci n).
-generalize (IHN n) (IHN (S n)); omega.
-Qed.
-
-Lemma fibonacci_incr : 
- forall n m, (n<=m)%nat -> fibonacci n <= fibonacci m.
-Proof.
-induction 1.
-auto with zarith.
-apply Zle_trans with (fibonacci m); auto.
-clear.
-destruct m.
-simpl; auto with zarith.
-change (fibonacci (S m) <= fibonacci (S m)+fibonacci m).
-generalize (fibonacci_pos m); omega.
-Qed.
-
-(** 3) We prove that fibonacci numbers are indeed worst-case: 
-   for a given number [n], if we reach a conclusion about [gcd(a,b)] in 
-   exactly [n+1] loops, then [fibonacci (n+1)<=a /\ fibonacci(n+2)<=b] *)
-
-Lemma Zgcdn_worst_is_fibonacci : forall n a b, 
- 0 < a < b -> 
- Zis_gcd a b (Zgcdn (S n) a b) -> 
- Zgcdn n a b <> Zgcdn (S n) a b -> 
- fibonacci (S n) <= a /\ 
- fibonacci (S (S n)) <= b.
-Proof.
-induction n.
-simpl; intros.
-destruct a; omega.
-intros.
-destruct a; [simpl in *; omega| | destruct H; discriminate].
-revert H1; revert H0.
-set (m:=S n) in *; (assert (m=S n) by auto); clearbody m.
-pattern m at 2; rewrite H0.
-simpl Zgcdn.
-unfold Zmod; generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
-destruct (Zdiv_eucl b (Zpos p)) as (q,r).
-intros (H1,H2).
-destruct H2.
-destruct (Zle_lt_or_eq _ _ H2).
-generalize (IHn _ _ (conj H4 H3)).
-intros H5 H6 H7. 
-replace (fibonacci (S (S m))) with (fibonacci (S m) + fibonacci m) by auto.
-assert (r = Zpos p * (-q) + b) by (rewrite H1; ring).
-destruct H5; auto.
-pattern r at 1; rewrite H8.
-apply Zis_gcd_sym.
-apply Zis_gcd_for_euclid2; auto.
-apply Zis_gcd_sym; auto.
-split; auto.
-rewrite H1.
-apply Zplus_le_compat; auto.
-apply Zle_trans with (Zpos p * 1); auto.
-ring (Zpos p * 1); auto.
-apply Zmult_le_compat_l.
-destruct q.
-omega.
-assert (0 < Zpos p0) by (compute; auto).
-omega.
-assert (Zpos p * Zneg p0 < 0) by (compute; auto).
-omega. 
-compute; intros; discriminate.
-(* r=0 *)
-subst r.
-simpl; rewrite H0.
-intros.
-simpl in H4.
-simpl in H5.
-destruct n.
-simpl in H5.
-simpl.
-omega.
-simpl in H5.
-elim H5; auto.
-Qed.
-
-(** 3b) We reformulate the previous result in a more positive way. *)
-
-Lemma Zgcdn_ok_before_fibonacci : forall n a b, 
- 0 < a < b -> a < fibonacci (S n) ->
- Zis_gcd a b (Zgcdn n a b).
-Proof.
-destruct a; [ destruct 1; elimtype False; omega | | destruct 1; discriminate].
-cut (forall k n b, 
-         k = (S (nat_of_P p) - n)%nat -> 
-         0 < Zpos p < b -> Zpos p < fibonacci (S n) -> 
-         Zis_gcd (Zpos p) b (Zgcdn n (Zpos p) b)).
-destruct 2; eauto.
-clear n; induction k. 
-intros.
-assert (nat_of_P p < n)%nat by omega.
-apply Zgcdn_linear_bound. 
-simpl.
-generalize (inj_le _ _ H2).
-rewrite inj_S.
-rewrite <- Zpos_eq_Z_of_nat_o_nat_of_P; auto.
-omega.
-intros.
-generalize (Zgcdn_worst_is_fibonacci n (Zpos p) b H0); intros.
-assert (Zis_gcd (Zpos p) b (Zgcdn (S n) (Zpos p) b)).
- apply IHk; auto.
- omega.
- replace (fibonacci (S (S n))) with (fibonacci (S n)+fibonacci n) by auto.
- generalize (fibonacci_pos n); omega.
-replace (Zgcdn n (Zpos p) b) with (Zgcdn (S n) (Zpos p) b); auto.
-generalize (H2 H3); clear H2 H3; omega.
-Qed.
-
-(** 4) The proposed bound leads to a fibonacci number that is big enough. *)
-
-Lemma Zgcd_bound_fibonacci : 
-  forall a, 0 < a -> a < fibonacci (Zgcd_bound a).
-Proof.
-destruct a; [omega| | intro H; discriminate].
-intros _.
-induction p.
-simpl Zgcd_bound in *.
-rewrite Zpos_xI.
-rewrite plus_comm; simpl plus.
-set (n:=S (Psiz p+Psiz p)) in *.
-change (2*Zpos p+1 < 
- fibonacci (S n) + fibonacci n + fibonacci (S n)).
-generalize (fibonacci_pos n).
-omega.
-simpl Zgcd_bound in *.
-rewrite Zpos_xO.
-rewrite plus_comm; simpl plus.
-set (n:= S (Psiz p +Psiz p)) in *.
-change (2*Zpos p < 
- fibonacci (S n) + fibonacci n + fibonacci (S n)).
-generalize (fibonacci_pos n).
-omega.
-simpl; auto with zarith.
-Qed.
-
-(* 5) the end: we glue everything together and take care of 
-   situations not corresponding to [0<a<b]. *)
-
-Lemma Zgcd_is_gcd : 
-  forall a b, Zis_gcd a b (Zgcd a b).
-Proof.
-unfold Zgcd; destruct a; intros.
-simpl; generalize (Zis_gcd_0_abs b); intuition.
-(*Zpos*)
-generalize (Zgcd_bound_fibonacci (Zpos p)).
-simpl Zgcd_bound.
-set (n:=S (Psiz p+Psiz p)); (assert (n=S (Psiz p+Psiz p)) by auto); clearbody n.
-simpl Zgcdn.
-unfold Zmod.
-generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
-destruct (Zdiv_eucl b (Zpos p)) as (q,r).
-intros (H1,H2) H3.
-rewrite H1.
-apply Zis_gcd_for_euclid2.
-destruct H2.
-destruct (Zle_lt_or_eq _ _ H0).
-apply Zgcdn_ok_before_fibonacci; auto; omega.
-subst r n; simpl.
-apply Zis_gcd_sym; apply Zis_gcd_0.
-(*Zneg*)
-generalize (Zgcd_bound_fibonacci (Zpos p)).
-simpl Zgcd_bound.
-set (n:=S (Psiz p+Psiz p)); (assert (n=S (Psiz p+Psiz p)) by auto); clearbody n.
-simpl Zgcdn.
-unfold Zmod.
-generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
-destruct (Zdiv_eucl b (Zpos p)) as (q,r).
-intros (H1,H2) H3.
-rewrite H1.
-apply Zis_gcd_minus.
-apply Zis_gcd_sym.
-apply Zis_gcd_for_euclid2.
-destruct H2.
-destruct (Zle_lt_or_eq _ _ H0).
-apply Zgcdn_ok_before_fibonacci; auto; omega.
-subst r n; simpl.
-apply Zis_gcd_sym; apply Zis_gcd_0.
-Qed.
-
-(** A generalized gcd: it additionnally keeps track of the divisors. *)
-
-Fixpoint Zggcdn (n:nat) : Z -> Z -> (Z*(Z*Z)) := fun a b => 
- match n with 
-  | O => (1,(a,b))   (*(Zabs b,(0,Zsgn b))*)
-  | S n => match a with 
-     | Z0 => (Zabs b,(0,Zsgn b))
-     | Zpos _ => 
-               let (q,r) := Zdiv_eucl b a in   (* b = q*a+r *)
-               let (g,p) := Zggcdn n r a in 
-               let (rr,aa) := p in        (* r = g *rr /\ a = g * aa *)
-               (g,(aa,q*aa+rr))
-     | Zneg a => 
-               let (q,r) := Zdiv_eucl b (Zpos a) in  (* b = q*(-a)+r *)
-               let (g,p) := Zggcdn n r (Zpos a) in 
-               let (rr,aa) := p in         (* r = g*rr /\ (-a) = g * aa *)  
-               (g,(-aa,q*aa+rr))
-     end
-   end.
-
-Definition Zggcd a b : Z * (Z * Z) := Zggcdn (Zgcd_bound a) a b.
-
-(** The first component of [Zggcd] is [Zgcd] *) 
-
-Lemma Zggcdn_gcdn : forall n a b, 
- fst (Zggcdn n a b) = Zgcdn n a b.
-Proof.
-induction n; simpl; auto.
-destruct a; unfold Zmod; simpl; intros; auto; 
- destruct (Zdiv_eucl b (Zpos p)) as (q,r); 
- rewrite <- IHn; 
- destruct (Zggcdn n r (Zpos p)) as (g,(rr,aa)); simpl; auto.
-Qed.
-
-Lemma Zggcd_gcd : forall a b, fst (Zggcd a b) = Zgcd a b.
-Proof.
-intros; unfold Zggcd, Zgcd; apply Zggcdn_gcdn; auto.
-Qed.
-
-(** [Zggcd] always returns divisors that are coherent with its 
- first output. *)
-
-Lemma Zggcdn_correct_divisors : forall n a b, 
-  let (g,p) := Zggcdn n a b in 
-  let (aa,bb):=p in 
-  a=g*aa /\ b=g*bb.
-Proof.
-induction n.
-simpl.
-split; [destruct a|destruct b]; auto.
-intros.
-simpl.
-destruct a.
-rewrite Zmult_comm; simpl.
-split; auto.
-symmetry; apply Zabs_Zsgn.
-generalize (Z_div_mod b (Zpos p)); 
-destruct (Zdiv_eucl b (Zpos p)) as (q,r).
-generalize (IHn r (Zpos p)); 
-destruct (Zggcdn n r (Zpos p)) as (g,(rr,aa)).
-intuition.
-destruct H0.
-compute; auto.
-rewrite H; rewrite H1; rewrite H2; ring.
-generalize (Z_div_mod b (Zpos p)); 
-destruct (Zdiv_eucl b (Zpos p)) as (q,r).
-destruct 1.
-compute; auto.
-generalize (IHn r (Zpos p)); 
-destruct (Zggcdn n r (Zpos p)) as (g,(rr,aa)).
-intuition.
-destruct H0.
-replace (Zneg p) with (-Zpos p) by compute; auto.
-rewrite H4; ring.
-rewrite H; rewrite H4; rewrite H0; ring.
-Qed.
-
-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.
-unfold Zggcd; intros; apply Zggcdn_correct_divisors; auto.
-Qed.
- 
-(** Due to the use of an explicit measure, the extraction of [Zgcd] 
-  isn't optimal. We propose here another version [Zgcd_spec] that 
-  doesn't suffer from this problem (but doesn't compute in Coq). *)
-
-Definition Zgcd_spec_pos :
-  forall a:Z,
-    0 <= a -> forall b:Z, {g : Z | 0 <= a -> Zis_gcd a b g /\ g >= 0}.
-Proof.
-intros a Ha.
-apply
- (Zlt_0_rec
-    (fun a:Z => forall b:Z, {g : Z | 0 <= a -> Zis_gcd a b g /\ g >= 0}));
- try assumption.
-intro x; case x.
-intros _ _ b; exists (Zabs b).
-generalize (Zis_gcd_0_abs b); intuition.
-  
-intros p Hrec _ b.
-generalize (Z_div_mod b (Zpos p)).
-case (Zdiv_eucl b (Zpos p)); intros q r Hqr.
-elim Hqr; clear Hqr; intros; auto with zarith.
-elim (Hrec r H0 (Zpos p)); intros g Hgkl.
-inversion_clear H0.
-elim (Hgkl H1); clear Hgkl; intros H3 H4.
-exists g; intros.
-split; auto.
-rewrite H.
-apply Zis_gcd_for_euclid2; auto.
-
-intros p _ H b.
-elim H; auto.
-Defined.
-
-Definition Zgcd_spec : forall a b:Z, {g : Z | Zis_gcd a b g /\ g >= 0}.
-Proof.
-intros a; case (Z_gt_le_dec 0 a).
-intros; assert (0 <= - a).
-omega.
-elim (Zgcd_spec_pos (- a) H b); intros g Hgkl.
-exists g.
-intuition.
-intros Ha b; elim (Zgcd_spec_pos a Ha b); intros g; exists g; intuition.
-Defined.
-
-(** A last version aimed at extraction that also returns the divisors. *)
-
-Definition Zggcd_spec_pos :
-  forall a:Z,
-    0 <= a -> forall b:Z, {p : Z*(Z*Z) | let (g,p):=p in let (aa,bb):=p in 
-                                       0 <= a -> Zis_gcd a b g /\ g >= 0 /\ a=g*aa /\ b=g*bb}.
-Proof.
-intros a Ha.
-pattern a; apply Zlt_0_rec; try assumption.
-intro x; case x.
-intros _ _ b; exists (Zabs b,(0,Zsgn b)).
-intros _; apply Zis_gcd_0_abs.
-  
-intros p Hrec _ b.
-generalize (Z_div_mod b (Zpos p)).
-case (Zdiv_eucl b (Zpos p)); intros q r Hqr.
-elim Hqr; clear Hqr; intros; auto with zarith.
-destruct (Hrec r H0 (Zpos p)) as ((g,(rr,pp)),Hgkl).
-destruct H0.
-destruct (Hgkl H0) as (H3,(H4,(H5,H6))).
-exists (g,(pp,pp*q+rr)); intros.
-split; auto.
-rewrite H.
-apply Zis_gcd_for_euclid2; auto.
-repeat split; auto.
-rewrite H; rewrite H6; rewrite H5; ring.
-
-intros p _ H b.
-elim H; auto.
-Defined.
-
-Definition Zggcd_spec : 
- forall a b:Z, {p : Z*(Z*Z) |  let (g,p):=p in let (aa,bb):=p in 
-                                          Zis_gcd a b g /\ g >= 0 /\ a=g*aa /\ b=g*bb}.
-Proof.
-intros a; case (Z_gt_le_dec 0 a).
-intros; assert (0 <= - a).
-omega.
-destruct (Zggcd_spec_pos (- a) H b) as ((g,(aa,bb)),Hgkl).
-exists (g,(-aa,bb)).
-intuition.
-rewrite <- Zopp_mult_distr_r.
-rewrite <- H2; auto with zarith.
-intros Ha b; elim (Zggcd_spec_pos a Ha b); intros p; exists p.
- repeat destruct p; intuition.
-Defined.
-
 (** * Relative primality *)
 
 Definition rel_prime (a b:Z) : Prop := Zis_gcd a b 1.
@@ -920,32 +470,25 @@ assert (g <> 0).
  elim H4; intros.
  rewrite H2 in H6; subst b; omega.
 unfold rel_prime in |- *.
-elim (Zgcd_spec (a / g) (b / g)); intros g' [H3 H4].
-assert (H5 := Zis_gcd_mult _ _ g _ H3).
-rewrite <- Z_div_exact_2 in H5; auto with zarith.
-rewrite <- Z_div_exact_2 in H5; auto with zarith.
-elim (Zis_gcd_uniqueness_apart_sign _ _ _ _ H1 H5).
-intros; rewrite (Zmult_reg_l 1 g' g); auto with zarith.
-intros; rewrite (Zmult_reg_l 1 (- g') g); auto with zarith.
-pattern g at 1 in |- *; rewrite H6; ring.
-
-elim H1; intros.
-elim H7; intros.
-rewrite H9.
-replace (q * g) with (0 + q * g).
-rewrite Z_mod_plus.
-compute in |- *; auto.
-omega.
-ring.
-
-elim H1; intros.
-elim H6; intros.
-rewrite H9.
-replace (q * g) with (0 + q * g).
-rewrite Z_mod_plus.
-compute in |- *; auto.
-omega.
-ring.
+destruct H1.
+destruct H1 as (a',H1).
+destruct H3 as (b',H3).
+replace (a/g) with a'; 
+ [|rewrite H1; rewrite Z_div_mult; auto with zarith].
+replace (b/g) with b'; 
+ [|rewrite H3; rewrite Z_div_mult; auto with zarith].
+constructor.
+exists a'; auto with zarith.
+exists b'; auto with zarith.
+intros x (xa,H5) (xb,H6).
+destruct (H4 (x*g)).
+exists xa; rewrite Zmult_assoc; rewrite <- H5; auto.
+exists xb; rewrite Zmult_assoc; rewrite <- H6; auto.
+replace g with (1*g) in H7; auto with zarith.
+do 2 rewrite Zmult_assoc in H7.
+generalize (Zmult_reg_r _ _ _ H2 H7); clear H7; intros.
+rewrite Zmult_1_r in H7.
+exists q; auto with zarith.
 Qed.
 
 (** * Primality *)
@@ -1045,3 +588,345 @@ case (Zdivide_dec p a); intuition.
 right; apply Gauss with a; auto with zarith.
 Qed.
 
+
+(** 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) { struct n } : 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.
+
+Fixpoint Pggcdn (n: nat) (a b : positive) { struct n } : (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 Pgcd (a b: positive) := Pgcdn (Psize a + Psize b)%nat a b.
+Definition Pggcd (a b: positive) := Pggcdn (Psize a + Psize b)%nat a b.
+
+Open Scope Z_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.
+
+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 Zgcd_is_pos : forall a b, 0 <= Zgcd a b.
+Proof.
+unfold Zgcd; destruct a; destruct b; auto with zarith.
+Qed.
+
+Lemma Psize_monotone : forall p q, Pcompare p q Eq = Lt -> (Psize p <= Psize q)%nat.
+Proof.
+induction p; destruct q; simpl; auto with arith; intros; try discriminate.
+intros; generalize (Pcompare_Gt_Lt _ _ H); auto with arith.
+intros; destruct (Pcompare_Lt_Lt _ _ H); auto with arith; subst; auto.
+Qed.
+
+Lemma Pminus_Zminus : forall a b, Pcompare a b Eq = Lt -> 
+ Zpos (b-a) = Zpos b - Zpos a.
+Proof.
+intros.
+repeat rewrite Zpos_eq_Z_of_nat_o_nat_of_P.
+rewrite nat_of_P_minus_morphism.
+apply inj_minus1.
+apply lt_le_weak.
+apply nat_of_P_lt_Lt_compare_morphism; auto.
+rewrite ZC4; rewrite H; auto.
+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 (Pminus_Zminus _ _ H1).
+ assert (0 < Zpos a) by (compute; auto).
+ omega.
+omega.  
+rewrite Zpos_xO; do 2 rewrite Zpos_xI.
+rewrite Pminus_Zminus; 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 (Pminus_Zminus 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 Pminus_Zminus; 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.
+
+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.
+Qed.
+
+
+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.
+
+Open Scope Z_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.
+
+(** A version of [Zgcd] that doesn't use an explicit measure can be found 
+  in users's contribution [Orsay/QArith]. It is slightly more efficient after 
+  extraction, but cannot be used to compute within Coq. *)
+