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
|
(* Here we give some functions that compute non-rational reals,
to measure the computation speed. *)
(* Expected time < 5.00s *)
Require Import QArith Qabs.
Require Import ConstructiveCauchyRealsMult.
Require Import Lqa.
Require Import Lia.
Local Open Scope CReal_scope.
(* We would need a shift instruction on positives to do this properly *)
Definition CReal_sqrt_Q_seq (q : Q) (n : Z) : Q
:= let (k,j) := q in
match k with
| Z0 => 0
| Z.pos i => match n with
| Z0
| Z.pos _ => Z.pos (Pos.sqrt (i*j)) # (j)
| Z.neg n' => Z.pos (Pos.sqrt (i*j*2^(2*n'))) # (j*2^n')
end
| Z.neg i => 0 (* unused *)
end.
Local Lemma Pos_pow_twice_r a b : (a^(2*b) = a^b * a^b)%positive.
Proof.
apply Pos2Z.inj.
rewrite Pos2Z.inj_mul.
do 2 rewrite Pos2Z.inj_pow.
rewrite Pos2Z.inj_mul.
apply Z.pow_twice_r.
Qed.
(* Approximation of the square root from below,
improves the convergence modulus. *)
Lemma CReal_sqrt_Q_le_below : forall (q : Q) (n : Z),
(0<=q)%Q -> (CReal_sqrt_Q_seq q n * CReal_sqrt_Q_seq q n <= q)%Q.
Proof.
intros q n Hqpos. destruct q as [k j]. unfold CReal_sqrt_Q_seq.
destruct k as [|i|i].
- apply Z.le_refl.
- destruct n as [|n|n].
+ pose proof (Pos.sqrt_spec (i * j)) as H. simpl in H.
destruct H as [H _].
unfold Qle, Qmult, Qnum, Qden.
rewrite <- Pos2Z.inj_mul, <- Pos2Z.inj_mul, <- Pos2Z.inj_mul.
apply Pos2Z.pos_le_pos. rewrite (Pos.mul_assoc i j j).
apply Pos.mul_le_mono_r; exact H.
+ pose proof (Pos.sqrt_spec (i * j)) as H. simpl in H.
destruct H as [H _].
unfold Qle, Qmult, Qnum, Qden.
rewrite <- Pos2Z.inj_mul, <- Pos2Z.inj_mul, <- Pos2Z.inj_mul.
apply Pos2Z.pos_le_pos. rewrite (Pos.mul_assoc i j j).
apply Pos.mul_le_mono_r; exact H.
+ pose proof (Pos.sqrt_spec (i * j * 2^(2*n))) as H. simpl in H.
destruct H as [H _].
unfold Qle, Qmult, Qnum, Qden.
rewrite <- Pos2Z.inj_mul, <- Pos2Z.inj_mul, <- Pos2Z.inj_mul.
apply Pos2Z.pos_le_pos. rewrite (Pos.mul_comm j (2^n)) at 2.
do 3 rewrite Pos.mul_assoc.
apply Pos.mul_le_mono_r.
simpl.
rewrite Pos_pow_twice_r in H at 3.
rewrite Pos.mul_assoc in H.
exact H.
- exact Hqpos.
Qed.
Lemma CReal_sqrt_Q_lt_above : forall (q : Q) (n : Z),
(0 <= q)%Q -> (q < ((CReal_sqrt_Q_seq q n + 2^n) * (CReal_sqrt_Q_seq q n + 2^n)))%Q.
Proof.
intros. destruct q as [k j]. unfold CReal_sqrt_Q_seq.
destruct k as [|i|i].
- ring_simplify.
setoid_rewrite <- Qpower.Qpower_mult.
setoid_rewrite QExtra.Qzero_eq.
pose proof QExtra.Qpower_pos_lt 2 (n*2)%Z ltac:(lra).
lra.
- destruct n as [|n|n].
+ pose proof (Pos.sqrt_spec (i * j)). simpl in H0.
destruct H0 as [_ H0].
change (2^0)%Q with 1%Q.
unfold Qlt, Qplus, Qmult, Qnum, Qden.
rewrite Pos.mul_1_r, Z.mul_1_r, Z.mul_1_l.
repeat rewrite <- Pos2Z.inj_add, <- Pos2Z.inj_mul.
apply Pos2Z.pos_lt_pos.
rewrite Pos.mul_assoc.
apply Pos.mul_lt_mono_r.
apply (Pos.lt_le_trans _ _ _ H0).
apply Pos.mul_le_mono; lia.
+ pose proof (Pos.sqrt_spec (i * j)). simpl in H0.
destruct H0 as [_ H0].
rewrite QExtra.Qpower_decomp'.
unfold Qlt, Qplus, Qmult, Qnum, Qden.
rewrite PosExtra.Pos_pow_1_r.
rewrite Pos.mul_1_r, Z.mul_1_r.
rewrite <- Pos2Z.inj_pow; do 2 rewrite <- Pos2Z.inj_mul; rewrite <- Pos2Z.inj_add.
apply Pos2Z.pos_lt_pos.
rewrite Pos.mul_assoc.
apply Pos.mul_lt_mono_r.
apply (Pos.lt_le_trans _ _ _ H0).
apply Pos.mul_le_mono;
pose proof Pos.le_1_l (2 ^ n * j)%positive; lia.
+ pose proof (Pos.sqrt_spec (i * j * 2 ^ (2 * n))). simpl in H0.
destruct H0 as [_ H0].
rewrite <- Pos2Z.opp_pos, Qpower.Qpower_opp.
rewrite QExtra.Qpower_decomp'.
rewrite <- Pos2Z.inj_pow, PosExtra.Pos_pow_1_r, <- QExtra.Qinv_swap_pos.
unfold Qlt, Qplus, Qmult, Qnum, Qden.
repeat rewrite Pos2Z.inj_mul.
ring_simplify.
replace (Z.pos i * Z.pos j ^ 2 * Z.pos (2 ^ n) ^ 4)%Z
with ((Z.pos i * Z.pos j * Z.pos (2 ^ n) ^ 2) * (Z.pos j * Z.pos (2 ^ n) ^ 2))%Z by ring.
replace (
Z.pos j ^ 3 * Z.pos (2 ^ n) ^ 2 +
2 * Z.pos j ^ 2 * Z.pos (2 ^ n) ^ 2 * Z.pos (Pos.sqrt (i * j * 2 ^ (2 * n))) +
Z.pos j * Z.pos (2 ^ n) ^ 2 * Z.pos (Pos.sqrt (i * j * 2 ^ (2 * n))) ^ 2)%Z
with (
(Z.pos j + Z.pos (Pos.sqrt (i * j * 2 ^ (2 * n))))^2 *
(Z.pos j * Z.pos (2 ^ n) ^ 2))%Z by ring.
repeat rewrite Pos2Z.inj_pow.
rewrite <- Z.pow_mul_r by lia.
repeat rewrite <- Pos2Z.inj_mul.
repeat rewrite <- Pos2Z.inj_pow.
repeat rewrite <- Pos2Z.inj_mul.
repeat rewrite <- Pos2Z.inj_add.
apply Pos2Z.pos_lt_pos.
rewrite (Pos.mul_comm n 2); change (2*n)%positive with (n~0)%positive.
apply Pos.mul_lt_mono_r.
apply (Pos.lt_le_trans _ _ _ H0).
apply Pos.mul_le_mono;
pose proof Pos.le_1_l (2 ^ n * j)%positive; lia.
- exfalso; unfold Qle, Z.le in H; simpl in H; exact (H eq_refl).
Qed.
Lemma CReal_sqrt_Q_pos : forall (q : Q) (n : Z),
(0 <= (CReal_sqrt_Q_seq q n))%Q.
Proof.
intros. unfold CReal_sqrt_Q_seq. destruct q, Qnum.
- apply Qle_refl.
- destruct n as [|n|n]; discriminate.
- apply Qle_refl.
Qed.
Lemma Qsqrt_lt : forall q r :Q,
(0 <= r -> q*q < r*r -> q < r)%Q.
Proof.
intros. destruct (Q_dec q r). destruct s. exact q0.
- exfalso. apply (Qlt_not_le _ _ H0). apply (Qle_trans _ (q * r)).
apply Qmult_le_compat_r. apply Qlt_le_weak, q0. exact H.
rewrite Qmult_comm.
apply Qmult_le_compat_r. apply Qlt_le_weak, q0.
apply (Qle_trans _ r _ H). apply Qlt_le_weak, q0.
- exfalso. rewrite q0 in H0. exact (Qlt_irrefl _ H0).
Qed.
Lemma CReal_sqrt_Q_cauchy :
forall q:Q, QCauchySeq (CReal_sqrt_Q_seq q).
Proof.
intro q. destruct q as [k j]. destruct k.
- intros n a b H H0.
change (Qabs _) with 0%Q.
apply QExtra.Qpower_pos_lt; reflexivity.
- assert (forall n a b, (b<=n)%Z ->
(CReal_sqrt_Q_seq (Z.pos p # j) a - CReal_sqrt_Q_seq (Z.pos p # j) b
< 2^n)%Q).
{ intros.
pose proof QExtra.Qpower_pos_lt 2 n eq_refl as Hpow.
rewrite <- (Qplus_lt_r _ _ (CReal_sqrt_Q_seq (Z.pos p # j) b)).
ring_simplify. apply Qsqrt_lt.
{ apply (Qle_trans _ (0+2^n)). lra.
apply Qplus_le_l. apply CReal_sqrt_Q_pos. }
apply (Qle_lt_trans _ (Z.pos p # j)).
{ apply CReal_sqrt_Q_le_below. discriminate. }
apply (Qlt_le_trans _ ((CReal_sqrt_Q_seq (Z.pos p # j) b + (2^b)) *
(CReal_sqrt_Q_seq (Z.pos p # j) b + (2^b)))).
{ apply CReal_sqrt_Q_lt_above. discriminate. }
apply (Qle_trans _ ((CReal_sqrt_Q_seq (Z.pos p # j) b + (2^n)) *
(CReal_sqrt_Q_seq (Z.pos p # j) b + (2^b)))).
{ apply Qmult_le_r.
- apply (Qlt_le_trans _ (0+(2^b))).
+ rewrite Qplus_0_l. apply QExtra.Qpower_pos_lt. reflexivity.
+ apply Qplus_le_l. apply CReal_sqrt_Q_pos.
- apply Qplus_le_r. apply QExtra.Qpower_le_compat.
exact H. discriminate. }
apply QExtra.Qmult_le_compat_nonneg.
- split.
+ pose proof CReal_sqrt_Q_pos (Z.pos p # j) b.
lra.
+ apply Qle_refl.
- split.
+ pose proof CReal_sqrt_Q_pos (Z.pos p # j) b.
pose proof QExtra.Qpower_pos_lt 2 b eq_refl as Hpowb.
lra.
+ apply Qplus_le_r.
apply QExtra.Qpower_le_compat.
exact H. discriminate.
}
intros n a b H0 H1. apply Qabs_case.
intros. apply H, H1.
intros.
setoid_replace (- (CReal_sqrt_Q_seq (Z.pos p # j) a - CReal_sqrt_Q_seq (Z.pos p # j) b))%Q
with (CReal_sqrt_Q_seq (Z.pos p # j) b - CReal_sqrt_Q_seq (Z.pos p # j) a)%Q.
2: ring. apply H, H0.
- intros n a b H H0.
change (Qabs _) with 0%Q.
apply QExtra.Qpower_pos_lt; reflexivity.
Qed.
Definition CReal_sqrt_Q_scale (q : Q) : Z
:= ((QExtra.Qbound_lt_ZExp2 q + 1)/2)%Z.
Lemma CReal_sqrt_Q_bound : forall (q : Q),
QBound (CReal_sqrt_Q_seq q) (CReal_sqrt_Q_scale q).
Proof.
intros q k.
unfold CReal_sqrt_Q_scale.
rewrite Qabs_pos.
2: apply CReal_sqrt_Q_pos.
apply Qsqrt_lt.
1: apply Qpower.Qpower_pos; discriminate.
destruct (Qlt_le_dec q 0) as [Hq|Hq].
- destruct q as [[|n|n] d].
+ discriminate Hq.
+ discriminate Hq.
+ reflexivity.
- apply (Qle_lt_trans _ _ _ (CReal_sqrt_Q_le_below _ _ Hq)).
rewrite <- Qpower.Qpower_plus.
2: discriminate.
rewrite Z.add_diag, Z.mul_comm.
pose proof Zdiv.Zmod_eq (QExtra.Qbound_lt_ZExp2 q + 1) 2 eq_refl as Hmod.
assert (forall a b c : Z, c=b-a -> a=b-c)%Z as H by (intros a b c H'; rewrite H'; ring).
apply H in Hmod; rewrite Hmod; clear H Hmod.
apply (Qlt_le_trans _ _ _ (QExtra.Qbound_lt_ZExp2_spec q)).
apply QExtra.Qpower_le_compat. 2: discriminate.
pose proof Z.mod_pos_bound (QExtra.Qbound_lt_ZExp2 q + 1)%Z 2%Z eq_refl.
lia.
Qed.
Definition CReal_sqrt_Q (q : Q) : CReal :=
{|
seq := CReal_sqrt_Q_seq q;
scale := CReal_sqrt_Q_scale q;
cauchy := CReal_sqrt_Q_cauchy q;
bound := CReal_sqrt_Q_bound q
|}.
Time Eval vm_compute in (seq (CReal_sqrt_Q 2) (-1000)%Z).
|
