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
|
(**
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2010-2018 Sylvie Boldo
#<br />#
Copyright (C) 2010-2018 Guillaume Melquiond
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Helper functions and theorems for computing the rounded square root of a floating-point number. *)
Require Import Raux Defs Digits Generic_fmt Float_prop Bracket.
Set Implicit Arguments.
Set Strongly Strict Implicit.
Section Fcalc_sqrt.
Variable beta : radix.
Notation bpow e := (bpow beta e).
Variable fexp : Z -> Z.
(** Computes a mantissa of precision p, the corresponding exponent,
and the position with respect to the real square root of the
input floating-point number.
The algorithm performs the following steps:
- Shift the mantissa so that it has at least 2p-1 digits;
shift it one digit more if the new exponent is not even.
- Compute the square root s (at least p digits) of the new
mantissa, and its remainder r.
- Compute the position according to the remainder:
-- r == 0 => Eq,
-- r <= s => Lo,
-- r >= s => Up.
Complexity is fine as long as p1 <= 2p-1.
*)
Lemma mag_sqrt_F2R :
forall m1 e1,
(0 < m1)%Z ->
mag beta (sqrt (F2R (Float beta m1 e1))) = Z.div2 (Zdigits beta m1 + e1 + 1) :> Z.
Proof.
intros m1 e1 Hm1.
rewrite <- (mag_F2R_Zdigits beta m1 e1) by now apply Zgt_not_eq.
apply mag_sqrt.
now apply F2R_gt_0.
Qed.
Definition Fsqrt_core m1 e1 e :=
let d1 := Zdigits beta m1 in
let m1' := (m1 * Zpower beta (e1 - 2 * e))%Z in
let (q, r) := Z.sqrtrem m1' in
let l :=
if Zeq_bool r 0 then loc_Exact
else loc_Inexact (if Zle_bool r q then Lt else Gt) in
(q, l).
Theorem Fsqrt_core_correct :
forall m1 e1 e,
(0 < m1)%Z ->
(2 * e <= e1)%Z ->
let '(m, l) := Fsqrt_core m1 e1 e in
inbetween_float beta m e (sqrt (F2R (Float beta m1 e1))) l.
Proof.
intros m1 e1 e Hm1 He.
unfold Fsqrt_core.
set (m' := Zmult _ _).
assert (0 <= m')%Z as Hm'.
{ apply Z.mul_nonneg_nonneg.
now apply Zlt_le_weak.
apply Zpower_ge_0. }
assert (sqrt (F2R (Float beta m1 e1)) = sqrt (IZR m') * bpow e)%R as Hf.
{ rewrite <- (sqrt_Rsqr (bpow e)) by apply bpow_ge_0.
rewrite <- sqrt_mult.
unfold Rsqr, m'.
rewrite mult_IZR, IZR_Zpower by omega.
rewrite Rmult_assoc, <- 2!bpow_plus.
now replace (_ + _)%Z with e1 by ring.
now apply IZR_le.
apply Rle_0_sqr. }
generalize (Z.sqrtrem_spec m' Hm').
destruct Z.sqrtrem as [q r].
intros [Hq Hr].
rewrite Hf.
unfold inbetween_float, F2R. simpl Fnum.
apply inbetween_mult_compat.
apply bpow_gt_0.
rewrite Hq.
case Zeq_bool_spec ; intros Hr'.
(* .. r = 0 *)
rewrite Hr', Zplus_0_r, mult_IZR.
fold (Rsqr (IZR q)).
rewrite sqrt_Rsqr.
now constructor.
apply IZR_le.
clear -Hr ; omega.
(* .. r <> 0 *)
constructor.
split.
(* ... bounds *)
apply Rle_lt_trans with (sqrt (IZR (q * q))).
rewrite mult_IZR.
fold (Rsqr (IZR q)).
rewrite sqrt_Rsqr.
apply Rle_refl.
apply IZR_le.
clear -Hr ; omega.
apply sqrt_lt_1.
rewrite mult_IZR.
apply Rle_0_sqr.
rewrite <- Hq.
now apply IZR_le.
apply IZR_lt.
omega.
apply Rlt_le_trans with (sqrt (IZR ((q + 1) * (q + 1)))).
apply sqrt_lt_1.
rewrite <- Hq.
now apply IZR_le.
rewrite mult_IZR.
apply Rle_0_sqr.
apply IZR_lt.
ring_simplify.
omega.
rewrite mult_IZR.
fold (Rsqr (IZR (q + 1))).
rewrite sqrt_Rsqr.
apply Rle_refl.
apply IZR_le.
clear -Hr ; omega.
(* ... location *)
rewrite Rcompare_half_r.
generalize (Rcompare_sqr (2 * sqrt (IZR (q * q + r))) (IZR q + IZR (q + 1))).
rewrite 2!Rabs_pos_eq.
intros <-.
replace ((2 * sqrt (IZR (q * q + r))) * (2 * sqrt (IZR (q * q + r))))%R
with (4 * Rsqr (sqrt (IZR (q * q + r))))%R by (unfold Rsqr ; ring).
rewrite Rsqr_sqrt.
rewrite <- plus_IZR, <- 2!mult_IZR.
rewrite Rcompare_IZR.
replace ((q + (q + 1)) * (q + (q + 1)))%Z with (4 * (q * q) + 4 * q + 1)%Z by ring.
generalize (Zle_cases r q).
case (Zle_bool r q) ; intros Hr''.
change (4 * (q * q + r) < 4 * (q * q) + 4 * q + 1)%Z.
omega.
change (4 * (q * q + r) > 4 * (q * q) + 4 * q + 1)%Z.
omega.
rewrite <- Hq.
now apply IZR_le.
rewrite <- plus_IZR.
apply IZR_le.
clear -Hr ; omega.
apply Rmult_le_pos.
now apply IZR_le.
apply sqrt_ge_0.
Qed.
Definition Fsqrt (x : float beta) :=
let (m1, e1) := x in
let e' := (Zdigits beta m1 + e1 + 1)%Z in
let e := Z.min (fexp (Z.div2 e')) (Z.div2 e1) in
let '(m, l) := Fsqrt_core m1 e1 e in
(m, e, l).
Theorem Fsqrt_correct :
forall x,
(0 < F2R x)%R ->
let '(m, e, l) := Fsqrt x in
(e <= cexp beta fexp (sqrt (F2R x)))%Z /\
inbetween_float beta m e (sqrt (F2R x)) l.
Proof.
intros [m1 e1] Hm1.
apply gt_0_F2R in Hm1.
unfold Fsqrt.
set (e := Z.min _ _).
assert (2 * e <= e1)%Z as He.
{ assert (e <= Z.div2 e1)%Z by apply Z.le_min_r.
rewrite (Zdiv2_odd_eqn e1).
destruct Z.odd ; omega. }
generalize (Fsqrt_core_correct m1 e1 e Hm1 He).
destruct Fsqrt_core as [m l].
apply conj.
apply Z.le_trans with (1 := Z.le_min_l _ _).
unfold cexp.
rewrite (mag_sqrt_F2R m1 e1 Hm1).
apply Z.le_refl.
Qed.
End Fcalc_sqrt.
|
