1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
|
(************************************************************************)
(* v * The Coq Proof Assistant / The Coq Development Team *)
(* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2017 *)
(* \VV/ **************************************************************)
(* // * This file is distributed under the terms of the *)
(* * GNU Lesser General Public License Version 2.1 *)
(************************************************************************)
(** * Zgcd_alt : an alternate version of Z.gcd, based on Euclid's algorithm *)
(**
Author: Pierre Letouzey
*)
(** The alternate [Zgcd_alt] given here used to be the main [Z.gcd]
function (see file [Znumtheory]), but this main [Z.gcd] is now
based on a modern binary-efficient algorithm. This earlier
version, based on Euclid's algorithm of iterated modulo, is kept
here due to both its intrinsic interest and its use as reference
point when proving gcd on Int31 numbers *)
Require Import ZArith_base.
Require Import ZArithRing.
Require Import Zdiv.
Require Import Znumtheory.
Require Import Omega.
Open Scope Z_scope.
(** In Coq, we need to control the number of iteration of modulo.
For that, we use an explicit measure in [nat], and we prove later
that using [2*d] 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 => Z.abs b
| Zpos _ => Zgcdn n (Z.modulo b a) a
| Zneg a => Zgcdn n (Z.modulo b (Zpos a)) (Zpos a)
end
end.
Definition Zgcd_bound (a:Z) :=
match a with
| Z0 => S O
| Zpos p => let n := Pos.size_nat p in (n+n)%nat
| Zneg p => let n := Pos.size_nat p in (n+n)%nat
end.
Definition Zgcd_alt a b := Zgcdn (Zgcd_bound a) a b.
(** A first obvious fact : [Z.gcd a b] is positive. *)
Lemma Zgcdn_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_alt_pos : forall a b, 0 <= Zgcd_alt a b.
Proof.
intros; unfold Z.gcd; apply Zgcdn_pos; auto.
Qed.
(** We now prove that Z.gcd is indeed a gcd. *)
(** 1) We prove a weaker & easier bound. *)
Lemma Zgcdn_linear_bound : forall n a b,
Z.abs a < Z.of_nat n -> Zis_gcd a b (Zgcdn n a b).
Proof.
induction n.
simpl; intros.
exfalso; generalize (Z.abs_nonneg a); omega.
destruct a; intros; simpl;
[ generalize (Zis_gcd_0_abs b); intuition | | ];
unfold Z.modulo;
generalize (Z_div_mod b (Zpos p) (eq_refl Gt));
destruct (Z.div_eucl b (Zpos p)) as (q,r);
intros (H0,H1);
rewrite Nat2Z.inj_succ in H; simpl Z.abs in H;
(assert (H2: Z.abs r < Z.of_nat n) by
(rewrite Z.abs_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.
enough (forall N n, (n<N)%nat -> 0<=fibonacci n) by 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 Z.le_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.
intros [|a|a]; intros; simpl; omega.
intros [|a|a] b (Ha,Ha'); [simpl; omega | | easy ].
remember (S n) as m.
rewrite Heqm at 2. simpl Zgcdn.
unfold Z.modulo; generalize (Z_div_mod b (Zpos a) eq_refl).
destruct (Z.div_eucl b (Zpos a)) as (q,r).
intros (EQ,(Hr,Hr')).
Z.le_elim Hr.
- (* r > 0 *)
replace (fibonacci (S (S m))) with (fibonacci (S m) + fibonacci m) by auto.
intros.
destruct (IHn r (Zpos a) (conj Hr Hr')); auto.
+ assert (EQ' : r = Zpos a * (-q) + b) by (rewrite EQ; ring).
rewrite EQ' at 1.
apply Zis_gcd_sym.
apply Zis_gcd_for_euclid2; auto.
apply Zis_gcd_sym; auto.
+ split; auto.
rewrite EQ.
apply Z.add_le_mono; auto.
apply Z.le_trans with (Zpos a * 1); auto.
now rewrite Z.mul_1_r.
apply Z.mul_le_mono_nonneg_l; auto with zarith.
change 1 with (Z.succ 0). apply Z.le_succ_l.
destruct q; auto with zarith.
assert (Zpos a * Zneg p < 0) by now compute. omega.
- (* r = 0 *)
clear IHn EQ Hr'; intros _.
subst r; simpl; rewrite Heqm.
destruct n.
+ simpl. omega.
+ now destruct 1.
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; exfalso; omega | | destruct 1; discriminate].
cut (forall k n b,
k = (S (Pos.to_nat 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 (Pos.to_nat p < n)%nat by omega.
apply Zgcdn_linear_bound.
simpl.
generalize (inj_le _ _ H2).
rewrite Nat2Z.inj_succ.
rewrite positive_nat_Z; 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; [ | | compute; auto ];
simpl Zgcd_bound in *;
rewrite plus_comm; simpl plus;
set (n:= (Pos.size_nat p+Pos.size_nat p)%nat) in *; simpl;
assert (n <> O) by (unfold n; destruct p; simpl; auto).
destruct n as [ |m]; [elim H; auto| ].
generalize (fibonacci_pos m); rewrite Pos2Z.inj_xI; omega.
destruct n as [ |m]; [elim H; auto| ].
generalize (fibonacci_pos m); rewrite Pos2Z.inj_xO; omega.
Qed.
(* 5) the end: we glue everything together and take care of
situations not corresponding to [0<a<b]. *)
Lemma Zgcd_bound_opp a : Zgcd_bound (-a) = Zgcd_bound a.
Proof.
now destruct a.
Qed.
Lemma Zgcdn_opp n a b : Zgcdn n (-a) b = Zgcdn n a b.
Proof.
induction n; simpl; auto.
destruct a; simpl; auto.
Qed.
Lemma Zgcdn_is_gcd_pos n a b : (Zgcd_bound (Zpos a) <= n)%nat ->
Zis_gcd (Zpos a) b (Zgcdn n (Zpos a) b).
Proof.
intros.
generalize (Zgcd_bound_fibonacci (Zpos a)).
simpl Zgcd_bound in *.
remember (Pos.size_nat a+Pos.size_nat a)%nat as m.
assert (1 < m)%nat.
{ rewrite Heqm; destruct a; simpl; rewrite 1?plus_comm;
auto with arith. }
destruct m as [ |m]; [inversion H0; auto| ].
destruct n as [ |n]; [inversion H; auto| ].
simpl Zgcdn.
unfold Z.modulo.
generalize (Z_div_mod b (Zpos a) (eq_refl Gt)).
destruct (Z.div_eucl b (Zpos a)) as (q,r).
intros (->,(H1,H2)) H3.
apply Zis_gcd_for_euclid2.
Z.le_elim H1.
+ apply Zgcdn_ok_before_fibonacci; auto.
apply Z.lt_le_trans with (fibonacci (S m));
[ omega | apply fibonacci_incr; auto].
+ subst r; simpl.
destruct m as [ |m]; [exfalso; omega| ].
destruct n as [ |n]; [exfalso; omega| ].
simpl; apply Zis_gcd_sym; apply Zis_gcd_0.
Qed.
Lemma Zgcdn_is_gcd n a b :
(Zgcd_bound a <= n)%nat -> Zis_gcd a b (Zgcdn n a b).
Proof.
destruct a.
- simpl; intros.
destruct n; [exfalso; omega | ].
simpl; generalize (Zis_gcd_0_abs b); intuition.
- apply Zgcdn_is_gcd_pos.
- rewrite <- Zgcd_bound_opp, <- Zgcdn_opp.
intros. apply Zis_gcd_minus, Zis_gcd_sym. simpl Z.opp.
now apply Zgcdn_is_gcd_pos.
Qed.
Lemma Zgcd_is_gcd :
forall a b, Zis_gcd a b (Zgcd_alt a b).
Proof.
unfold Zgcd_alt; intros; apply Zgcdn_is_gcd; auto.
Qed.
|
