1
(************************************************************************)
2
(* v * The Coq Proof Assistant / The Coq Development Team *)
3
(* <O___,, * CNRS-Ecole Polytechnique-INRIA Futurs-Universite Paris Sud *)
4
(* \VV/ **************************************************************)
5
(* // * This file is distributed under the terms of the *)
6
(* * GNU Lesser General Public License Version 2.1 *)
7
(************************************************************************)
9
(*i $Id: Zgcd_alt.v 10997 2008-05-27 15:16:40Z letouzey $ i*)
11
(** * Zgcd_alt : an alternate version of Zgcd, based on Euler's algorithm *)
14
Author: Pierre Letouzey
17
(** The alternate [Zgcd_alt] given here used to be the main [Zgcd]
18
function (see file [Znumtheory]), but this main [Zgcd] is now
19
based on a modern binary-efficient algorithm. This earlier
20
version, based on Euler's algorithm of iterated modulo, is kept
21
here due to both its intrinsic interest and its use as reference
22
point when proving gcd on Int31 numbers *)
24
Require Import ZArith_base.
25
Require Import ZArithRing.
27
Require Import Znumtheory.
31
(** In Coq, we need to control the number of iteration of modulo.
32
For that, we use an explicit measure in [nat], and we prove later
33
that using [2*d] is enough, where [d] is the number of binary
34
digits of the first argument. *)
36
Fixpoint Zgcdn (n:nat) : Z -> Z -> Z := fun a b =>
38
| O => 1 (* arbitrary, since n should be big enough *)
41
| Zpos _ => Zgcdn n (Zmod b a) a
42
| Zneg a => Zgcdn n (Zmod b (Zpos a)) (Zpos a)
46
Definition Zgcd_bound (a:Z) :=
49
| Zpos p => let n := Psize p in (n+n)%nat
50
| Zneg p => let n := Psize p in (n+n)%nat
53
Definition Zgcd_alt a b := Zgcdn (Zgcd_bound a) a b.
55
(** A first obvious fact : [Zgcd a b] is positive. *)
57
Lemma Zgcdn_pos : forall n a b,
61
simpl; auto with zarith.
62
destruct a; simpl; intros; auto with zarith; auto.
65
Lemma Zgcd_alt_pos : forall a b, 0 <= Zgcd_alt a b.
67
intros; unfold Zgcd; apply Zgcdn_pos; auto.
70
(** We now prove that Zgcd is indeed a gcd. *)
72
(** 1) We prove a weaker & easier bound. *)
74
Lemma Zgcdn_linear_bound : forall n a b,
75
Zabs a < Z_of_nat n -> Zis_gcd a b (Zgcdn n a b).
79
elimtype False; generalize (Zabs_pos a); omega.
80
destruct a; intros; simpl;
81
[ generalize (Zis_gcd_0_abs b); intuition | | ];
83
generalize (Z_div_mod b (Zpos p) (refl_equal Gt));
84
destruct (Zdiv_eucl b (Zpos p)) as (q,r);
86
rewrite inj_S in H; simpl Zabs in H;
87
(assert (H2: Zabs r < Z_of_nat n) by
88
(rewrite Zabs_eq; auto with zarith));
89
assert (IH:=IHn r (Zpos p) H2); clear IHn;
92
apply Zis_gcd_for_euclid2; auto.
93
apply Zis_gcd_minus; apply Zis_gcd_sym.
94
apply Zis_gcd_for_euclid2; auto.
97
(** 2) For Euclid's algorithm, the worst-case situation corresponds
98
to Fibonacci numbers. Let's define them: *)
100
Fixpoint fibonacci (n:nat) : Z :=
104
| S (S n as p) => fibonacci p + fibonacci n
107
Lemma fibonacci_pos : forall n, 0 <= fibonacci n.
109
cut (forall N n, (n<N)%nat -> 0<=fibonacci n).
115
simpl; auto with zarith.
117
simpl; auto with zarith.
118
change (0 <= fibonacci (S n) + fibonacci n).
119
generalize (IHN n) (IHN (S n)); omega.
122
Lemma fibonacci_incr :
123
forall n m, (n<=m)%nat -> fibonacci n <= fibonacci m.
127
apply Zle_trans with (fibonacci m); auto.
130
simpl; auto with zarith.
131
change (fibonacci (S m) <= fibonacci (S m)+fibonacci m).
132
generalize (fibonacci_pos m); omega.
135
(** 3) We prove that fibonacci numbers are indeed worst-case:
136
for a given number [n], if we reach a conclusion about [gcd(a,b)] in
137
exactly [n+1] loops, then [fibonacci (n+1)<=a /\ fibonacci(n+2)<=b] *)
139
Lemma Zgcdn_worst_is_fibonacci : forall n a b,
141
Zis_gcd a b (Zgcdn (S n) a b) ->
142
Zgcdn n a b <> Zgcdn (S n) a b ->
143
fibonacci (S n) <= a /\
144
fibonacci (S (S n)) <= b.
150
destruct a; [simpl in *; omega| | destruct H; discriminate].
151
revert H1; revert H0.
152
set (m:=S n) in *; (assert (m=S n) by auto); clearbody m.
153
pattern m at 2; rewrite H0.
155
unfold Zmod; generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
156
destruct (Zdiv_eucl b (Zpos p)) as (q,r).
159
destruct (Zle_lt_or_eq _ _ H2).
160
generalize (IHn _ _ (conj H4 H3)).
162
replace (fibonacci (S (S m))) with (fibonacci (S m) + fibonacci m) by auto.
163
assert (r = Zpos p * (-q) + b) by (rewrite H1; ring).
165
pattern r at 1; rewrite H8.
167
apply Zis_gcd_for_euclid2; auto.
168
apply Zis_gcd_sym; auto.
171
apply Zplus_le_compat; auto.
172
apply Zle_trans with (Zpos p * 1); auto.
173
ring_simplify (Zpos p * 1); auto.
174
apply Zmult_le_compat_l.
177
assert (0 < Zpos p0) by (compute; auto).
179
assert (Zpos p * Zneg p0 < 0) by (compute; auto).
181
compute; intros; discriminate.
196
(** 3b) We reformulate the previous result in a more positive way. *)
198
Lemma Zgcdn_ok_before_fibonacci : forall n a b,
199
0 < a < b -> a < fibonacci (S n) ->
200
Zis_gcd a b (Zgcdn n a b).
202
destruct a; [ destruct 1; elimtype False; omega | | destruct 1; discriminate].
204
k = (S (nat_of_P p) - n)%nat ->
205
0 < Zpos p < b -> Zpos p < fibonacci (S n) ->
206
Zis_gcd (Zpos p) b (Zgcdn n (Zpos p) b)).
208
clear n; induction k.
210
assert (nat_of_P p < n)%nat by omega.
211
apply Zgcdn_linear_bound.
213
generalize (inj_le _ _ H2).
215
rewrite <- Zpos_eq_Z_of_nat_o_nat_of_P; auto.
218
generalize (Zgcdn_worst_is_fibonacci n (Zpos p) b H0); intros.
219
assert (Zis_gcd (Zpos p) b (Zgcdn (S n) (Zpos p) b)).
222
replace (fibonacci (S (S n))) with (fibonacci (S n)+fibonacci n) by auto.
223
generalize (fibonacci_pos n); omega.
224
replace (Zgcdn n (Zpos p) b) with (Zgcdn (S n) (Zpos p) b); auto.
225
generalize (H2 H3); clear H2 H3; omega.
228
(** 4) The proposed bound leads to a fibonacci number that is big enough. *)
230
Lemma Zgcd_bound_fibonacci :
231
forall a, 0 < a -> a < fibonacci (Zgcd_bound a).
233
destruct a; [omega| | intro H; discriminate].
235
induction p; [ | | compute; auto ];
236
simpl Zgcd_bound in *;
237
rewrite plus_comm; simpl plus;
238
set (n:= (Psize p+Psize p)%nat) in *; simpl;
239
assert (n <> O) by (unfold n; destruct p; simpl; auto).
241
destruct n as [ |m]; [elim H; auto| ].
242
generalize (fibonacci_pos m); rewrite Zpos_xI; omega.
244
destruct n as [ |m]; [elim H; auto| ].
245
generalize (fibonacci_pos m); rewrite Zpos_xO; omega.
248
(* 5) the end: we glue everything together and take care of
249
situations not corresponding to [0<a<b]. *)
252
forall n a b, (Zgcd_bound a <= n)%nat ->
253
Zis_gcd a b (Zgcdn n a b).
257
destruct n; [elimtype False; omega | ].
258
simpl; generalize (Zis_gcd_0_abs b); intuition.
260
generalize (Zgcd_bound_fibonacci (Zpos p)).
261
simpl Zgcd_bound in *.
262
remember (Psize p+Psize p)%nat as m.
264
rewrite Heqm; destruct p; simpl; rewrite 1? plus_comm;
266
destruct m as [ |m]; [inversion H0; auto| ].
267
destruct n as [ |n]; [inversion H; auto| ].
270
generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
271
destruct (Zdiv_eucl b (Zpos p)) as (q,r).
274
apply Zis_gcd_for_euclid2.
276
destruct (Zle_lt_or_eq _ _ H1).
277
apply Zgcdn_ok_before_fibonacci; auto.
278
apply Zlt_le_trans with (fibonacci (S m)); [ omega | apply fibonacci_incr; auto].
280
destruct m as [ |m]; [elimtype False; omega| ].
281
destruct n as [ |n]; [elimtype False; omega| ].
282
simpl; apply Zis_gcd_sym; apply Zis_gcd_0.
284
generalize (Zgcd_bound_fibonacci (Zpos p)).
285
simpl Zgcd_bound in *.
286
remember (Psize p+Psize p)%nat as m.
288
rewrite Heqm; destruct p; simpl; rewrite 1? plus_comm;
290
destruct m as [ |m]; [inversion H0; auto| ].
291
destruct n as [ |n]; [inversion H; auto| ].
294
generalize (Z_div_mod b (Zpos p) (refl_equal Gt)).
295
destruct (Zdiv_eucl b (Zpos p)) as (q,r).
300
apply Zis_gcd_for_euclid2.
302
destruct (Zle_lt_or_eq _ _ H2).
303
apply Zgcdn_ok_before_fibonacci; auto.
304
apply Zlt_le_trans with (fibonacci (S m)); [ omega | apply fibonacci_incr; auto].
306
destruct m as [ |m]; [elimtype False; omega| ].
307
destruct n as [ |n]; [elimtype False; omega| ].
308
simpl; apply Zis_gcd_sym; apply Zis_gcd_0.
312
forall a b, Zis_gcd a b (Zgcd_alt a b).
314
unfold Zgcd_alt; intros; apply Zgcdn_is_gcd; auto.
added
removed
theories/ZArith/Zgcd_alt.v
