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
|
(************************************************************************)
(* v * The Coq Proof Assistant / The Coq Development Team *)
(* <O___,, * INRIA - CNRS - LIX - LRI - PPS - Copyright 1999-2010 *)
(* \VV/ **************************************************************)
(* // * This file is distributed under the terms of the *)
(* * GNU Lesser General Public License Version 2.1 *)
(************************************************************************)
Require Import BinPos.
Open Scope positive_scope.
Definition Pleb x y :=
match Pcompare x y Eq with Gt => false | _ => true end.
(** A Square Root function for positive numbers *)
(** We procede by blocks of two digits : if p is written qbb'
then sqrt(p) will be sqrt(q)~0 or sqrt(q)~1.
For deciding easily in which case we are, we store the remainder
(as a positive_mask, since it can be null).
Instead of copy-pasting the following code four times, we
factorize as an auxiliary function, with f and g being either
xO or xI depending of the initial digits.
NB: (Pminus_mask (g (f 1)) 4) is a hack, morally it's g (f 0).
*)
Definition Psqrtrem_step (f g:positive->positive) p :=
match p with
| (s, IsPos r) =>
let s' := s~0~1 in
let r' := g (f r) in
if Pleb s' r' then (s~1, Pminus_mask r' s')
else (s~0, IsPos r')
| (s,_) => (s~0, Pminus_mask (g (f 1)) 4)
end.
Fixpoint Psqrtrem p : positive * positive_mask :=
match p with
| 1 => (1,IsNul)
| 2 => (1,IsPos 1)
| 3 => (1,IsPos 2)
| p~0~0 => Psqrtrem_step xO xO (Psqrtrem p)
| p~0~1 => Psqrtrem_step xO xI (Psqrtrem p)
| p~1~0 => Psqrtrem_step xI xO (Psqrtrem p)
| p~1~1 => Psqrtrem_step xI xI (Psqrtrem p)
end.
Definition Psqrt p := fst (Psqrtrem p).
(** An inductive relation for specifying sqrt results *)
Inductive PSQRT : positive*positive_mask -> positive -> Prop :=
| PSQRT_exact : forall s x, x=s*s -> PSQRT (s,IsNul) x
| PSQRT_approx : forall s r x, x=s*s+r -> r <= s~0 -> PSQRT (s,IsPos r) x.
(** Correctness proofs *)
Lemma Psqrtrem_step_spec : forall f g p x,
(f=xO \/ f=xI) -> (g=xO \/ g=xI) ->
PSQRT p x -> PSQRT (Psqrtrem_step f g p) (g (f x)).
Proof.
intros f g _ _ Hf Hg [ s _ -> | s r _ -> Hr ].
(* exact *)
unfold Psqrtrem_step.
destruct Hf,Hg; subst; simpl Pminus_mask;
constructor; try discriminate; now rewrite Psquare_xO.
(* approx *)
assert (Hfg : forall p q, g (f (p+q)) = p~0~0 + g (f q))
by (intros; destruct Hf, Hg; now subst).
unfold Psqrtrem_step. unfold Pleb.
case Pcompare_spec; [intros EQ | intros LT | intros GT].
(* - EQ *)
rewrite <- EQ. rewrite Pminus_mask_diag.
destruct Hg; subst g; try discriminate.
destruct Hf; subst f; try discriminate.
injection EQ; clear EQ; intros <-.
constructor. now rewrite Psquare_xI.
(* - LT *)
destruct (Pminus_mask_Gt (g (f r)) (s~0~1)) as (y & -> & H & _).
change Eq with (CompOpp Eq). now rewrite <- Pcompare_antisym, LT.
constructor.
rewrite Hfg, <- H. now rewrite Psquare_xI, Pplus_assoc.
apply Ple_lteq, Pcompare_p_Sq in Hr; change (r < s~1) in Hr.
apply Ple_lteq, Pcompare_p_Sq; change (y < s~1~1).
apply Pplus_lt_mono_l with (s~0~1).
rewrite H. simpl. rewrite Pplus_carry_spec, Pplus_diag. simpl.
set (s1:=s~1) in *; clearbody s1.
destruct Hf,Hg; subst; red; simpl;
apply Hr || apply Pcompare_eq_Lt, Hr.
(* - GT *)
constructor.
rewrite Hfg. now rewrite Psquare_xO.
apply Ple_lteq, Pcompare_p_Sq, GT.
Qed.
Lemma Psqrtrem_spec : forall p, PSQRT (Psqrtrem p) p.
Proof.
fix 1.
destruct p; try destruct p; try (constructor; easy);
apply Psqrtrem_step_spec; auto.
Qed.
Lemma Psqrt_spec : forall p,
let s := Psqrt p in s*s <= p < (Psucc s)*(Psucc s).
Proof.
intros p. simpl.
assert (H:=Psqrtrem_spec p).
unfold Psqrt in *. destruct Psqrtrem as (s,rm); simpl.
inversion_clear H; subst.
(* exact *)
split. red. rewrite Pcompare_refl. discriminate.
apply Pmult_lt_mono; apply Pcompare_p_Sp.
(* approx *)
split.
apply Ple_lteq; left. apply Plt_plus_r.
rewrite (Pplus_one_succ_l).
rewrite Pmult_plus_distr_r, !Pmult_plus_distr_l.
rewrite !Pmult_1_r. simpl (1*s).
rewrite <- Pplus_assoc, (Pplus_assoc s s), Pplus_diag, Pplus_assoc.
rewrite (Pplus_comm (_+_)). apply Pplus_lt_mono_l.
rewrite <- Pplus_one_succ_l. now apply Pcompare_p_Sq, Ple_lteq.
Qed.
|
