(** First, we show the following theorem: if p is a prime number and
gcd(p,k)=1 then sqrt(p*k) is not rational *)

(** Then, we strengthen the result to the n-th root of p^r*k where 0 <
r < n obtaining the theorem: if p is a prime number, gcd(p,k)=1 and 0
< r < n then the n-th root of p^r*k is not rational *)

(** To achieve this goal, we define the notion of divisibility in N,
then we show the (constructive) existence and uniqueness of the gcd
(still in N) We show also Bezout and Gauss's theorems. We have
volontarily made the choice to not use existing definitions of the
standard library that work in Z *) 

(** author: sebastien.briais@epfl.ch *)
Require Import Arith.
Require Import ArithRing.
Require Import Omega.

(** We first begin with some lemmas that relates +, * and - that are not in the standard library *)
Lemma mult_lemma1 : forall (n m:nat),(n <> O)->(m <> 0)->(n <= n*m).
intros.
rewrite mult_comm.
induction m;simpl;auto with arith.
elim H0;trivial.
Qed.

Lemma mult_lemma2 : forall (n m:nat),(n*m = O)->(n=O)\/(m=O).
intros.
induction n.
tauto.
simpl in H.
right.
assert (m <= O);try omega.
rewrite <- H.
auto with arith.
Qed.

Lemma mult_lemma3 : forall (n m:nat),(n <> O)->(m > 1)->(n < n*m).
intros.
rewrite mult_comm.
induction m.
inversion H0.
simpl.
assert (O < m*n);try omega.
inversion H0;try omega.
assert (1 <= n);try omega.
assert (m > 1);try omega.
generalize (IHm H4);omega.
Qed.

Lemma mult_lemma4 : forall (n m:nat),n=n*m -> n=O \/ m=1.
intros n m.
case n.
left;trivial.
intros.
right.
destruct m.
rewrite mult_comm in H.
discriminate.
destruct m;trivial.
assert ((S n0)<(S n0)*(S (S m))).
apply mult_lemma3;intros;auto with arith.
rewrite <- H in H0.
elim (lt_irrefl (S n0) H0).
Qed.

Lemma mult_lemma5 : forall (n m:nat),((n * m) =1)->(n=1)/\(m=1).
induction n;simpl;intros;try discriminate.
induction m.
rewrite mult_comm in H.
simpl in H;discriminate.
assert ((S n)<=((S n)*(S m))).
apply mult_lemma1;discriminate.
assert (((S n)*(S m))=((S m)+n*(S m))).
reflexivity.
rewrite H1 in H0.
rewrite H in H0.
assert ((S n)=1).
omega.
split;trivial.
inversion H2.
rewrite H4 in H.
simpl in H.
omega.
Qed.

Lemma plus_minus_lemma1 : forall (y x:nat),(x+y-y=x).
induction y;intros;rewrite plus_comm;simpl.
auto with arith.
rewrite plus_comm.
apply IHy.
Qed.

Lemma mult_minus_lemma1 : forall (a n:nat),a*n-n = (a-1)*n.
intros.
induction a.
simpl.
trivial.
replace (S a*n) with (n+a*n);try (auto with arith).
rewrite plus_comm.
rewrite plus_minus_lemma1.
simpl.
rewrite <- minus_n_O;trivial.
Qed.

Lemma mult_lemma6 : forall (a b n:nat),(n <> O)->(n*a=n*b)->(a=b).
induction a.
intros;rewrite <- mult_n_O in H0; generalize (mult_lemma2 n b); intros Hl2; elim Hl2; intros; (auto || elim H ; auto).
intros b n H.
rewrite mult_comm;simpl;rewrite mult_comm;intro.
assert (n*a = n*b-n).
apply plus_minus;auto.
assert (a*n=(b-1)*n).
rewrite <- mult_minus_lemma1;rewrite mult_comm;rewrite (mult_comm b n);trivial.
assert (a=(b-1)).
apply (IHa (b-1) n);trivial.
rewrite mult_comm;rewrite (mult_comm n (b-1));trivial.
destruct b;simpl in H3.
rewrite H3 in H0;rewrite (mult_comm n 0) in H0;rewrite plus_comm in H0;simpl in H0;elim H;trivial.
rewrite <- minus_n_O in H3;auto.
Qed.

Lemma mult_lemma7 : forall (x y z t:nat),x*y*(z*t)=z*(x*y*t).
intros.
ring.
Qed.

Lemma minus_lemma1 : forall (a b:nat),(S a-S b)<S a.
intros.
omega.
Qed.

Lemma minus_lemma2 : forall (n m:nat),(n<=m)->(n-m=O).
intros.
omega.
Qed.

Lemma mult_minus_lemma2 : forall (x y z:nat),(x*(y-z))=(x*y-x*z).
intros.
case (le_lt_dec y z);intro.
rewrite (minus_lemma2 y z l);rewrite mult_comm;simpl;rewrite minus_lemma2;trivial;auto with arith.
assert (y=z+(y-z)).
rewrite <- (le_plus_minus z y);try (auto with arith).
replace (x*y) with (x*(z+(y-z))).
rewrite mult_plus_distr_l;rewrite minus_plus;trivial.
rewrite <- H;trivial.
Qed.

Lemma plus_minus_lemma2 : forall (x y z:nat),(y<=x)->(x-y+z)=(x+z-y).
intros.
rewrite (le_plus_minus y x);try (auto with arith).
rewrite minus_plus;rewrite <- plus_assoc;rewrite minus_plus;trivial.
Qed.

Lemma minus_minus_lemma1 : forall (x y z:nat),(z<=y)->(x-(y-z))=(x+z-y).
intros.
rewrite (le_plus_minus z y);trivial.
rewrite minus_plus;rewrite plus_comm;rewrite <- minus_plus_simpl_l_reverse;trivial.
Qed.

Lemma minus_minus_lemma2 : forall (x y z:nat),(x-y-z)=(x-(y+z)).
induction x;simpl;trivial.
intros.
case y;simpl;trivial.
Qed.

Lemma minus_lt_lemma1 : forall (b a:nat),(a<b)->(0<b-a).
intros.
omega.
Qed.

(** b | a if there is q such that a = b * q*)
Definition divides (a b:nat) := exists q:nat,a = (b*q).

(** 1 divides every natural number *)
Lemma one_min_div : forall (n:nat),(divides n 1).
intros.
red.
exists n.
auto with arith.
Qed.

(** 0 is divides by every natural number *)
Lemma zero_max_div : forall (n:nat),(divides O n).
intros.
red.
exists O.
auto with arith.
Qed.

(** the relation of divisibility is reflexive *)
Lemma divides_refl : forall (a:nat),(divides a a).
intros.
red.
exists 1.
auto with arith.
Qed.

(** the relation of divisibility is transitive *)
Lemma divides_trans : forall (a b c:nat),(divides a b)->(divides b c)->(divides a c).
unfold divides.
intros.
elim H;intro q;intro.
elim H0;intro q';intro.
rewrite H2 in H1.
exists (q' * q).
rewrite H1.
auto with arith.
Qed.

(** the relation of divisibility is antisymmetric *)
Lemma divides_antisym : forall (a b:nat),(divides a b)->(divides b a)->a=b.
unfold divides.
intros.
elim H;intro q;intro.
elim H0;intro q';intro.
rewrite H2 in H1.
assert ((a = 0) \/ (q' * q)=1).
apply mult_lemma4.
replace (a*(q'*q)) with (a*q'*q);try (auto with arith).
case H3;intro.
rewrite H4 in H2;simpl in H2;rewrite H2;trivial.
elim (mult_lemma5 q' q H4);intros.
rewrite H5 in H2;rewrite mult_comm in H2;simpl in H2;rewrite plus_comm in H2;simpl in H2;symmetry;trivial.
Qed.

(** corollary: forall a<>1, not(a | 1) *)
Lemma non_div_1 : forall (a:nat),(a<>1)->~(divides 1 a).
intros.
red.
intro.
apply H.
apply divides_antisym;trivial.
apply one_min_div.
Qed.

(** if d | a and d | b then d | (a+b) *)
Lemma divides_plus : forall (d a b:nat),(divides a d)->(divides b d)->(divides (plus a b) d).
unfold divides.
intros.
elim H;intro q;intro.
elim H0;intro q';intro.
exists (q+q').
rewrite H1;rewrite H2.
ring.
Qed.

(** if d | a then d | a*b *)
Lemma divides_mult : forall (d a b:nat),(divides a d)->(divides (a*b) d).
unfold divides.
intros.
elim H;intro q;intro.
exists (b * q).
rewrite H0.
ring.
Qed.

(** if d | a and d | b then d | (b-a) *)
Lemma divides_minus : forall (d a b:nat),(divides a d)->(divides b d)->(divides (b-a) d).
unfold divides.
intros.
elim H;intro q;intro.
elim H0;intro q';intro.
rewrite H1;rewrite H2.
exists (q'-q).
rewrite (mult_comm d q');rewrite (mult_comm d q);rewrite (mult_comm d (q'-q));auto with arith.
Qed.

(** d is a common divisor of a and b if d | a and d | b *)
Definition is_cd (d a b : nat) := (divides a d)/\(divides b d).

(** d is a greatest common divisor of a and b if it is a maximal common divisor *)
Definition is_gcd (d a b:nat) := (is_cd d a b)/\(forall (d':nat),(is_cd d' a b)->(divides d d')).

(** there is at most one gcd of a and b *)
Theorem gcd_unique : forall (d d' a b:nat),(is_gcd d a b)->(is_gcd d' a b)->d=d'.
unfold is_gcd.
intros.
elim H;elim H0;intros.
apply divides_antisym;auto.
Qed.

(** gcd(a,b) = gcd(b,a) *)
Lemma gcd_sym : forall (d a b:nat),(is_gcd d a b)->(is_gcd d b a).
unfold is_gcd.
intros.
elim H;intros.
split.
red;red in H0;tauto.
intros.
apply H1.
red;red in H2;tauto.
Qed.

(** gcd(0,a)=a *)
Lemma gcd_zero : forall (a:nat),(is_gcd a O a).
unfold is_gcd.
intro.
split.
red;split;[apply zero_max_div | apply divides_refl].
unfold is_cd;tauto.
Qed.

(** gcd(1,a)=1 *)
Lemma gcd_one : forall (a:nat),(is_gcd 1 1 a).
unfold is_gcd.
intros.
split.
red;split;[apply divides_refl | apply one_min_div].
unfold is_cd;tauto.
Qed.

(** if a <= b then gcd(a,b)=gcd(a,b-a) *)
Lemma gcd_minus : forall (d a b:nat),(a <= b)->((is_gcd d a b)<->(is_gcd d a (b-a))).
intros.
unfold is_gcd.
split;intro.
elim H0;intros.
split.
red in H1;red.
elim H1;intros.
split;try tauto.
apply divides_minus;trivial.
unfold is_cd;intros.
apply H2;red;elim H3;intros.
split;[tauto | rewrite (le_plus_minus a b H);apply divides_plus;trivial].
elim H0;unfold is_cd;intros.
split.
split;[tauto | elim H1;intros;rewrite (le_plus_minus a b H);apply divides_plus;trivial].
intros.
elim H3;intros;apply H2.
split;try (apply divides_minus);trivial.
Qed.

(** gcd(a,a) = a *)
Lemma gcd_refl : forall (a:nat),(is_gcd a a a).
unfold is_gcd.
intros.
unfold is_cd.
split;try tauto.
split;apply divides_refl.
Qed.

(** two trivial lemmas: gcd(a,b) | a *)
Lemma gcd_div_l : forall (d a b:nat),(is_gcd d a b)->(divides a d).
unfold is_gcd;unfold is_cd;intros;tauto.
Qed.

(** gcd(a,b) | b *)
Lemma gcd_div_r : forall (d a b:nat),(is_gcd d a b)->(divides b d).
unfold is_gcd;unfold is_cd;intros;tauto.
Qed.

(** we now show that gcd(a,b) exists for all a and b (we even have an algorithm) *)
Require Import Wf_nat.

Definition f (x:nat*nat) := (fst x)+(snd x).

Definition R (x y:nat*nat) := (f x)<(f y).

Lemma Rwf : well_founded R.
unfold R.
apply (well_founded_ltof (nat*nat) f).
Qed.

(** proof of existence of gcd(a,b): it relies on the relation gcd(a,b)=gcd(a,b-a) if a<=b *)
Lemma gcd_exists_prod : forall (x:nat*nat),{d:nat | (is_gcd d (fst x) (snd x))}.
apply (induction_ltof2 (nat*nat) f (fun x:nat*nat => {d:nat | (is_gcd d (fst x) (snd x))})).
unfold ltof.
unfold f.
intros.
case (lt_eq_lt_dec (fst x) (snd x));intro.
case s;intro.
destruct (fst x).
exists (snd x);apply gcd_zero.
elim (H (S n,snd x-S n)).
simpl;intro d;intro.
exists d.
elim (gcd_minus d (S n) (snd x));try (auto with arith).
simpl.
omega.
rewrite e;exists (snd x);apply gcd_refl.
destruct (snd x).
exists (fst x);apply gcd_sym;apply gcd_zero.
elim (H (S n,fst x-S n)).
simpl;intro d;intro.
exists d.
apply gcd_sym.
elim (gcd_minus d (S n) (fst x));try (auto with arith).
simpl.
omega.
Qed.

(** here we are: the gcd exists *)
Theorem gcd_exists : forall (a b:nat),{d:nat | (is_gcd d a b)}.
intros.
elim (gcd_exists_prod (a,b)).
simpl;intro d;intro;exists d;trivial.
Qed.

(** take the first projection of the proof *)
Definition gcd (a b:nat) := let (d,_):=(gcd_exists a b) in d.

(** the gcd is the gcd! *)
Lemma gcd_is_gcd : forall (a b:nat),(is_gcd (gcd a b) a b).
intros.
unfold gcd.
generalize (gcd_exists a b).
intro;elim s;intro d;trivial.
Qed.

(** a and b are relatively prime if gcd(a,b)=1 *)
Definition rel_prime (a b:nat) := (is_gcd 1 a b). 

(** if a and b are relatively prime then so are b and a *)
Lemma rel_prime_sym : forall (a b:nat),(rel_prime a b)->(rel_prime b a).
unfold rel_prime.
intros;apply gcd_sym;trivial.
Qed.

(** for all a, a and 1 are relatively prime *)
Lemma rel_prime_1 : forall (a:nat),(rel_prime a 1).
unfold rel_prime.
intros;apply gcd_sym;apply gcd_one.
Qed.

(** here we show that if b | a then it is possible to compute q such that a = b*q *)
Lemma quo_dec : forall (a b:nat),(divides a b)->{q:nat | a=b*q}.
intros.
apply (lt_wf_rec a (fun x:nat => (divides x b)->{q:nat | x = b*q}));trivial.
intro.
case n;intros.
exists 0;auto with arith.
elim (H0 ((S n0)-b)).
intro q;intro.
exists (S q).
replace (S n0) with (b+(S n0-b)).
rewrite p;rewrite plus_comm;auto with arith.
symmetry.
apply le_plus_minus.
elim H1;intros.
rewrite H2.
replace (b <= b*x) with (1*b <= b*x);rewrite (mult_comm b x).
apply mult_le_compat_r.
destruct x;[rewrite mult_comm in H2;discriminate | auto with arith].
simpl;auto with arith.
destruct b.
elim H1;simpl;intros;discriminate.
omega.
apply (divides_minus b b (S n0));[apply divides_refl | trivial].
Qed.

(** we can now define the quotient of a by b in case of b | a *)
Definition quo (a b:nat) (H:(divides a b)) := let (q,_):=(quo_dec a b H) in q.

(** the quotient is the quotient! *)
Lemma quo_is_quo : forall (a b:nat)(H:divides a b),a=(mult b (quo a b H)).
intros.
unfold quo.
generalize (quo_dec a b H);intro;elim s;trivial.
Qed.

(** if b | a then (n*a/b) = n*(a/b) *) 
Lemma quo_mult : forall (a b:nat)(H:divides a b),forall (n:nat),(b<>O)->(quo (a*n) b (divides_mult b a n H))=n*(quo a b H).
intros.
generalize (quo_is_quo (a*n) b (divides_mult b a n H));intro.
generalize (quo_is_quo a b H);intro.
replace (a*n = b * quo (a * n) b (divides_mult b a n H)) with (b*(quo a b H)*n = b * quo (a * n) b (divides_mult b a n H)) in H1.
symmetry;rewrite mult_comm.
apply mult_lemma6 with b;trivial.
rewrite mult_assoc;trivial.
rewrite <- H2;trivial.
Qed.

(** we have that a/gcd(a,b) and b/gcd(a,b) are relatively prime *)
Lemma gcd_rel_prime : forall (d a b:nat)(H:(is_gcd d a b)),(d <> O)->(rel_prime (quo a d (gcd_div_l d a b H)) (quo b d (gcd_div_r d a b H))).
unfold rel_prime.
intros.
generalize (quo_is_quo a d (gcd_div_l d a b H));intro.
generalize (quo_is_quo b d (gcd_div_r d a b H));intro.
unfold is_gcd;split;unfold is_cd.
split;apply one_min_div.
intros.
elim H3;intros.
elim H4;intro q;intro.
elim H5;intro q';intro.
rewrite H6 in H1.
rewrite H7 in H2.
assert (divides d (d*d')).
red in H;elim H;intros.
apply H9;red;split;[exists q;rewrite H1;ring | exists q';rewrite H2;ring].
elim H8;intros.
exists x.
apply mult_lemma6 with d;trivial.
rewrite mult_assoc;rewrite <- H9;auto with arith.
Qed.

(** definition of square *)
Definition square (x:nat) := x*x.

(** (x*y)^2 = x^2*y^2 *)
Lemma square_mult_lemma : forall (a b:nat),(square (a*b))=((square a)*(square b)).
unfold square.
intros.
ring.
Qed.

(** if q<>0 then gcd(p,q)<>0 *)
Lemma gcd_non_zero : forall (d p q:nat),(q<>O)->(is_gcd d p q)->(d<>O).
unfold is_gcd.
intros.
elim H0;intros.
intro.
elim H1;intros.
elim H5;intros.
rewrite H3 in H6;simpl in H6;auto.
Qed.

(** we now exhibit an algorithm that computes Bezout coefficient: for all a b, there is u and v such that a*u-b*v = gcd(a,b) or b*v-a*u = gcd(a,b) *)
(** the 4 lemmae gives the idea of the algorithm *)
Lemma bezout_aux1 : forall (x y u v:nat),(x<=y)->(is_gcd (u*x-v*(y-x)) x (y-x))->(is_gcd ((u+v)*x-v*y) x y).
intros.
elim (gcd_minus ((u+v)*x-v*y) x y H);intros.
apply H2.
rewrite mult_plus_distr_r;rewrite <- minus_minus_lemma1;try (auto with arith);rewrite <- mult_minus_lemma2;trivial.
Qed.

Lemma bezout_aux2 : forall (x y u v:nat),(x<=y)->(is_gcd (v*(y-x)-u*x) x (y-x))->(is_gcd (v*y-(u+v)*x) x y).
intros.
elim (gcd_minus (v*y-(u+v)*x) x y H);intros.
apply H2.
rewrite mult_plus_distr_r;rewrite plus_comm;rewrite <- minus_minus_lemma2;rewrite <- mult_minus_lemma2;trivial.
Qed.

(** Bezout coefficient *)
Lemma bezout_exists_prod : forall (x:nat*nat),{y:nat*nat | (is_gcd ((fst y)*(fst x)-(snd y)*(snd x)) (fst x) (snd x))}+{y:nat*nat | (is_gcd ((snd y)*(snd x)-(fst y)*(fst x)) (fst x) (snd x))}.
apply (induction_ltof2 (nat*nat) f (fun x:nat*nat => ({y:nat*nat | (is_gcd ((fst y)*(fst x)-(snd y)*(snd x)) (fst x) (snd x))}+{y:nat*nat | (is_gcd ((snd y)*(snd x)-(fst y)*(fst x)) (fst x) (snd x))})%type)).
unfold ltof.
unfold f.
intros.
case (lt_eq_lt_dec (fst x) (snd x));intro.
case s;intro.
destruct (fst x).
right;exists (0,1);simpl;rewrite <- minus_n_O;rewrite plus_comm;simpl;apply gcd_zero.
elim (H (S n,snd x-S n));try (intro;simpl).
elim a;intro y;intro.
left;exists ((fst y)+(snd y),(snd y)).
simpl;apply bezout_aux1;try (auto with arith).
elim b;intro y;intro.
right;exists ((fst y)+(snd y),(snd y)).
simpl;apply bezout_aux2;try (auto with arith).
simpl;omega.
rewrite e;left;exists (1,0);simpl;rewrite <- minus_n_O;rewrite plus_comm;simpl;apply gcd_refl.
destruct (snd x).
left;exists (1,0);simpl;rewrite <- minus_n_O;rewrite plus_comm;simpl;apply gcd_sym;apply gcd_zero.
elim (H (S n,fst x-S n));try (intro;simpl).
elim a;intro y;intro.
right;exists ((snd y),(fst y)+(snd y));apply gcd_sym.
simpl;apply bezout_aux1;try (auto with arith).
elim b;intro y;intro.
left;exists ((snd y),(fst y)+(snd y));apply gcd_sym.
simpl;apply bezout_aux2;try (auto with arith).
simpl;omega.
Qed.

(** Bezout' theorem *)
Theorem bezout_exists : forall (a b:nat),{u:nat & {v:nat | (is_gcd (a*u-b*v) a b)}}+{u:nat & {v:nat | (is_gcd (b*v-a*u) a b)}}.
intros.
elim (bezout_exists_prod (a,b));intro.
elim a0;destruct x;simpl;intros.
left;exists n;exists n0;rewrite mult_comm;rewrite (mult_comm b);trivial.
elim b0;destruct x;simpl;intros.
right;exists n;exists n0;rewrite mult_comm;rewrite (mult_comm a);trivial.
Qed.

(** Bezout' theorem reformulated *)
Theorem bezout : forall (d a b:nat),(is_gcd d a b)->exists u:nat,exists v:nat,d=a*u-b*v \/ d=b*v-a*u.
intros.
elim (bezout_exists a b);intro.
elim a0;intro u;intro;elim p;intro v;intro;exists u;exists v;left;apply (gcd_unique d (a*u-b*v) a b);trivial.
elim b0;intro u;intro;elim p;intro v;intro;exists u;exists v;right;apply (gcd_unique d (b*v-a*u) a b);trivial.
Qed.

(** Bezout' theorem and relatively prime numbers *)
Theorem bezout_rel_prime : forall (a b:nat),(rel_prime a b)<->(exists u:nat, exists v:nat, 1=a*u-b*v \/ 1 = b*v-a*u).
intros.
unfold rel_prime.
split;intro.
apply bezout;trivial.
elim H;intro u;intro H0.
elim H0;intro v;intro.
unfold is_gcd;unfold is_cd.
split.
split;apply one_min_div.
intros.
elim H2;intros.
elim H3;intro q;intro.
elim H4;intro q';intro.
rewrite H5 in H1;rewrite H6 in H1.
case H1;intro.
exists (q*u-q'*v);rewrite mult_minus_lemma2;rewrite mult_assoc;rewrite mult_assoc;trivial.
exists (q'*v-q*u);rewrite mult_minus_lemma2;rewrite mult_assoc;rewrite mult_assoc;trivial.
Qed.

(** gcd(n*a,n*b) = n*gcd(a,b) *)
Lemma gcd_mult : forall (d a b:nat),(is_gcd d a b)->(forall (n:nat),(is_gcd (n*d) (n*a) (n*b))).
unfold is_gcd;unfold is_cd.
intros.
elim H;intros.
elim H0;intros.
split.
elim H2;intro q;intro.
elim H3;intro q';intro.
rewrite H4;rewrite mult_assoc.
rewrite H5;rewrite mult_assoc.
split;[exists q;trivial | exists q';trivial].
intros.
elim H4;intros.
elim (bezout d a b);try (unfold is_gcd;unfold is_cd;trivial).
intro u;intro.
elim H7;intro v;intro.
elim H5;intro q;intro.
elim H6;intro q';intro.
case H8;intro;[exists (q*u-q'*v) | exists (q'*v-q*u)];rewrite mult_minus_lemma2;rewrite mult_assoc;rewrite mult_assoc;rewrite <- H9;rewrite <- H10;rewrite H11;rewrite mult_minus_lemma2;rewrite mult_assoc;rewrite mult_assoc;trivial.
Qed.

(** Gauss' theorem (use Bezout) *)
Theorem gauss : forall (d a b:nat),(rel_prime a d)->(divides (a*b) d)->(divides b d).
unfold rel_prime.
intros.
elim (bezout 1 a d H);intro u;intro.
elim H1;intro v;intro.
elim H0;intro q;intro.
case H2;intro;[exists (q*u-b*v) | exists (b*v-q*u)];rewrite mult_minus_lemma2;rewrite mult_assoc;rewrite mult_assoc;rewrite <- H3;rewrite (mult_comm a b);rewrite (mult_comm d b);rewrite <- mult_assoc;rewrite <- mult_assoc;rewrite <- mult_minus_lemma2;rewrite <- H4;auto with arith.
Qed.

(** definition of a "prime number" *)
Definition is_prime (p:nat) := (p<>1)/\(forall (d:nat),(divides p d)->(d=1)\/(d=p)).

Lemma not_prime_zero : ~(is_prime O).
unfold is_prime.
intro.
elim H;intros.
case (H1 2);try (apply zero_max_div);intro;discriminate.
Qed.

(** if p is prime and not(p | a) then a and p are relatively prime *)
Lemma prime_div_gcd : forall (p a:nat),(is_prime p)->~(divides a p)->(rel_prime p a).
unfold is_prime.
unfold rel_prime.
intros.
unfold is_gcd;unfold is_cd.
split.
split;apply one_min_div.
intros.
elim H;intros.
elim H1;intros.
case (H3 d' H4);intro.
rewrite H6;apply divides_refl.
rewrite H6 in H5;tauto.
Qed.

(** We show some lemmae about prime numbers *)
(** if p is prime then gcd(a,p)=1 or gcd(a,p)=p *)
Lemma prime_gcd : forall (d p a:nat),(is_prime p)->(is_gcd d a p)->(d=1)\/(d=p).
unfold is_prime.
intros.
elim H;intros.
apply H2.
elim H0;intros.
elim H3;trivial.
Qed.

(** if p is prime and gcd(a,p)<>1 then p | a *)
Lemma prime_rel_prime : forall (p a:nat),(is_prime p)->~(rel_prime a p)->(divides a p).
unfold rel_prime.
intros.
generalize (gcd_is_gcd a p);intros.
case (prime_gcd (gcd a p) p a H H1);intro;rewrite H2 in H1;try tauto.
elim H1;intros.
elim H3;trivial.
Qed.

(** if p is prime and p | a*b then p | a or p | b *)
Lemma prime_mult : forall (p a b:nat),(is_prime p)->(divides (a*b) p)->((divides a p)\/(divides b p)).
intros.
generalize (gcd_is_gcd a p);intro.
case (prime_gcd (gcd a p) p a H H1);intro;rewrite H2 in H1.
right;apply gauss with a;trivial.
red in H1;elim H1;intros.
red in H3;tauto.
Qed.

(** corollary: if p is prime and p | a^2 then p | a *)
Lemma prime_square : forall (p a:nat),(is_prime p)->(divides (square a) p)->(divides a p).
unfold square.
intros;case (prime_mult p a a H H0);trivial.
Qed.

(** now, we show the result claimed in the header *)
Lemma sqrt_prime_irrat_aux : forall (p k a b:nat),(is_prime p)->(rel_prime p k)->(rel_prime a b)->(p*k*(square b) <> (square a)).
intros.
intro.
assert (divides a p).
apply prime_square;trivial.
exists (k*(square b)).
rewrite <- H2;ring.
elim H3;intro n_a;intro.
rewrite H4 in H2;rewrite square_mult_lemma in H2;unfold square in H2.
assert (k*(b*b)=p*(n_a*n_a)).
apply mult_lemma6 with p.
intro H5;rewrite H5 in H;apply not_prime_zero;trivial.
rewrite mult_assoc;rewrite H2;ring.
assert (divides b p).
apply prime_square;trivial;unfold square.
apply gauss with k.
apply rel_prime_sym;trivial.
exists (n_a*n_a);trivial.
assert (p=1).
unfold rel_prime in H1.
elim H1;intros.
apply divides_antisym;try (apply one_min_div).
apply H8;red;tauto.
elim H;tauto.
Qed.

(** Theorem: if p is prime, p and k are relatively prime, then sqrt(p*k) is not rationnal *)
Theorem sqrt_prime_irrat : forall (p k a b:nat),(is_prime p)->(rel_prime p k)->(b<>O)->(p*k*(square b) <> (square a)).
intros.
generalize (gcd_is_gcd a b);intro.
generalize (quo_is_quo a (gcd a b) (gcd_div_l (gcd a b) a b H2));intro.
generalize (quo_is_quo b (gcd a b) (gcd_div_r (gcd a b) a b H2));intro.
intro.
rewrite H3 in H5.
replace (square b) with (square (gcd a b * quo b (gcd a b) (gcd_div_r (gcd a b) a b H2))) in H5;auto.
rewrite square_mult_lemma in H5;rewrite square_mult_lemma in H5.
assert (p*k*(square (quo b (gcd a b) (gcd_div_r (gcd a b) a b H2)))=(square (quo a (gcd a b) (gcd_div_l (gcd a b) a b H2)))).
apply mult_lemma6 with (square (gcd a b)).
unfold square.
generalize (gcd_non_zero (gcd a b) a b H1 H2);intro.
intro;apply H6.
case (mult_lemma2 (gcd a b) (gcd a b) H7);trivial.
rewrite <- H5;ring.
apply (sqrt_prime_irrat_aux p k (quo a (gcd a b) (gcd_div_l (gcd a b) a b H2)) (quo b (gcd a b) (gcd_div_r (gcd a b) a b H2)));auto.
apply gcd_rel_prime;apply (gcd_non_zero (gcd a b) a b);trivial.
Qed.

(** if p is prime then sqrt(p) is not rationnal *)
Fact sqrt_prime : forall (p:nat),(is_prime p)->forall (a b:nat),(b<>O)->(p*(square b)<>(square a)).
intros.
replace p with (p*1);try (auto with arith).
apply sqrt_prime_irrat;trivial;apply rel_prime_1.
Qed.

(** We now deduce from this theorem that sqrt(2) is not rationnal *)
Lemma is_prime_2 : (is_prime 2).
unfold is_prime.
split.
intro;discriminate.
intros.
elim H;destruct x;rewrite mult_comm.
intro;discriminate.
simpl.
case d.
simpl.
rewrite mult_comm;simpl;intro;discriminate.
intros.
inversion H0.
symmetry in H2.
case (plus_is_one n (x*(S n)) H2);intro.
elim a;intros.
left;rewrite H1;trivial.
elim a;intros.
right;rewrite H1;trivial.
Qed.

(** here is it! *)
Fact sqrt_2_irrat : forall (p q:nat),(q<>O)->(2*(square q)<>(square p)).
intros.
apply sqrt_prime;trivial.
apply is_prime_2.
Qed.

(*****************************************************************************)

(** we now generalize the theorem to the nth-root *)
Fixpoint power (x n:nat) {struct n} : nat :=
match n with
O => 1
| (S n) => (x*(power x n))
end.

Lemma power_mult_lemma1 : forall (n x y:nat),(power (x*y) n)=(power x n)*(power y n).
induction n;simpl;trivial.
intros;rewrite (IHn x y);ring.
Qed.

Lemma power_plus_lemma1 : forall (n m x:nat),(power x (n+m))=(power x n)*(power x m).
induction n;simpl;intros.
auto with arith.
rewrite IHn;ring.
Qed.

Lemma power_divides_lemma1 : forall (n x:nat),(0<n)->(divides (power x n) x).
induction n;simpl;intros.
inversion H.
exists (power x n);trivial.
Qed.

Lemma power_power_lemma1 : forall (n m x:nat),(power (power x n) m)=(power x (n*m)).
induction n;simpl;intros.
induction m;simpl;auto with arith.
rewrite IHm;ring.
rewrite power_mult_lemma1;rewrite IHn;rewrite <- power_plus_lemma1;trivial.
Qed.

Lemma power_zero : forall (n x:nat),(power x n)=O->x=O.
induction n;simpl;intros.
discriminate.
case (mult_lemma2 x (power x n) H);auto.
Qed.

Lemma prime_power : forall (p n x:nat),(is_prime p)->(divides (power x n) p)->(divides x p).
induction n;simpl;intros.
elim H;intros.
elim H1;apply divides_antisym;trivial;apply one_min_div.
case (prime_mult p x (power x n) H H0);trivial.
intro;apply IHn;trivial.
Qed.

Lemma prime_power_qn : forall (p n q x:nat),(is_prime p)->(divides (power x n) (power p (q*n)))->(1<=n)->(divides x (power p q)).
induction q;simpl;intros.
apply one_min_div.
rewrite power_plus_lemma1 in H0.
assert (divides x (power p q)).
apply IHq;trivial.
elim H0;intros;exists ((power p n)*x0).
rewrite H2;ring.
elim H2;intros.
rewrite H3 in H0.
rewrite power_mult_lemma1 in H0;rewrite power_power_lemma1 in H0;rewrite (mult_comm (power p n)) in H0.
elim H0;intros.
assert ((power p (q*n))<>0).
intro.
generalize (power_zero (q*n) p H5).
intro.
apply not_prime_zero.
rewrite H6 in H;trivial.
rewrite <- mult_assoc in H4.
generalize (mult_lemma6 (power x0 n) ((power p n)*x1) (power p (q*n)) H5 H4).
intro.
assert (exists n':nat,n=(S n')).
inversion H1;[exists 0 | exists m];trivial.
elim H7;intro n';intro.
rewrite H8 in H6;simpl in H6.
assert (divides x0 p).
case (prime_mult p x0 (power x0 n'));trivial.
rewrite H6.
exists ((power p n')*x1);ring.
intro.
apply prime_power with n';trivial.
elim H9;intros.
rewrite H10 in H3.
rewrite H3.
exists x2;ring.
Qed.

Lemma nth_root_irrat_aux : forall (p k a b n r:nat),(is_prime p)->(rel_prime p k)->(0<r)->(r<n)->(rel_prime a b)->((power p r)*k*(power b n) <> (power a n)).
intros.
intro.
assert (divides a p).
apply prime_power with n;trivial.
generalize (power_divides_lemma1 r p H1);intro.
elim H5;intro q;intros.
rewrite H6 in H4.
rewrite <- H4;exists (q*k*(power b n));ring.
assert (divides b p).
elim H5;intro q;intros.
rewrite H6 in H4.
rewrite power_mult_lemma1 in H4.
assert ((power p n)=(power p (r+(n-r)))).
rewrite <- le_plus_minus;try (auto with arith).
rewrite H7 in H4;rewrite power_plus_lemma1 in H4.
assert ((power p r)<>O).
intro.
apply not_prime_zero.
assert (p=O).
apply power_zero with r;trivial.
rewrite H9 in H;trivial.
rewrite <- mult_assoc in H4;rewrite <- mult_assoc in H4;generalize (mult_lemma6 (k*(power b n)) ((power p (n-r))*(power q n)) (power p r) H8 H4);intro.
assert (divides (power p (n-r)) p).
apply power_divides_lemma1;apply minus_lt_lemma1;trivial.
apply prime_power with n;trivial.
apply gauss with k;try (apply rel_prime_sym;trivial).
rewrite H9;apply divides_mult;trivial.
elim H3;intros.
elim H;intros.
apply H9;apply divides_antisym;try (apply one_min_div).
apply H8;red;tauto.
Qed.

(** generalization of the theorem: if p is a prime number, 0 < r < n and gcd(p,k)=1 then the n-th root of p^r*k is not rationnal! *)
Theorem nth_root_irrat : forall (p k a b n r:nat),(is_prime p)->(rel_prime p k)->(0<r)->(r<n)->(b<>0)->((power p r)*k*(power b n) <> (power a n)).
intros.
intro.
generalize (gcd_is_gcd a b);intro.
generalize (quo_is_quo a (gcd a b) (gcd_div_l (gcd a b) a b H5));intro.
generalize (quo_is_quo b (gcd a b) (gcd_div_r (gcd a b) a b H5));intro.
assert ((power a n)=(power (gcd a b * quo a (gcd a b) (gcd_div_l (gcd a b) a b H5)) n));try (rewrite <- H6;trivial).
assert ((power b n)=(power (gcd a b * quo b (gcd a b) (gcd_div_r (gcd a b) a b H5)) n));try (rewrite <- H7;trivial).
rewrite power_mult_lemma1 in H8;rewrite H8 in H4.
rewrite power_mult_lemma1 in H9;rewrite H9 in H4.
rewrite mult_lemma7 in H4.
assert ((power (gcd a b) n)<>O).
intro.
generalize (power_zero n (gcd a b) H10);intro.
apply (gcd_non_zero (gcd a b) a b);trivial.
generalize (mult_lemma6 (power p r * k * power (quo b (gcd a b) (gcd_div_r (gcd a b) a b H5)) n) (power (quo a (gcd a b) (gcd_div_l (gcd a b) a b H5)) n) (power (gcd a b) n) H10 H4).
fold ((power p r * k * power (quo b (gcd a b) (gcd_div_r (gcd a b) a b H5)) n)<>(power (quo a (gcd a b) (gcd_div_l (gcd a b) a b H5)) n)).
apply nth_root_irrat_aux;trivial.
apply gcd_rel_prime;apply (gcd_non_zero (gcd a b) a b);trivial.
Qed.

(** Generalization of the previous theorem *)
Theorem nth_root_irrational : forall (p k a b n q r:nat),(is_prime p)->(rel_prime p k)->(0<r)->(r<n)->(b<>0)->((power p (q*n+r))*k*(power b n) <> (power a n)).
intros.
intro.
rewrite power_plus_lemma1 in H4.
assert (divides a (power p q)).
apply prime_power_qn with n;try (auto with arith);try omega.
exists ((power p r)*k*(power b n)).
rewrite <- H4;ring.
assert (0<n);try omega.
elim H5;intro a';intro.
rewrite H7 in H4.
rewrite power_mult_lemma1 in H4;rewrite power_power_lemma1 in H4.
assert ((power p (q*n))<>0).
intro;apply not_prime_zero;generalize (power_zero (q*n) p H8);intro;rewrite H9 in H;trivial.
rewrite <- (mult_assoc (power p (q*n))) in H4;rewrite <- (mult_assoc (power p (q*n))) in H4.
generalize (mult_lemma6 (power p r*k*power b n) (power a' n) (power p (q*n)) H8 H4).
fold (power p r * k * power b n <> power a' n).
apply nth_root_irrat;trivial.
Qed.

(******************************************************************************)

(** lemmae about divisibility *)
Lemma divides_le : forall (a b:nat),(a<>O)->(divides a b)->(b<=a).
intros.
elim H0;intro q;intro.
replace b with (b*1);try ring.
rewrite H1.
apply mult_le_compat;try omega.
destruct q;omega.
Qed.

(** Euclide theorem (existence) *)
Theorem euclide : forall (a b:nat),(b<>O)->{q:nat & { r:nat | (a=b*q+r) /\ (r < b)}}.
intros.
apply (lt_wf_rec a (fun a:nat =>{q : nat &  {r : nat | a = b * q + r /\ r < b}})).
intros.
case (le_lt_dec b n);intro.
elim (H0 (n-b)).
intro q;intro.
elim p;intro r;intro.
exists (q+1);exists r.
split;try tauto.
rewrite (le_plus_minus b n);trivial.
elim p0;intros.
rewrite H1;ring.
omega.
exists 0;exists n.
split;try tauto.
ring.
Qed.

Definition quotient_euclide (a b:nat)(H:(b<>O)) := let (q,_) := (euclide a b H) in q.

Definition remainder_euclide (a b:nat)(H:(b<>O)) := let (_,e0) := (euclide a b H) in let (r,_) := e0 in r.

(** a div b where b<>0 *)
Lemma quo_rem_euclide : forall (a b:nat)(H:(b<>O)),a=b*(quotient_euclide a b H)+(remainder_euclide a b H).
unfold quotient_euclide;unfold remainder_euclide;intros.
generalize (euclide a b H);intros.
elim s;intro q;intro.
elim p;intro r;intro.
tauto.
Qed.

(** a mod b where b<>0 *)
Lemma rem_euclide : forall (a b:nat)(H:(b<>O)),(remainder_euclide a b H)<b.
unfold remainder_euclide;intros.
generalize (euclide a b H);intros.
elim s;intro q;intro.
elim p;intro r;intro.
tauto.
Qed.

(** Euclide division is unique *)
Lemma euclide_unique : forall (a b q r q' r':nat),(b<>O)->a=b*q+r->a=b*q'+r'->r<b->r'<b->(q=q')/\(r=r').
intros.
rewrite H1 in H0.
case (lt_eq_lt_dec q q');intro.
case s;intro.
rewrite (le_plus_minus q q') in H0;try (auto with arith).
rewrite mult_plus_distr_l in H0.
assert (b*(q'-q)+r' = r).
apply plus_reg_l with (b*q).
rewrite plus_assoc;trivial.
assert (0<(q'-q));try omega.
assert (b<=b*(q'-q));try omega.
case (mult_O_le b (q'-q));intro;try omega.
rewrite mult_comm;trivial.
split;try tauto.
rewrite <- e in H0.
symmetry;apply plus_reg_l with (b*q);trivial.
rewrite (le_plus_minus q' q) in H0;try (auto with arith).
rewrite mult_plus_distr_l in H0.
assert (r'=(b*(q-q')+r)).
apply plus_reg_l with (b*q').
rewrite plus_assoc;trivial.
assert (0<(q-q'));try omega.
assert (b<=b*(q-q'));try omega.
case (mult_O_le b (q-q'));intro;try omega.
rewrite mult_comm;trivial.
Qed.

(** if b<>0, then b | a iff a mod b = 0 *) 
Lemma divides_euclide : forall (a b:nat)(H:(b<>O)),((divides a b)<->((remainder_euclide a b H)=O)).
intros.
red.
split;intro.
generalize (quo_rem_euclide a b H);intro.
generalize (rem_euclide a b H);intro.
elim H0;intro q;intro.
assert (a=b*q+0).
rewrite plus_comm;simpl;trivial.
assert (0<b);try omega.
generalize (euclide_unique a b (quotient_euclide a b H) (remainder_euclide a b H) q 0 H H1 H4 H2 H5).
intros;tauto.
generalize (quo_rem_euclide a b H).
rewrite H0;rewrite plus_comm;simpl.
intro;exists (quotient_euclide a b H);trivial.
Qed.

(** divisibility is decidable *)
Lemma divides_dec : forall (a b:nat),{divides a b}+{~(divides a b)}.
intros.
case (eq_nat_dec b 0).
case (eq_nat_dec a 0);intros.
rewrite e;left;apply zero_max_div.
right;rewrite e;intro.
elim H;intro q;intro.
simpl in H0;apply n;trivial.
intro.
case (eq_nat_dec (remainder_euclide a b n) 0);[left | right];intros;elim (divides_euclide a b n);auto.
Qed.

(** we show that if b<>0, then gcd(a,b)=gcd(b,a mod b) *)
Lemma gcd_euclide : forall (d a b:nat)(H:(b<>0)),(is_gcd d a b)<->(is_gcd d b (remainder_euclide a b H)).
intros.
generalize (quo_rem_euclide a b H);intro.
red;split;intro.
rewrite H0 in H1.
elim H1;intros.
unfold is_gcd;unfold is_cd.
elim H2;intros.
split.
split;try tauto.
elim H4;intro q;intro.
elim H5;intro q';intro.
replace (b*(quotient_euclide a b H)) with (d*q'*(quotient_euclide a b H)) in H6.
assert ((remainder_euclide a b H)=(d*q-d*q'*(quotient_euclide a b H))).
rewrite <- H6;rewrite minus_plus;trivial.
rewrite <- mult_assoc in H8;rewrite <- mult_minus_lemma2 in H8.
exists (q-q'*(quotient_euclide a b H));trivial.
rewrite <- H7;trivial.
intros.
elim H6;intros.
apply H3.
unfold is_cd;split;try tauto.
elim H7;intro q;intro.
elim H8;intro q';intro.
rewrite H10.
replace (b*(quotient_euclide a b H)) with (d'*q*(quotient_euclide a b H)).
rewrite <- mult_assoc;rewrite <- mult_plus_distr_l.
exists (q*(quotient_euclide a b H)+q');trivial.
rewrite <- H9;trivial.
unfold is_gcd;unfold is_cd.
unfold is_gcd in H1;unfold is_cd in H1.
elim H1;intros.
elim H2;intros.
rewrite H0.
split.
split;try tauto.
elim H4;intro q;intro.
elim H5;intro q';intro.
rewrite H7.
replace (b*(quotient_euclide a b H)) with (d*q*(quotient_euclide a b H)).
rewrite <- mult_assoc;rewrite <- mult_plus_distr_l.
exists (q*(quotient_euclide a b H)+q');trivial.
rewrite <- H6;trivial.
intros.
apply H3.
split;try tauto.
elim H6;intros.
elim H7;intro q;intro.
elim H8;intro q';intro.
assert ((remainder_euclide a b H)=b*(quotient_euclide a b H)+(remainder_euclide a b H)-b*(quotient_euclide a b H)).
rewrite minus_plus;trivial.
rewrite H9 in H11.
exists (q-q'*(quotient_euclide a b H)).
rewrite mult_minus_lemma2;rewrite mult_assoc.
rewrite <- H10;trivial.
Qed.

(** we give a "more efficient" algorithm to compute gcd(a,b) *)
Lemma gcd_exists_prod_bis : forall (x:nat*nat),{d:nat | (is_gcd d (fst x) (snd x))}.
apply (induction_ltof2 (nat*nat) f (fun x:nat*nat => {d:nat | (is_gcd d (fst x) (snd x))})).
unfold ltof;unfold f;intros.
case (lt_eq_lt_dec (fst x) (snd x));intro.
case s;intro.
case (eq_nat_dec (fst x) 0);intro.
rewrite e;exists (snd x);apply gcd_zero.
elim (H ((fst x),(remainder_euclide (snd x) (fst x) n)));simpl.
intro d;intro.
exists d.
apply gcd_sym.
elim (gcd_euclide d (snd x) (fst x) n);auto.
generalize (rem_euclide (snd x) (fst x) n);try omega.
rewrite e;exists (snd x);apply gcd_refl.
case (eq_nat_dec (snd x) 0);intro.
rewrite e;exists (fst x);apply gcd_sym;apply gcd_zero.
elim (H ((snd x),(remainder_euclide (fst x) (snd x) n)));simpl.
intro d;intro.
exists d.
elim (gcd_euclide d (fst x) (snd x) n);auto.
generalize (rem_euclide (fst x) (snd x) n);try omega.
Qed.

(** efficient algorithm to compute gcd(a,b) *)
Theorem gcd_exists_bis : forall (a b:nat),{d:nat | (is_gcd d a b)}.
intros.
elim (gcd_exists_prod_bis (a,b));intro d;simpl;intros.
exists d;trivial.
Qed.

(** study of prime numbers *)
Lemma divides_prime : forall (a:nat),(exists p:nat,(p<>a)/\(is_prime p)/\(divides a p)) -> ~(is_prime a).
intros;intro.
elim H;intro p;intro.
elim H1;intros.
elim H3;intros.
unfold is_prime in H0.
elim H0;intros.
unfold is_prime in H4.
elim H4;intros.
case (H7 p H5);auto.
Qed.

(** p<>1 is prime if forall n, gcd(n,p)=1 or gcd(n,p)=p *)
Lemma gcd_prime : forall (p:nat),(p<>1)->(forall (d a:nat),(is_gcd d a p)->(d=1)\/(d=p))->(is_prime p).
intros.
split;try tauto.
intro d';intro.
assert (is_gcd d' d' p).
unfold is_gcd;unfold is_cd.
split;[split;[apply divides_refl | tauto] | tauto].
eapply H0;apply H2.
Qed.

(** p is prime iff forall n, n<>1 /\ n<>p -> not(a | p) *)
Lemma prime_cond : forall (p:nat),((p<>1)/\(forall (a:nat),(a<>1)->(a<>p)->~(divides p a))<->(is_prime p)).
split;intros.
elim H;intros.
split;try tauto.
intros.
case (eq_nat_dec d 1);intro;try tauto.
case (eq_nat_dec d p);intro;try tauto.
elim (H1 d n n0 H2).
elim H;intros.
split;try tauto.
intros;intro.
elim (H1 a H4);auto.
Qed.

(** if a property about integer is decidable then it is decidable if there is an integer less than n that satisfies this property *)
Lemma dec_impl_lt_dec : forall (P:nat->Prop),(forall (n:nat),{(P n)}+{~(P n)})->(forall (m:nat),{n:nat | (n<m)/\(P(n))}+{(forall (n:nat),(n<m)->~(P n))}).
intros.
induction m.
right;intros;inversion H0.
case (H m);intro.
left;exists m;split;try (auto with arith).
case IHm;intro.
elim s;intro n0;intro.
left;exists n0;split;[omega | tauto].
right;intros.
inversion H0;trivial.
apply n0;omega.
Qed.

(** forall n, either forall p, p<>1 /\ p<>n -> not(p | n) or there is p such that p<>1 and p<>n and p | n *) 
Lemma divides_nat : forall (n:nat),{p:nat | (p<>1)/\(p<>n)/\(divides n p)}+{forall (p:nat),(p<>1)->(p<>n)->~(divides n p)}.
intros.
case (dec_impl_lt_dec (fun p => (p<>1)/\(divides n p))) with n;intros.
case (divides_dec n n0);intro.
case (eq_nat_dec n0 1);intros.
right;intro;tauto.
left;tauto.
right;tauto.
elim s;intros.
left;exists x.
split;try tauto.
split;try tauto.
omega.
case (eq_nat_dec n 0);intro.
rewrite e;left;exists 2.
split;try (intro;discriminate).
split;try (intro;discriminate).
apply zero_max_div.
right;intros.
case (lt_eq_lt_dec p n);intro.
case s;intro;[red in n0;intro;apply n0 with p;tauto | auto].
intro;generalize (divides_le n p n1 H1);omega.
Qed.

(** a number n is either prime or not (it is decidable) *)
Lemma prime_dec : forall (n:nat),{is_prime n}+{~(is_prime n)}.
intro.
case (divides_nat n);intro.
elim s;intros.
right;intro.
unfold is_prime in H.
elim H;intros.
elim (H1 x);try tauto.
case (eq_nat_dec n 1);intro.
right;unfold is_prime;tauto.
left;unfold is_prime.
split;trivial.
intros.
case (eq_nat_dec d 1);try tauto.
case (eq_nat_dec d n);try tauto;intros.
elim (n0 d n3 n2 H).
Qed.

(** if n is not prime then either n = 1 or there is a prime number p such that p | n *)
Lemma not_prime_impl_prime_divides : forall (n:nat),(~(is_prime n)->({p:nat | (is_prime p)/\(divides n p)}+{n=1})).
intro.
apply (lt_wf_rec n (fun n:nat => ~(is_prime n)->({p:nat | (is_prime p)/\(divides n p)}+{n=1})));intros.
case (eq_nat_dec n0 1);try tauto;intro.
case (eq_nat_dec n0 0);intro.
left;exists 2.
split;[apply is_prime_2 | rewrite e;apply zero_max_div].
case (divides_nat n0);intro.
elim s;intro d;intro.
elim p;intros.
elim H2;intros.
assert (d<n0).
generalize (divides_le n0 d n2 H4);omega.
case (prime_dec d);intro.
left;exists d;tauto.
elim (H d H5 n3);try tauto.
intro.
elim a;intro q;intro.
left;exists q.
split;try tauto.
apply divides_trans with d;try tauto.
elim (prime_cond n0);intros;elim H0;apply H1;auto.
Qed.

(** if n<>0 and p is a prime number then there is m such that p^m | n and not(p^(m+1) | n) *)
Lemma nat_factor : forall (n p:nat),(is_prime p)->(n<>0)->{m:nat | (divides n (power p m))/\~(divides n (power p (m+1)))}.
intros n p H.
apply (lt_wf_rec n (fun n:nat => n <> 0 -> {m : nat | divides n (power p m) /\ ~ divides n (power p (m + 1))}));intros.
case (divides_dec n0 p);intro.
generalize (quo_is_quo n0 p d);intro.
elim (H0 (quo n0 p d)).
intro m;intros.
exists (m+1).
elim p0;intros.
elim H3;intros.
rewrite H5 in H2;rewrite mult_assoc in H2.
rewrite plus_comm.
split;simpl.
exists x;trivial.
rewrite plus_comm;simpl.
rewrite (mult_comm p (power p m));rewrite mult_assoc;intro.
elim H6;intros.
rewrite H2 in H7.
assert (p<>0).
intro.
rewrite H8 in H.
apply not_prime_zero;trivial.
assert ((power p m)*x=(power p m)*p*x0).
apply mult_lemma6 with p;trivial.
rewrite mult_assoc;rewrite H7;ring.
rewrite <- H5 in H9;rewrite (mult_comm (power p m) p) in H9.
apply H4.
rewrite plus_comm;simpl.
exists x0;trivial.
rewrite mult_comm in H2;rewrite H2;apply mult_lemma3.
intro.
apply H1;rewrite H2;rewrite H3;trivial.
elim H.
intros.
destruct p;omega.
intro;apply H1.
rewrite H2;rewrite H3;ring.
exists 0;simpl.
split.
apply one_min_div.
rewrite mult_comm;simpl;rewrite plus_comm;simpl;trivial.
Qed.

(** if 1<p and 0<m then p^m>1 *)
Lemma power_lt : forall (p m:nat),(1<p)->(0<m)->1<(power p m).
induction m;simpl;try omega;intros.
destruct m;simpl;try omega.
simpl in IHm.
assert (1 < p*(power p m)).
apply IHm;auto with arith.
rewrite mult_comm.
apply lt_trans with (1*p);try omega.
apply mult_lt_compat_r;try omega.
Qed.

(** if n>1 then there is a prime number p and two integers m and q with m>0 and q<n such that n=p^m*q and gcd(p,q)=1 *)
Lemma nat_factor_prime : forall (n:nat),(n<>0)->(n<>1)->{p:nat & {m:nat & {q:nat | (is_prime p)/\(m>0)/\(n=(power p m)*q)/\(is_gcd 1 p q)/\(q<n)}}}.
intros.
case (prime_dec n);intro.
exists n;exists 1;exists 1;simpl.
split;try tauto.
split;try omega.
split;try ring.
split;try omega.
apply gcd_sym;apply gcd_one.
elim (not_prime_impl_prime_divides n n0);intro;try tauto.
elim a;intro p;intro.
elim p0;intros.
elim (nat_factor n p H1).
intro m;intro.
elim p1;intros.
exists p;exists m;exists (quo n (power p m) H3).
split;trivial.
split.
destruct m;try omega.
elim H4;simpl;rewrite mult_comm;simpl;rewrite plus_comm;simpl;trivial.
split.
apply (quo_is_quo n (power p m) H3).
generalize (gcd_is_gcd p (quo n (power p m) H3));intro.
assert ((gcd p (quo n (power p m) H3))=1).
case (prime_gcd (gcd p (quo n (power p m) H3)) p (quo n (power p m) H3));trivial.
apply gcd_sym;trivial.
intro.
rewrite H6 in H5;elim H5;intros.
elim H7;intros.
elim H4;rewrite plus_comm;simpl.
generalize (quo_is_quo n (power p m) H3);intro.
elim H10;intros.
rewrite H12 in H11.
exists x;rewrite H11;ring.
rewrite H6 in H5;trivial.
generalize (quo_is_quo n (power p m) H3);intro.
split;trivial.
rewrite H7;rewrite mult_comm;apply mult_lemma3.
intro.
apply H;rewrite H7.
rewrite H8;ring.
red;apply power_lt.
destruct p.
elim (not_prime_zero H1).
destruct p;try omega.
elim H1;tauto.
destruct m;try omega.
elim H4;simpl;rewrite mult_comm;simpl;rewrite plus_comm;simpl;trivial.
trivial.
Qed.

(******************************************************************************)

(** we now concentrate our efforts to show that each natural number n >= 1 has a unique decomposition in prime factors upto permutations *)
Require Import List.

(** this compute p1^n1*p2^n2*...*pm^nm if l = (p1,n1)::(p2,n2)::...::(pm,nm) *)
Fixpoint refactor (l:(list (nat*nat))) {struct l} : nat :=
match l with
nil => 1
| (cons (p,n) tail) => (power p n)*(refactor tail)
end.

(** a list l is well-formed (is a factorisation) if whenever (p,n) is in l then p is prime and n>0 and p is relatively prime to (refactor (l\(p,n))) and l\(p,n) is well-formed *)
Inductive is_wf : (list (nat*nat))->Prop :=
nil_is_wf : (is_wf nil)
|cons_is_wf : forall (p n:nat)(tail:(list (nat*nat))),(is_prime p)->(n>0)->(is_wf tail)->(rel_prime p (refactor tail))->(is_wf (cons (p,n) tail)).

(** this gives an inefficient algorithm to compute a factorisation of n *)
Lemma factorisation : forall (n:nat),{l:(list (nat*nat)) | (is_wf l)/\n=(refactor l)}+{n=0}.
intro.
case (eq_nat_dec n 0);intro.
right;trivial.
case (eq_nat_dec n 1).
intro;left;exists (nil (A:=nat*nat)).
split;[apply nil_is_wf | simpl;trivial].
generalize n0.
apply (lt_wf_rec n (fun n:nat => n<>0 -> n <> 1 -> {l : list (nat * nat) | is_wf l /\ n = refactor l}+{n=0}));intros.
elim (nat_factor_prime n1 H0 H1).
intro p;intro.
elim p0;intro m;intro.
elim p1;intro q;intro.
elim p2;intros.
elim H3;intros.
elim H5;intros.
elim H7;intros.
case (eq_nat_dec q 1);intro.
left;exists (cons (p,m) nil);simpl;rewrite e in H6.
split;trivial.
apply cons_is_wf;auto;try (apply nil_is_wf).
unfold rel_prime;simpl;rewrite e in H8;trivial.
assert (q<>0).
intro;rewrite H10 in H6;rewrite mult_comm in H6;simpl in H6;auto.
elim (H q H9 H10 n2).
intro.
elim a;intro l;intro.
elim p3;intros.
left;exists (cons (p,m) l);simpl;rewrite H12 in H6;split;trivial.
apply cons_is_wf;auto.
rewrite <- H12;unfold rel_prime;trivial.
intro;tauto.
Qed.

(** we now show that a factorisation is unique upto permutations *)
(** if (p,n) is in l then p^n | (refactor l) *)
Lemma factor_divides_refactor : forall (x:nat*nat)(l:list (nat*nat)),(In x l)->(divides (refactor l) (power (fst x) (snd x))).
induction l;simpl;try tauto.
intro.
case H;intro.
destruct a.
rewrite <- H0;simpl.
exists (refactor l);trivial.
destruct a.
elim (IHl H0);intros.
rewrite H1.
exists ((power n n0)*x0);ring.
Qed.

(** if p is prime and p | (refactor l) then there is n such that (p,n) is in l *)
Lemma prime_divides_refactor : forall (p:nat)(l:list (nat*nat)),(is_prime p)->(is_wf l)->(divides (refactor l) p)->(exists m:nat,(In (p,m) l)).
induction l;simpl;intros.
assert (p=1).
apply divides_antisym;trivial.
apply one_min_div.
rewrite H2 in H;elim H;tauto.
destruct a.
case (divides_dec (power n n0) p);intro.
generalize (prime_power p n0 n H d);intro.
assert (n=p).
inversion H0.
elim H6.
intros.
elim H;intros.
case (H11 p H2);try tauto;try omega.
exists n0;rewrite <- H3;left;trivial.
inversion H0.
elim (IHl H H7).
intros;exists x;tauto.
apply gauss with (power n n0);trivial.
apply rel_prime_sym;apply prime_div_gcd;trivial.
Qed.

(** 1^n = 1 *)
Lemma power_one : forall (n:nat),(power 1 n)=1.
induction n;simpl;trivial.
rewrite IHn;ring.
Qed.

(** it is decidable to say if a and b are relatively prime *)
Lemma rel_prime_dec : forall (a b:nat),{rel_prime a b}+{~(rel_prime a b)}.
intros.
unfold rel_prime.
generalize (gcd_is_gcd a b);intro.
case (eq_nat_dec (gcd a b) 1);intro.
left;rewrite e in H;trivial.
right;intro;apply n;apply (gcd_unique (gcd a b) 1 a b);trivial.
Qed.

(** if gcd(a,b)=1 and gcd(a,c)=1 then gcd(a,b*c)=1 *)
Lemma rel_prime_mult : forall (a b c:nat),(rel_prime a b)->(rel_prime a c)->(rel_prime a (b*c)).
intros.
split.
split;try (apply one_min_div).
intros.
elim H1;intros.
case (rel_prime_dec b d');intro.
assert (divides c d').
apply gauss with b;trivial.
elim H0;intros.
apply H6;unfold is_cd;tauto.
generalize (gcd_is_gcd b d');intro.
assert ((gcd b d')<>1).
intro;apply n.
unfold rel_prime;rewrite <- H5;trivial.
generalize (gcd_div_l (gcd b d') b d' H4);intro.
generalize (gcd_div_r (gcd b d') b d' H4);intro.
assert (divides a (gcd b d')).
apply divides_trans with d';[apply H2 | apply H7].
elim H5.
apply divides_antisym.
apply one_min_div.
elim H;intros;apply H10;unfold is_cd;tauto.
Qed.

(** if gcd(a,b*c)=1 then gcd(a,b)=1 and gcd(a,c)=1 *)
Lemma mult_rel_prime : forall (a b c:nat),(rel_prime a (b*c))->((rel_prime a b)/\(rel_prime a c)).
intros.
split;split;[split | intros | split | intros];try (apply one_min_div);elim H0;intros;elim H;intros;apply H4;split;trivial;elim H2;intro q;intro;rewrite H5;[exists (q*c) | exists (q*b)];ring.
Qed.

(** if gcd(a,d)=1 then gcd(a,d^n)=1 *)
Lemma rel_prime_power : forall (d a n:nat),(rel_prime a d)->(rel_prime a (power d n)).
induction n;simpl;intros.
unfold rel_prime;apply gcd_sym;apply gcd_one.
generalize (IHn H);intro.
apply rel_prime_mult;trivial.
Qed.

(** if n>0 and gcd(a,d^n)=1 then gcd(a,d)=1 *)
Lemma power_rel_prime : forall (d a n:nat),(n>0)->(rel_prime a (power d n))->(rel_prime a d).
destruct n;simpl;intros.
inversion H.
elim (mult_rel_prime a d (power d n));auto.
Qed.

(** if n>0 and m>0 then gcd(a^n,b^m)=1 iff gcd(a,b)=1 *)
Lemma power_power_rel_prime : forall (a n b m:nat),(n>0)->(m>0)->((rel_prime (power a n) (power b m))<->(rel_prime a b)).
split;intro.
apply power_rel_prime with m;trivial;apply rel_prime_sym;apply power_rel_prime with n;trivial;apply rel_prime_sym;trivial.
apply rel_prime_power;apply rel_prime_sym;apply rel_prime_power;apply rel_prime_sym;trivial.
Qed.

(** if x>1 and x^m | x^n then m<=n *)
Lemma power_divides_power : forall (x n m:nat),(x>1)->(divides (power x n) (power x m))->(m<=n).
intros.
case (le_lt_dec m n);trivial.
intro.
generalize (le_plus_minus n m);intro.
rewrite H1 in H0;try omega.
elim H0;intro q;rewrite power_plus_lemma1;intro.
assert (1=(power x (m-n))*q).
apply mult_lemma6 with (power x n).
intro;generalize (power_zero n x H3);omega.
rewrite mult_assoc;rewrite <- H2;ring.
symmetry in H3;elim (mult_lemma5 (power x (m-n)) q H3);intros.
case (eq_nat_dec (m-n) 0);intro;try omega.
assert (x=1);try omega.
apply divides_antisym;[apply one_min_div | rewrite <- H4;apply power_divides_lemma1;omega].
Qed.

(** if (p,n) is in l and l is well-formed then p is prime and n>0 *)
Lemma in_wf : forall (l:list (nat*nat))(p n:nat),(In (p,n) l)->(is_wf l)->(is_prime p)/\(n>0).
induction l;simpl;try tauto.
intros;destruct a.
inversion H0.
case (in_inv H);intros.
inversion H8;rewrite <- H10;rewrite <- H11;try tauto.
apply IHl;trivial.
Qed.

(** if (p,n) is in l and (q,m)::l is well-formed then gcd(p,q)=1 *)
Lemma rel_prime_wf : forall (l:list (nat*nat))(p n q m:nat),(In (p,n) l)->(is_wf ((q,m)::l))->(rel_prime p q).
induction l;simpl;intros;try tauto.
destruct a.
inversion H0;case (in_inv H);intro.
rewrite H8 in H7;simpl in H7.
elim (in_wf ((n0,n1)::l) p n H);trivial;intros.
elim (mult_rel_prime q (power p n) (refactor l) H7);intros.
apply rel_prime_sym;apply power_rel_prime with n;trivial.
apply (IHl p n q m);trivial.
apply cons_is_wf;trivial.
inversion H6;trivial.
simpl in H7;elim (mult_rel_prime q (power n0 n1) (refactor l) H7);auto.
Qed.

Inductive is_pwd : list (nat*nat) -> Prop :=
nil_is_pwd : (is_pwd nil)
|cons_is_pwd : forall (p n:nat)(tail:list (nat*nat)),(is_pwd tail)->(forall (n:nat),~(In (p,n) tail))->(is_pwd ((p,n)::tail)).

(** if l is well-formed then the first projection of l is pairwise distinct (pwd) *)
Lemma wf_impl_pwd : forall (l:list (nat*nat)),(is_wf l)->(is_pwd l).
induction l;intro.
apply nil_is_pwd.
destruct a.
inversion H.
apply cons_is_pwd;auto.
intros;intro.
assert (rel_prime n n).
eapply rel_prime_wf;[apply H7 | apply H].
generalize (gcd_refl n);intro.
unfold rel_prime in H8.
assert (1=n).
eapply gcd_unique;eauto.
rewrite <- H10 in H3;elim H3;tauto.
Qed.

(** if p is prime, n>0, l is well-formed and p^n | (refactor l) and (p,m) is in l then n<=m *)
Lemma prime_power_divides_refactor : forall (p n m:nat)(l:list (nat*nat)),(is_prime p)->(n>0)->(is_wf l)->(divides (refactor l) (power p n))->(In (p,m) l)->(n<=m).
induction l;simpl;intros;try tauto.
case H3;intro.
rewrite H4 in H2.
rewrite H4 in H1;inversion H1.
generalize (rel_prime_power p (refactor l) n (rel_prime_sym p (refactor l) H11));intro.
rewrite mult_comm in H2.
generalize (gauss (power p n) (refactor l) (power p m) H12 H2);intro.
apply power_divides_power with p;trivial.
destruct p.
elim (not_prime_zero H8).
elim H8;omega.
destruct a.
inversion H1.
apply IHl;trivial.
apply gauss with (power n0 n1);trivial.
apply rel_prime_power;apply rel_prime_sym;apply rel_prime_power.
apply (rel_prime_wf l p m n0 n1);trivial.
Qed.

(** we define some notions on lists such as being a permutation of *)
Section Permutation.

(** insertion x l l' iff l' is l where x has been inserted somewhere *)
Inductive insertion (A:Set) : A -> list A -> list A -> Prop :=
head_insertion : forall (x:A)(l:list A),(insertion A x l (x::l))
|tail_insertion : forall (x y:A)(l l':list A),(insertion A x l l')->(insertion A x (y::l) (y::l')).

(** if (insertion x l l') then x is in l' *)
Lemma insertion_in : forall (A:Set)(x:A)(l l':list A),(insertion A x l l')->(In x l').
intros.
induction H;simpl;tauto.
Qed.

(** if (insertion x l l') then l is included in l' *)
Lemma insertion_inclusion : forall (A:Set)(x:A)(l l':list A),(insertion A x l l')->forall (y:A),(In y l)->(In y l').
induction l;simpl;try tauto;intros.
inversion H;simpl;try tauto.
case H0;try tauto.
right;auto.
Qed.

(** if x is in l, then there is l' such that (insertion x l' l) *)
Lemma in_insertion : forall (A:Set)(x:A)(l:list A),(In x l)->exists l':list A,(insertion A x l' l).
induction l;simpl;try tauto;intros.
case H;intro.
rewrite H0;exists l;apply head_insertion.
elim (IHl H0);intro l';intro.
exists (a::l');apply tail_insertion;trivial.
Qed.

(** if (insertion x l l') and y is in l' then y=x or y is in l *)
Lemma in_insertion_inv : forall (A:Set)(x y:A)(l l':list A),(insertion A y l l')->(In x l')->(x=y)\/(In x l).
intros.
induction H;simpl in H0.
case H0;intro H1;try (symmetry in H1);tauto.
case H0;simpl;intro;tauto.
Qed.

(** a list is a set iff all the elements are pairwise distinct *)
Inductive is_set (A:Set) : list A->Prop :=
nil_is_set : (is_set A nil)
|cons_is_set : forall (x:A)(l:list A),(is_set A l)->~(In x l)->(is_set A (x::l)).

(** if (insertion x l l') and l' is a set then l is a set *)
Lemma is_set_insertion : forall (A:Set)(l l':list A)(x:A),(insertion A x l l')->(is_set A l')->(is_set A l).
induction 1;intros.
inversion H;trivial.
inversion H0.
apply cons_is_set.
apply IHinsertion;trivial.
intro;apply H4;apply (insertion_inclusion A x l l');trivial.
Qed.

(** if (insertion x l l') and l' is a set then x is not in l *)
Lemma is_set_insertion_in : forall (A:Set)(l l':list A)(x:A),(insertion A x l l')->(is_set A l')->~(In x l).
induction l;simpl;try tauto;intros.
inversion H;rewrite <- H3 in H0;inversion H0.
simpl in H7;trivial.
intro.
case H10;intro.
apply H9;rewrite H11;eapply insertion_in;apply H5.
elim (IHl l'0 x H5 H8 H11).
Qed.

(** l' is a permutation of l *)
Inductive is_permutation (A:Set) : list A->list A->Prop :=
nil_is_permutation : (is_permutation A nil nil)
|cons_is_permutation : forall (l l':list A),(is_permutation A l l')->forall (x:A)(l'':(list A)),(insertion A x l' l'')->(is_permutation A (x::l) l'').

(** if l and l' have the same content and are pairwise distinct then l' is a permutation of l *)
Lemma is_set_eq_impl_permutation : forall (A:Set)(l l':list A),(forall (x:A),(In x l)<->(In x l'))->(is_set A l)->(is_set A l')->(is_permutation A l l').
induction l;intros;simpl in H.
destruct l'.
apply nil_is_permutation.
elim (H a);intros.
elim H3;simpl;tauto.
inversion H0.
symmetry in H2;rewrite H2 in H;elim (H a);intros.
rewrite H2 in H6;rewrite H2.
assert (In x l');auto.
elim (in_insertion A x l' H8).
intro l'';intro.
apply cons_is_permutation with l'';trivial.
apply IHl;trivial.
split;intro.
elim (H x0);intros.
elim (in_insertion_inv A x0 x l'' l');auto.
intro;rewrite H13 in H10;rewrite H2 in H5;tauto.
elim (H x0);intros.
case H12;try tauto.
apply (insertion_inclusion A x l'' l');trivial.
intro;rewrite <- H13 in H10.
elim (is_set_insertion_in A l'' l' x);trivial.
eapply is_set_insertion;eauto.
Qed.

(** is_permutation is reflexive *)
Lemma is_permutation_refl : forall (A:Set)(l:list A),(is_permutation A l l).
induction l.
apply nil_is_permutation.
eapply cons_is_permutation;[apply IHl | apply head_insertion].
Qed.

(** if l' is l where x has been inserted then l' is a permutation of x::l *)
Lemma insertion_is_permutation : forall (A:Set)(l l':list A)(x:A),(insertion A x l l')->(is_permutation A (x::l) l').
induction 1.
apply cons_is_permutation with l;[apply is_permutation_refl | apply head_insertion].
apply cons_is_permutation with (y::l);[apply is_permutation_refl | apply tail_insertion;trivial].
Qed.

(** if l1 is l0 where x has been inserted and l2 is l1 where y has been inserted then there is l3 such that l3 is l0 where y has been inserted and l2 is l3 where x has been inserted *)
Lemma insertion_trans : forall (A:Set)(l0 l1:list A)(x:A),(insertion A x l0 l1)->forall (l2:list A)(y:A),(insertion A y l1 l2)->exists l3:list A,(insertion A y l0 l3)/\(insertion A x l3 l2).
induction 1;intros.
inversion H.
exists (y::l);split;[apply head_insertion | apply tail_insertion;apply head_insertion].
exists l';split;[trivial | apply head_insertion].
inversion H0.
exists (y0::y::l);split;[apply head_insertion | apply tail_insertion;apply tail_insertion;trivial].
elim (IHinsertion l'0 y0 H5);intro l3;intro.
elim H6;intros.
exists (y::l3);split;[apply tail_insertion | apply tail_insertion];trivial.
Qed.

(** if l1 is a permutation of l0 and then l1 where x has been inserted is a permutation of l0 where x has been inserted *)
Lemma permutation_insertion : forall (A:Set)(l0 l1:list A),(is_permutation A l0 l1)->forall (x:A)(l2 l3:list A),(insertion A x l0 l2)->(insertion A x l1 l3)->(is_permutation A l2 l3).
induction 1;intros.
inversion H;inversion H0;apply is_permutation_refl.
inversion H1.
apply cons_is_permutation with l'';trivial.
apply cons_is_permutation with l';trivial.
elim (insertion_trans A l' l'' x H0 l3 x0 H2).
intro l4;intro.
elim H8;intros.
apply cons_is_permutation with l4;trivial.
eapply IHis_permutation;eauto.
Qed.

(** is_permutation is symmetric *)
Lemma is_permutation_sym : forall (A:Set)(l l':list A),(is_permutation A l l')->(is_permutation A l' l).
induction 1;[apply nil_is_permutation | eapply permutation_insertion;eauto;apply head_insertion].
Qed.

Lemma permutation_in : forall (A:Set)(l l':list A),(is_permutation A l l')->forall (x:A),(In x l)<->(In x l').
induction l;simpl;intros.
inversion H;simpl;tauto.
inversion H;simpl.
split;intro.
case H5;intro.
eapply insertion_in;rewrite H6 in H4;apply H4.
elim (IHl l'0 H2 x);intros.
eapply insertion_inclusion;eauto.
case (in_insertion_inv A x a l'0 l' H4 H5);intro.
rewrite H6;tauto.
elim (IHl l'0 H2 x);intros.
right;auto.
Qed.

Lemma permutation_insertion_aux : forall (A:Set)(l l' l'':list A)(x:A),(insertion A x l l')->(insertion A x l l'')->(is_permutation A l' l'').
intros.
eapply permutation_insertion;eauto.
apply is_permutation_refl.
Qed.

Lemma length_recursion : forall (A:Set),forall (P:list A->Prop),(forall (x:list A),(forall (y:list A),(length y)<(length x)->(P y))->(P x))->(forall (a:list A),(P a)).
intros.
apply (well_founded_ind (well_founded_ltof (list A) (fun l:list A => length l)));unfold ltof;auto.
Qed.

Lemma insertion_length : forall (A:Set)(l l':list A)(x:A),(insertion A x l l')->(length l')=(S (length l)).
induction 1;simpl;trivial.
rewrite IHinsertion;trivial.
Qed.

Lemma permutation_length : forall (A:Set)(l l':list A),(is_permutation A l l')->(length l)=(length l').
induction 1;simpl;trivial.
generalize (insertion_length A l' l'' x H0);intro;congruence.
Qed.

Lemma insertion_permutation_eq : forall (A:Set)(l l':list A)(x:A),(insertion A x l' l)->forall (l'':list A),(insertion A x l'' l)->(is_permutation A l' l'').
induction l;intros;inversion H.
inversion H0.
apply is_permutation_refl.
rewrite <- H4;destruct l.
inversion H8.
generalize (head_insertion A a0 l);intro.
assert (In x (a0::l)).
eapply insertion_in;apply H8.
case (in_insertion_inv A x a0 l (a0::l) H10 H11);intro.
rewrite H12;rewrite <- H12 in H10;rewrite <- H12 in H8;rewrite <- H12 in IHl.
assert (is_permutation A l l1).
eapply IHl;eauto.
eapply cons_is_permutation;eauto;apply head_insertion.
elim (in_insertion A x l H12);intro l2;intro.
generalize (tail_insertion A x a0 l2 l H13);intro.
assert (is_permutation A (a0::l2) l1).
eapply IHl;eauto.
apply is_permutation_sym;auto.
eapply cons_is_permutation;eauto.
rewrite H1 in H3.
inversion H0.
rewrite <- H9;apply insertion_is_permutation;trivial.
assert (is_permutation A l0 l1).
eapply IHl;eauto.
eapply cons_is_permutation;eauto;apply head_insertion.
Qed.

Lemma permutation_insertion_comm : forall (A:Set)(l1 l2:list A)(x:A),(insertion A x l1 l2)->forall (l4:list A),(is_permutation A l2 l4)->(exists l3:list A,(is_permutation A l1 l3) /\ (insertion A x l3 l4)).
intros A l1 l2.
generalize l1;clear l1.
induction l2;intros.
inversion H.
inversion H.
subst a.
subst x0.
subst l.
subst l2.
inversion H0.
subst x0.
subst l.
subst l''.
exists l'.
tauto.
subst x0.
subst l'.
subst a.
subst l1.
inversion H0.
subst x0.
subst l''.
subst l0.
elim (IHl2 l x H4 l' H3).
intro l3;intros.
elim H1;clear H1;intros.
elim (insertion_trans A l3 l' x H2 l4 y H6).
intro l5;intros.
elim H5;clear H5;intros.
exists l5.
split;trivial.
eapply cons_is_permutation.
apply H1.
trivial.
Qed.

Lemma permutation_insertion_permutation : forall (A:Set)(l l':list A),(is_permutation A l l')->forall (x:A)(l'':list A),(insertion A x l' l'')->forall (l''':list A),(is_permutation A l'' l''')->(is_permutation A (x::l) l''').
induction 1;intros.
inversion H.
rewrite <- H3 in H0.
trivial.
elim (permutation_insertion_comm A l'' l''0 x0 H1 l''' H2).
intro l1;intro.
elim H3;clear H3;intros.
eapply cons_is_permutation.
eapply IHis_permutation.
apply H0.
apply H3.
trivial.
Qed.

(** is_permutation is transitive *)
Lemma is_permutation_trans : forall (A:Set)(l l':list A),(is_permutation A l l')->forall (l'':list A),(is_permutation A l' l'')->(is_permutation A l l'').
induction l.
intros.
inversion H.
rewrite <- H2 in H0;trivial.
intros.
inversion H.
induction H5;inversion H0.
eapply cons_is_permutation;try (apply IHl with l1;eauto);trivial.
eapply permutation_insertion_permutation.
apply H3.
apply tail_insertion.
apply H5.
eapply cons_is_permutation.
apply H8.
apply H10.
Qed.
End Permutation.

(** if l and l' are well-formed and (refactor l)=(refactor l') then l is included in l' *)
Lemma factorisation_unique_upto_equiv_aux : forall (l l':list (nat*nat)),(is_wf l)->(is_wf l')->(refactor l)=(refactor l')->(forall (x n:nat),(In (x,n) l)->(In (x,n) l')).
intros l l';intro;intro;intro;intros p n;intro.
elim (in_wf l p n H2 H);intros.
generalize (factor_divides_refactor (p,n) l H2);simpl;intro.
rewrite H1 in H5.
assert (divides (refactor l') p).
apply divides_trans with (power p n);[trivial | apply power_divides_lemma1;auto with arith].
elim (prime_divides_refactor p l' H3 H0 H6);intro m;intro.
cut (n=m).
intro;rewrite H8;trivial.
apply le_antisym.
eapply prime_power_divides_refactor;eauto.
generalize (factor_divides_refactor (p,m) l' H7);simpl;intro.
rewrite <- H1 in H8.
apply prime_power_divides_refactor with p l;auto.
elim (in_wf l' p m);trivial.
Qed.

Lemma pwd_impl_set : forall (l:list (nat*nat)),(is_pwd l)->(is_set (nat*nat) l).
induction 1;[apply nil_is_set | apply cons_is_set;auto].
Qed.

(** if l and l' are well-formed and (refactor l)=(refactor l') then l' is a permutation of l *)
Lemma factorisation_unique_upto_perm : forall (l l':list (nat*nat)),(is_wf l)->(is_wf l')->(refactor l)=(refactor l')->(is_permutation (nat*nat) l l').
intros.
assert (forall (x n:nat),(In (x,n) l)->(In (x,n) l')).
apply factorisation_unique_upto_equiv_aux;trivial.
assert (forall (x n:nat),(In (x,n) l')->(In (x,n) l)).
apply factorisation_unique_upto_equiv_aux;auto.
apply is_set_eq_impl_permutation;intros;(try (apply pwd_impl_set;apply wf_impl_pwd;trivial));destruct x;split;auto.
Qed.

(** a list is a factorisation of n if l is well-formed and (refactor l)=n *)
Definition is_factorisation (n:nat)(l:list (nat*nat)) := (is_wf l)/\(n=(refactor l)).

(** The fondamental theorem of arithmetic *)
(** forall n<>0, there exists a factorisation l *)
Theorem factorisation_exists : forall (n:nat),n<>0->{l:list (nat*nat) | (is_factorisation n l)}.
intros.
case (factorisation n);intros;try tauto.
Qed.

(** if l and l' are two factorisation of n, then l' is a permutation of l *)
Theorem factorisation_unique_upto_permutation : forall (n:nat)(l l':list (nat*nat)),(is_factorisation n l)->(is_factorisation n l')->(is_permutation (nat*nat) l l').
unfold is_factorisation;intros.
elim H;intros.
elim H0;intros.
apply factorisation_unique_upto_perm;auto;congruence.
Qed.

Lemma wf_power_dec : forall (n:nat)(l:list (nat*nat)),(is_wf l)->(n>0)->{x:nat | (refactor l)=(power x n)}+{p:nat & {q:nat & {r:nat & {k:nat | (is_prime p)/\(0<r)/\(r<n)/\(refactor l)=(power p (q*n+r))*k/\(rel_prime p k)}}}}.
intro.
induction l;simpl;intros.
left;exists 1;rewrite power_one;trivial.
destruct a.
assert (n<>0);try omega.
generalize (quo_rem_euclide n1 n H1);intro.
case (eq_nat_dec (remainder_euclide n1 n H1) 0);intro.
rewrite e in H2;rewrite plus_comm in H2;simpl in H2.
case IHl;intros;trivial.
inversion H;trivial.
elim s;intro y;intro.
rewrite H2.
left;rewrite p;rewrite (mult_comm n);rewrite <- power_power_lemma1;rewrite <- power_mult_lemma1;exists (power n0 (quotient_euclide n1 n H1)*y);trivial.
elim s;intro p;intro.
elim p0;intro q;intro.
elim p1;intro r;intro.
elim p2;intro k;intro.
elim p3;intros.
elim H4;intros.
elim H6;intros.
elim H8;intros.
right.
exists p;exists q;exists r.
rewrite H9;rewrite mult_comm;rewrite <- mult_assoc.
exists (k*(power n0 n1)).
split;trivial.
split;trivial.
split;trivial.
split;trivial.
apply rel_prime_mult;trivial.
inversion H.
rewrite H9 in H17.
elim (mult_rel_prime n0 (power p (q*n+r)) k H17);intros.
apply rel_prime_power;apply rel_prime_sym;apply power_rel_prime with (q*n+r);trivial.
rewrite plus_comm;auto with arith.
right.
exists n0;exists (quotient_euclide n1 n H1);exists (remainder_euclide n1 n H1).
rewrite (mult_comm (quotient_euclide n1 n H1));rewrite <- H2.
exists (refactor l).
elim (in_wf ((n0,n1)::l) n0 n1);intros.
split;trivial.
split;trivial.
destruct (remainder_euclide n1 n H1);try tauto;auto with arith.
split;trivial.
apply rem_euclide.
split;trivial.
inversion H;trivial.
simpl;tauto.
trivial.
Qed.

(** let n be a natural number and m>1 then either n=x^m or there is a prime number p and three numbers q,r and k such that n = p^(q*m+r)*k with 0<r<m and gcd(p,k)=1 *)
Lemma is_power_m_dec : forall (n m:nat),(m>0)->{x:nat | n=(power x m)}+{p:nat & {q:nat & {r:nat & {k:nat | (is_prime p)/\(0<r)/\(r<m)/\n=(power p (q*m+r))*k/\(rel_prime p k)}}}}.
intros n m;intro H.
case (eq_nat_dec n 0);intro.
left;exists 0.
destruct m;simpl;try omega;trivial.
generalize (factorisation_exists n n0);intro.
elim H0;intro l;intro.
elim p;intros.
rewrite H2;apply wf_power_dec;trivial.
Qed.

(** let x and n be two numbers such that n > 0, then either the n-th root of x is a natural number of it is not rationnal *)
Theorem nth_root : forall (x n:nat),(n>0)->{y:nat | x=(power y n)}+{forall (a b:nat),(b<>0)->x*(power b n)<>(power a n)}.
intros.
case (is_power_m_dec x n H);intro;try tauto.
elim s;intro p;intro.
elim p0;intro q;intro.
elim p1;intro r;intro.
elim p2;intro k;intro.
right;intros.
assert (x=(power p (q*n+r))*k);try tauto.
rewrite H1;apply nth_root_irrational;tauto.
Qed.