(* This program is free software; you can redistribute it and/or      *)
(* modify it under the terms of the GNU Lesser General Public License *)
(* as published by the Free Software Foundation; either version 2.1   *)
(* of the License, or (at your option) any later version.             *)
(*                                                                    *)
(* This program is distributed in the hope that it will be useful,    *)
(* but WITHOUT ANY WARRANTY; without even the implied warranty of     *)
(* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *)
(* GNU General Public License for more details.                       *)
(*                                                                    *)
(* You should have received a copy of the GNU Lesser General Public   *)
(* License along with this program; if not, write to the Free         *)
(* Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA *)
(* 02110-1301 USA                                                     *)



Require Import ZArith.
Require Import ZArithRing.
Require Import Zdiv.
Require Import Omega.

Require Import divide.

(** * Greatest common divisor (gcd). *)
   
(** There is no unicity of the gcd; hence we define the
    predicate [gcd a b d] expressing that [d] is a gcd of 
    [a] and [b]. (We show later that the [gcd] is actually
    unique if we discard its sign.) *)

Inductive gcd (a b d : Z) : Prop :=
    gcd_intro :
      (d | a)%Z ->
      (d | b)%Z ->
      (forall x : Z, (x | a)%Z -> (x | b)%Z -> (x | d)%Z) -> gcd a b d.

(** Trivial properties of [gcd] *)

Lemma gcd_sym : forall a b d : Z, gcd a b d -> gcd b a d.
Proof.
(* Goal: forall (a b c d : Z) (_ : gcd a b d), gcd (Z.mul c a) (Z.mul c b) (Z.mul c d) *)
simple induction 1; constructor; intuition.
Qed.

Lemma gcd_0 : forall a : Z, gcd a 0 a.
Proof.
(* Goal: forall a : Z, gcd a Z0 a *)
constructor; auto.
Qed.

Lemma gcd_minus : forall a b d : Z, gcd a (- b) d -> gcd b a d.
Proof.
(* Goal: forall (a b c d : Z) (_ : gcd a b d), gcd (Z.mul c a) (Z.mul c b) (Z.mul c d) *)
simple induction 1; constructor; intuition.
Qed.

Lemma gcd_opp : forall a b d : Z, gcd a b d -> gcd b a (- d).
Proof.
(* Goal: forall (a b c d : Z) (_ : gcd a b d), gcd (Z.mul c a) (Z.mul c b) (Z.mul c d) *)
simple induction 1; constructor; intuition.
Qed.

Hint Resolve gcd_sym gcd_0 gcd_minus gcd_opp.

(** * Extended Euclid algorithm. *)

(** Euclid's algorithm to compute the [gcd] mainly relies on
    the following property. *)

Lemma gcd_for_euclid : forall a b d q : Z, gcd b (a - q * b) d -> gcd a b d.
Proof.
(* Goal: forall (a b c d : Z) (_ : gcd a b d), gcd (Z.mul c a) (Z.mul c b) (Z.mul c d) *)
simple induction 1; constructor; intuition.
(* Goal: divide d a *)
replace a with (a - q * b + q * b)%Z. auto. ring.
Qed.

(** We implement the extended version of Euclid's algorithm,
    i.e. the one computing Bezout's coefficients as it computes
    the [gcd]. We follow the algorithm given in Knuth's
    "Art of Computer Programming", vol 2, page 325. *)

Section extended_euclid_algorithm.

Variable a b : Z.

(** The specification of Euclid's algorithm is the existence of
    [u], [v] and [d] such that [ua+vb=d] and [(gcd a b d)]. *)

Inductive Euclid : Set :=
    Euclid_intro :
      forall u v d : Z, (u * a + v * b)%Z = d -> gcd a b d -> Euclid.

(** 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)]. 
    *)

Lemma euclid_rec :
 forall v3 : Z,
 (0 <= v3)%Z ->
 forall u1 u2 u3 v1 v2 : Z,
 (u1 * a + u2 * b)%Z = u3 ->
 (v1 * a + v2 * b)%Z = v3 ->
 (forall d : Z, gcd u3 v3 d -> gcd a b d) -> Euclid.
Proof.
(* Goal: forall (v3 : Z) (_ : Z.le Z0 v3) (u1 u2 u3 v1 v2 : Z) (_ : @eq Z (Z.add (Z.mul u1 a) (Z.mul u2 b)) u3) (_ : @eq Z (Z.add (Z.mul v1 a) (Z.mul v2 b)) v3) (_ : forall (d : Z) (_ : gcd u3 v3 d), gcd a b d), Euclid *)
intros v3 Hv3; generalize Hv3; pattern v3 in |- *.
(* Goal: (fun z : Z => forall (_ : Z.le Z0 z) (u1 u2 u3 v1 v2 : Z) (_ : @eq Z (Z.add (Z.mul u1 a) (Z.mul u2 b)) u3) (_ : @eq Z (Z.add (Z.mul v1 a) (Z.mul v2 b)) z) (_ : forall (d : Z) (_ : gcd u3 z d), gcd a b d), Euclid) v3 *)
apply Z_lt_rec.
(* Goal: forall (x : Z) (_ : forall (y : Z) (_ : and (Z.le Z0 y) (Z.lt y x)) (_ : Z.le Z0 y) (u1 u2 u3 v1 v2 : Z) (_ : @eq Z (Z.add (Z.mul u1 a) (Z.mul u2 b)) u3) (_ : @eq Z (Z.add (Z.mul v1 a) (Z.mul v2 b)) y) (_ : forall (d : Z) (_ : gcd u3 y d), gcd a b d), Euclid) (_ : Z.le Z0 x) (u1 u2 u3 v1 v2 : Z) (_ : @eq Z (Z.add (Z.mul u1 a) (Z.mul u2 b)) u3) (_ : @eq Z (Z.add (Z.mul v1 a) (Z.mul v2 b)) x) (_ : forall (d : Z) (_ : gcd u3 x d), gcd a b d), Euclid *)
(* Goal: Z.le Z0 v3 *)
clear v3 Hv3; intros.
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
elim (Z_zerop x); intro.
(* Goal: Euclid *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
apply Euclid_intro with (u := u1) (v := u2) (d := u3).
(* Goal: Z.le Z0 v3 *)
assumption.
(* Goal: gcd a b u3 *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
apply H2.
(* Goal: divide d (Z.add (Z.sub a (Z.mul q b)) (Z.mul q b)) *)
(* Goal: @eq Z (Z.add (Z.sub a (Z.mul q b)) (Z.mul q b)) a *)
rewrite a0; auto.
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
set (q := (u3 / x)%Z) in *.
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
assert (Hq : (0 <= u3 - q * x < x)%Z).
(* Goal: and (Z.le Z0 (Z.sub u3 (Z.mul q x))) (Z.lt (Z.sub u3 (Z.mul q x)) x) *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
replace (u3 - q * x)%Z with (u3 mod x)%Z.
(* Goal: and (Z.le Z0 (Z.modulo u3 x)) (Z.lt (Z.modulo u3 x) x) *)
(* Goal: @eq Z (Z.modulo u3 x) (Z.sub u3 (Z.mul q x)) *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
apply Z_mod_lt; omega.
(* Goal: @eq Z (Z.modulo u3 x) (Z.sub u3 (Z.mul q x)) *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
assert (xpos : (x > 0)%Z). omega.
(* Goal: @eq Z (Z.modulo u3 x) (Z.sub u3 (Z.mul q x)) *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
generalize (Z_div_mod_eq u3 x xpos). 
(* Goal: forall _ : @eq Z u3 (Z.add (Z.mul x (Z.div u3 x)) (Z.modulo u3 x)), @eq Z (Z.modulo u3 x) (Z.sub u3 (Z.mul q x)) *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
unfold q in |- *. 
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
intro eq; pattern u3 at 2 in |- *; rewrite eq; ring.
apply
 (H (u3 - q * x)%Z Hq (proj1 Hq) v1 v2 x (u1 - q * v1)%Z (u2 - q * v2)%Z).
(* Goal: divide d (Z.add (Z.sub a (Z.mul q b)) (Z.mul q b)) *)
(* Goal: @eq Z (Z.add (Z.sub a (Z.mul q b)) (Z.mul q b)) a *)
tauto.
replace ((u1 - q * v1) * a + (u2 - q * v2) * b)%Z with
 (u1 * a + u2 * b - q * (v1 * a + v2 * b))%Z.
(* Goal: @eq Z (Z.sub (Z.add (Z.mul u1 a) (Z.mul u2 b)) (Z.mul q (Z.add (Z.mul v1 a) (Z.mul v2 b)))) (Z.sub u3 (Z.mul q x)) *)
(* Goal: @eq Z (Z.sub (Z.add (Z.mul u1 a) (Z.mul u2 b)) (Z.mul q (Z.add (Z.mul v1 a) (Z.mul v2 b)))) (Z.add (Z.mul (Z.sub u1 (Z.mul q v1)) a) (Z.mul (Z.sub u2 (Z.mul q v2)) b)) *)
(* Goal: forall (d : Z) (_ : gcd x (Z.sub u3 (Z.mul q x)) d), gcd a b d *)
(* Goal: Z.le Z0 v3 *)
rewrite H0; rewrite H1; trivial.
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
ring.
(* Goal: gcd a b u3 *)
(* Goal: Euclid *)
(* Goal: Z.le Z0 v3 *)
intros; apply H2.
(* Goal: Z.le Z0 v3 *)
apply gcd_for_euclid with q; assumption.
(* Goal: Z.le Z0 v3 *)
assumption.
Qed.

(** We get Euclid's algorithm by applying [euclid_rec] on
    [1,0,a,0,1,b] when [b>=0] and [1,0,a,0,-1,-b] when [b<0]. *)

Lemma euclid : Euclid.
Proof.
(* Goal: Euclid *)
case (Z_le_gt_dec 0 b); intro.
intros;
 apply
  euclid_rec
   with (u1 := 1%Z) (u2 := 0%Z) (u3 := a) (v1 := 0%Z) (v2 := 1%Z) (v3 := b);
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
 auto; ring.
intros;
 apply
  euclid_rec
   with
     (u1 := 1%Z)
     (u2 := 0%Z)
     (u3 := a)
     (v1 := 0%Z)
     (v2 := (-1)%Z)
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
     (v3 := (- b)%Z); auto; try ring.
(* Goal: Z.le Z0 (Z.opp b) *)
omega.
Qed.

End extended_euclid_algorithm.

Theorem gcd_uniqueness_apart_sign :
 forall a b d d' : Z, gcd a b d -> gcd a b d' -> d = d' \/ d = (- d')%Z.
Proof.
(* Goal: forall (a b d d' : Z) (_ : gcd a b d) (_ : gcd a b d'), or (@eq Z d d') (@eq Z d (Z.opp d')) *)
simple induction 1.
(* Goal: forall (_ : divide d a) (_ : divide d b) (_ : forall (x : Z) (_ : divide x a) (_ : divide x b), divide x d) (_ : gcd a b d'), or (@eq Z d d') (@eq Z d (Z.opp d')) *)
intros H1 H2 H3; simple induction 1; intros.
(* Goal: or (@eq Z d d') (@eq Z d (Z.opp d')) *)
generalize (H3 d' H4 H5); intro Hd'd.
(* Goal: or (@eq Z d d') (@eq Z d (Z.opp d')) *)
generalize (H6 d H1 H2); intro Hdd'.
(* Goal: or (@eq Z d d') (@eq Z d (Z.opp d')) *)
exact (divide_antisym d d' Hdd' Hd'd).
Qed.

(** * Bezout's coefficients *)

Inductive Bezout (a b d : Z) : Prop :=
    Bezout_intro : forall u v : Z, (u * a + v * b)%Z = d -> Bezout a b d.

(** Existence of Bezout's coefficients for the [gcd] of [a] and [b] *)

Lemma gcd_bezout : forall a b d : Z, gcd a b d -> Bezout a b d.
Proof.
(* Goal: forall (a b d : Z) (_ : gcd a b d), Bezout a b d *)
intros a b d Hgcd.
(* Goal: Bezout a b d *)
elim (euclid a b); intros.
(* Goal: Bezout a b d *)
generalize (gcd_uniqueness_apart_sign a b d d0 Hgcd g).
(* Goal: forall _ : or (@eq Z d d0) (@eq Z d (Z.opp d0)), Bezout a b d *)
intro H; elim H; clear H; intros.
(* Goal: Bezout a b d *)
(* Goal: Bezout a b d *)
apply Bezout_intro with u v.
(* Goal: Z.le Z0 v3 *)
rewrite H; assumption.
(* Goal: Bezout a b d *)
apply Bezout_intro with (- u)%Z (- v)%Z.
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
rewrite H; rewrite <- e; ring.
Qed.

(** gcd of [ca] and [cb] is [c gcd(a,b)]. *)

Lemma gcd_mult : forall a b c d : Z, gcd a b d -> gcd (c * a) (c * b) (c * d).
Proof.
(* Goal: forall (a b c d : Z) (_ : gcd a b d), gcd (Z.mul c a) (Z.mul c b) (Z.mul c d) *)
simple induction 1; constructor; intuition.
(* Goal: divide x (Z.mul c d) *)
elim (gcd_bezout a b d H); intros.
(* Goal: divide x (Z.mul c d) *)
elim H3; intros.
(* Goal: divide x (Z.mul c d) *)
elim H4; intros.
(* Goal: divide x (Z.mul c d) *)
apply divide_intro with (u * q + v * q0)%Z.
(* Goal: @eq Z (Z.mul c d) (Z.mul (Z.add (Z.mul u q) (Z.mul v q0)) x) *)
rewrite <- H5.
(* Goal: @eq Z (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) (Z.mul (Z.add (Z.mul u q) (Z.mul v q0)) x) *)
replace (c * (u * a + v * b))%Z with (u * (c * a) + v * (c * b))%Z.
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
rewrite H6; rewrite H7; ring.
(* Goal: @eq Z (Z.add (Z.mul u (Z.mul c a)) (Z.mul v (Z.mul c b))) (Z.mul c (Z.add (Z.mul u a) (Z.mul v b))) *)
ring.
Qed.
