Sqrt_sqr.v 57 KB
Newer Older
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/

Copyright (C) 2013-2018 Sylvie Boldo

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.
*)

18 19
Require Import Reals Psatz.
Require Import Flocq.Core.Core.
20
Require Import Interval.Interval_tactic.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
21 22

Section Sec1.
23

BOLDO Sylvie's avatar
BOLDO Sylvie committed
24 25 26 27 28 29 30
Open Scope R_scope.

Variable beta : radix.
Notation bpow e := (bpow beta e).

Variable prec : Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
31 32 33
Variable choice1 : Z -> bool.
Variable choice2 : Z -> bool.
Variable choice3 : Z -> bool.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
34 35

Notation format := (generic_format beta (FLX_exp prec)).
36 37 38
Notation round_flx1 :=(round beta (FLX_exp prec) (Znearest choice1)).
Notation round_flx2 :=(round beta (FLX_exp prec) (Znearest choice2)).
Notation round_flx3 :=(round beta (FLX_exp prec) (Znearest choice3)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
39 40 41 42 43 44 45 46 47
Notation ulp_flx :=(ulp beta (FLX_exp prec)).
Notation pred_flx := (pred beta (FLX_exp prec)).

Hypothesis pGt1: (1 < prec)%Z.

Variables x:R.
Hypothesis xPos: 0 < x.
Hypothesis Fx: format x.

48 49
Let y:=round_flx1(x*x).
Let z:=round_flx2(sqrt y).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
50 51


Guillaume Melquiond's avatar
Guillaume Melquiond committed
52
Theorem round_flx_sqr_sqrt_middle: x = sqrt (IZR beta) * bpow (mag beta x - 1) -> z=x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
53 54 55
Proof with auto with typeclass_instances.
intros Hx.
unfold z, y.
56 57
replace (round_flx1 (x * x)) with (bpow (2*mag beta x -1)).
replace (sqrt (bpow (2 * mag beta x - 1))) with x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
58 59 60 61
apply round_generic...
apply sym_eq, sqrt_lem_1.
apply bpow_ge_0.
now left.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
62
rewrite Hx at 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
63
rewrite Rmult_comm; rewrite Hx at 1.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
64
apply trans_eq with ((sqrt (IZR beta)*sqrt (IZR beta)) * (bpow (mag beta x-1)*bpow (mag beta x-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
65 66 67 68
ring.
rewrite <- bpow_plus, sqrt_def, <- bpow_1, <- bpow_plus.
apply f_equal; ring.
left; apply radix_pos.
69
replace (x*x) with (bpow (2 * mag beta x - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
70 71 72
apply sym_eq, round_generic...
apply generic_format_bpow.
unfold FLX_exp; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
73
apply sym_eq; rewrite Hx at 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
74
rewrite Rmult_comm; rewrite Hx at 1.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
75
apply trans_eq with ((sqrt (IZR beta)*sqrt (IZR beta)) * (bpow (mag beta x-1)*bpow (mag beta x-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
76 77 78 79 80 81 82 83 84
ring.
rewrite <- bpow_plus, sqrt_def, <- bpow_1, <- bpow_plus.
apply f_equal; ring.
left; apply radix_pos.
Qed.

Theorem round_flx_sqr_sqrt_le: (beta <= 4)%Z -> z <= x.
Proof with auto with typeclass_instances.
intros Hb.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
85
case (Req_dec x (sqrt (IZR beta) * bpow (mag beta x - 1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
86 87 88 89 90 91 92
intros L; right; now apply round_flx_sqr_sqrt_middle.
intros Hx'.
assert (0 < x*x) by now apply Rmult_lt_0_compat.
assert (0 <= 1 + ulp_flx (x * x) / 2 / (x * x)).
rewrite <- (Rplus_0_l 0).
apply Rplus_le_compat; auto with real.
unfold Rdiv; apply Rmult_le_pos.
93
apply Rmult_le_pos; auto with real; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
94 95 96 97 98
left; auto with real.
assert (0 <= 1 + ulp_flx x / 2 / x).
rewrite <- (Rplus_0_l 0).
apply Rplus_le_compat; auto with real.
unfold Rdiv; apply Rmult_le_pos.
99
apply Rmult_le_pos; auto with real; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
100 101
left; auto with real.
(* *)
102
apply round_N_le_midp; try assumption...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
103
apply Rlt_le_trans with (x*(1+ulp_flx x / 2 / x)).
104
2: right; rewrite succ_eq_pos; try now left.
105
2: field; auto with real.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
106 107 108 109 110 111
apply Rle_lt_trans with (sqrt ((x*x)*(1+ulp_flx (x*x)/2/(x*x)))).
apply sqrt_le_1_alt.
apply Rplus_le_reg_l with (-x*x).
apply Rle_trans with (y-x*x);[right; ring|idtac].
apply Rle_trans with (/2*ulp_flx (x*x));[idtac|right; field; auto with real].
apply Rle_trans with (1:=RRle_abs _).
112
apply error_le_half_ulp...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
rewrite sqrt_mult; auto with real.
rewrite sqrt_square; auto with real.
apply Rmult_lt_compat_l.
assumption.
rewrite <- (sqrt_square (1 + ulp_flx x / 2 / x)); try assumption.
apply sqrt_lt_1_alt.
split; try assumption.
apply Rmult_lt_reg_l with (2*x*x).
rewrite Rmult_assoc; apply Rmult_lt_0_compat; auto with real.
apply Rplus_lt_reg_r with (-2*x*x - ulp_flx x*ulp_flx x/2).
apply Rle_lt_trans with (ulp_flx (x*x) - ulp_flx x*ulp_flx x/2).
right; field.
auto with real.
apply Rlt_le_trans with (2*x*ulp_flx x).
2: right; field; auto with real.
apply Rle_lt_trans with (ulp_flx (x * x) + ulp_flx x * ulp_flx x / 2).
apply Rplus_le_compat_l.
apply Rplus_le_reg_r with (ulp_flx x * ulp_flx x / 2).
apply Rle_trans with 0;[right; ring|idtac].
apply Rle_trans with (ulp_flx x * ulp_flx x).
133
apply Rmult_le_pos; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
134
right; field.
135 136 137
rewrite 2!ulp_neq_0.
2: now apply Rgt_not_eq.
2: now apply Rgt_not_eq.
138
unfold cexp, FLX_exp.
139
destruct (mag beta x) as (e,He).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
140 141 142 143 144 145
simpl.
assert (bpow (e - 1) <= x < bpow e).
rewrite <- (Rabs_right x).
apply He; auto with real.
apply Rle_ge; auto with real.
rewrite <- bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
146
case (Rle_or_lt x (sqrt (IZR (radix_val beta))*bpow (e-1))); intros Hx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
147 148
(* . *)
destruct Hx.
149
replace (mag beta (x * x)-prec)%Z with (2*e-1-prec)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
apply Rlt_le_trans with (bpow (2 * e - 1 - prec) + bpow (2*e - 1 - prec) / 2).
apply Rplus_lt_compat_l.
unfold Rdiv; apply Rmult_lt_compat_r.
auto with real.
apply bpow_lt.
omega.
apply Rle_trans with (2*bpow (2 * e - 1 - prec)).
apply Rle_trans with (3/2*bpow (2 * e - 1 - prec)).
right; field.
apply Rmult_le_compat_r.
apply bpow_ge_0.
interval.
rewrite Rmult_assoc; apply Rmult_le_compat_l.
auto with real.
apply Rle_trans with (bpow (e-1)*bpow (e - prec)).
rewrite <- bpow_plus.
apply bpow_le.
omega.
apply Rmult_le_compat_r.
apply bpow_ge_0.
apply H2.
apply f_equal with (f:= fun z => (z-prec)%Z).
172
apply sym_eq, mag_unique.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
173 174 175 176 177 178
rewrite Rabs_right.
split.
apply Rle_trans with (bpow (e-1)*bpow (e-1)).
rewrite <- bpow_plus.
right; apply f_equal; ring.
apply Rmult_le_compat; try apply H2; try apply bpow_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
179
apply Rlt_le_trans with ((sqrt (IZR beta)*bpow (e-1))*(sqrt (IZR beta)*bpow (e-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
180
apply Rmult_le_0_lt_compat; try apply H3; now left.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
181
apply Rle_trans with ((sqrt (IZR beta)*sqrt (IZR beta)) * (bpow (e-1)*bpow (e-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
182 183 184 185 186 187 188 189
right; ring.
rewrite <- bpow_plus, sqrt_def, <- bpow_1, <- bpow_plus.
right; apply f_equal; ring.
left; apply radix_pos.
apply Rle_ge; auto with real.
(* . *)
simpl in Hx'; contradict H3; assumption.
(* . *)
190
replace (mag beta (x * x)-prec)%Z with (2*e-prec)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
191 192
case (Zle_lt_or_eq _ _ Hb); clear Hb; intros Hb.
(* .. *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
193
apply Rle_lt_trans with (2*(sqrt (IZR beta) * bpow (e - 1))* bpow (e - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
194 195 196 197
2: apply Rmult_lt_compat_r.
2: apply bpow_gt_0.
2: apply Rmult_lt_compat_l; try assumption.
2: auto with real.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
198
apply Rle_trans with ((2 * sqrt (IZR beta)) * bpow (2*e - 1- prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
199 200
2: replace (2*e)%Z with (e+e)%Z by ring; unfold Zminus.
2: repeat rewrite bpow_plus; right; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
201
apply Rle_trans with ((IZR beta + /4)*bpow (2 * e - 1 - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
rewrite Rmult_plus_distr_r.
apply Rplus_le_compat.
rewrite <- bpow_1, <- bpow_plus.
right; apply f_equal; ring.
apply Rmult_le_reg_l with 4%R.
apply Rmult_lt_0_compat; auto with real.
apply Rle_trans with (2*bpow (e - prec + (e - prec))).
right; field.
apply Rle_trans with (bpow (2 * e - 1 - prec)).
2: right; field.
apply Rle_trans with (bpow (1+(e - prec + (e - prec)))).
rewrite (bpow_plus _ 1%Z).
apply Rmult_le_compat_r.
apply bpow_ge_0.
rewrite bpow_1.
217
apply IZR_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
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
generalize (radix_gt_1 beta).
omega.
apply bpow_le.
omega.
apply Rmult_le_compat_r.
apply bpow_ge_0.
assert ((radix_val beta = 2%Z)%Z \/ radix_val beta = 3)%Z.
assert (1 < radix_val beta)%Z by apply radix_gt_1.
omega.
destruct H3; rewrite H3; simpl; interval.
(* .. *)
apply Rlt_le_trans with (2 * (2*bpow (e - 1)+bpow (e-prec)) * bpow (e - prec)).
apply Rlt_le_trans with (4* bpow (e - 1)* bpow (e - prec)
  + bpow (e - prec) * bpow (e - prec)*2).
2: right; ring.
replace (4 * bpow (e - 1) * bpow (e - prec)) with (bpow (2 * e - prec)).
apply Rplus_lt_compat_l.
rewrite <- bpow_plus.
unfold Rdiv; apply Rmult_lt_compat_l.
apply bpow_gt_0.
interval.
replace 4 with (bpow 1).
repeat rewrite <- bpow_plus.
apply f_equal; ring.
rewrite bpow_1, Hb; reflexivity.
apply Rmult_le_compat_r.
apply bpow_ge_0.
apply Rmult_le_compat_l.
auto with real.
apply Rle_trans with (2 * bpow (e - 1) + ulp_flx (2 * bpow (e - 1))).
apply Rplus_le_compat_l.
249 250 251 252
rewrite ulp_neq_0.
2: apply Rmult_integral_contrapositive_currified; apply Rgt_not_eq.
2: apply Rlt_0_2.
2: apply bpow_gt_0.
253
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
254 255
right; apply f_equal.
apply f_equal2; try reflexivity.
256
apply sym_eq, mag_unique.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
257 258 259 260 261 262
rewrite Rabs_right.
split.
rewrite <- (Rmult_1_l (bpow (e-1))) at 1.
apply Rmult_le_compat_r.
apply bpow_ge_0.
auto with real.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
263
apply Rlt_le_trans with (IZR beta*bpow (e - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
264 265 266 267 268 269 270 271
apply Rmult_lt_compat_r.
apply bpow_gt_0.
rewrite Hb; simpl; interval.
rewrite <- bpow_1, <- bpow_plus.
right; apply f_equal; ring.
apply Rle_ge, Rmult_le_pos.
auto with real.
apply bpow_ge_0.
272
rewrite <- succ_eq_pos.
273
apply succ_le_lt...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
274 275 276 277 278 279
apply generic_format_FLX.
exists (Float beta 2 (e-1)).
unfold F2R; now simpl.
apply Zlt_le_trans with (4^1)%Z.
simpl; unfold Z.pow_pos; simpl; omega.
rewrite Hb.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
280 281
apply (Zpower_le (Build_radix 4 eq_refl)).
now apply Zlt_le_weak.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
282
apply Rle_lt_trans with (2:=Hx).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
283
replace (sqrt (IZR beta)) with 2%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
284 285 286
now right.
rewrite Hb; simpl.
apply sym_eq, sqrt_square; auto with real.
287 288 289
apply Rmult_le_pos.
now auto with real.
apply bpow_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
290
apply f_equal with (f:= fun z => (z-prec)%Z).
291
apply sym_eq, mag_unique.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
292 293
rewrite Rabs_right.
split.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
294 295
apply Rle_trans with ((sqrt (IZR beta)*bpow (e-1))*(sqrt (IZR beta)*bpow (e-1))).
apply Rle_trans with ((sqrt (IZR beta)*sqrt (IZR beta)) * (bpow (e-1)*bpow (e-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
2: right; ring.
rewrite <- bpow_plus, sqrt_def, <- bpow_1, <- bpow_plus.
right; apply f_equal; ring.
left; apply radix_pos.
apply Rmult_le_compat.
apply Rmult_le_pos; try apply sqrt_pos; apply bpow_ge_0.
apply Rmult_le_pos; try apply sqrt_pos; apply bpow_ge_0.
now left.
now left.
apply Rlt_le_trans with (bpow (e)*bpow (e)).
2: rewrite <- bpow_plus.
2: right; apply f_equal; ring.
apply Rmult_le_0_lt_compat; try apply H2.
now left.
now left.
apply Rle_ge; auto with real.
Qed.

Guillaume Melquiond's avatar
Guillaume Melquiond committed
314
Lemma ulp_sqr_ulp_lt: forall u, 0 < u ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
315
  (u < sqrt (IZR (radix_val beta)) * bpow (mag beta u-1)) ->
316
  ulp_flx (u * u) + ulp_flx u * ulp_flx u / 2 < 2 * u * ulp_flx u.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
317
Proof with auto with typeclass_instances.
318 319 320 321
intros u Hu.
rewrite 2!ulp_neq_0.
2: now apply Rgt_not_eq.
2: now apply Rgt_not_eq, Rmult_lt_0_compat.
322
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
323
(* *)
324
destruct (mag beta u) as (e,He); simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
325 326 327 328 329 330
intros Cu.
assert (He2:(bpow (e - 1) <= u < bpow e)).
rewrite <- (Rabs_right u).
apply He; auto with real.
apply Rle_ge; now left.
clear He.
331
destruct (mag beta (u*u)) as (i,Hi); simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
332 333 334 335 336 337 338 339 340 341 342 343 344
assert (Hi2:(bpow (i - 1) <= u*u < bpow i)).
rewrite <- (Rabs_right (u*u)).
apply Hi; auto with real.
apply Rle_ge; apply Rmult_le_pos; auto with real.
clear Hi.
assert (i=2*e-1)%Z.
assert (2*e-2 < i /\ i-1 < 2*e-1)%Z;[idtac|omega].
split;apply lt_bpow with beta.
apply Rle_lt_trans with (2:=proj2 Hi2).
replace (2*e-2)%Z with ((e-1)+(e-1))%Z by ring.
rewrite bpow_plus.
apply Rmult_le_compat; try apply bpow_ge_0; apply He2.
apply Rle_lt_trans with (1:=proj1 Hi2).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
345 346
apply Rlt_le_trans with  ((sqrt (IZR beta) * bpow (e - 1))*(sqrt (IZR beta) * bpow (e - 1))).
apply Rlt_trans with ((sqrt (IZR beta) * bpow (e - 1))*u).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
347 348 349 350
now apply Rmult_lt_compat_r.
apply Rmult_lt_compat_l.
apply Rlt_trans with (1:=Hu); assumption.
assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
351
right; apply trans_eq with ((sqrt (IZR beta) * sqrt (IZR beta)) *(bpow (e - 1) * bpow (e - 1)));[ring|idtac].
BOLDO Sylvie's avatar
BOLDO Sylvie committed
352 353
rewrite <- bpow_plus.
rewrite sqrt_sqrt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
354
replace (IZR beta) with (bpow 1).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
355 356 357 358
rewrite <- bpow_plus.
apply f_equal; ring.
apply bpow_1.
generalize (radix_gt_0 beta); intros.
359
left; now apply IZR_lt.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
360 361
rewrite H.
apply Rlt_le_trans with (2 * bpow (e-1) * bpow (e - prec)).
362
change 2%R with (1 + 1)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381
rewrite Rmult_assoc, Rmult_plus_distr_r, Rmult_1_l.
rewrite <- 2!bpow_plus.
replace (e - 1 + (e - prec))%Z with (2 * e - 1 - prec)%Z by ring.
apply Rplus_lt_compat_l.
apply Rmult_lt_reg_l with 2%R; auto with real.
apply Rle_lt_trans with (bpow (e - prec + (e - prec))).
right; field.
apply Rlt_le_trans with (1*bpow (2 * e - 1 - prec)).
rewrite Rmult_1_l; apply bpow_lt.
omega.
apply Rmult_le_compat_r; auto with real.
apply bpow_ge_0.
apply Rmult_le_compat_r.
apply bpow_ge_0.
apply Rmult_le_compat_l; auto with real.
apply He2.
Qed.

Theorem round_flx_sqr_sqrt_exact_case1:
Guillaume Melquiond's avatar
Guillaume Melquiond committed
382
  (x < sqrt (IZR (radix_val beta)) * bpow (mag beta x-1)) ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
383 384 385
    z = x.
Proof with auto with typeclass_instances.
intros Cu.
386
case (Req_dec x (bpow (mag beta x - 1))); intros Hx.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
387 388 389
(* *)
unfold z, y.
rewrite Hx, <- bpow_plus.
390 391
rewrite (round_generic _ _ _ (bpow (mag beta x - 1 + (mag beta x - 1))))...
replace (sqrt (bpow (mag beta x - 1 + (mag beta x - 1)))) with (bpow (mag beta x - 1 )).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
apply round_generic...
apply generic_format_bpow.
unfold FLX_exp; auto with zarith.
apply sym_eq, sqrt_lem_1; try apply bpow_ge_0.
apply sym_eq, bpow_plus.
apply generic_format_bpow.
unfold FLX_exp; auto with zarith.
(* *)
assert (0 < x*x) by now apply Rmult_lt_0_compat.
assert (0 <= 1 - ulp_flx x / 2 / x).
apply Rplus_le_reg_l with (ulp_flx x / 2 / x); ring_simplify.
apply Rmult_le_reg_l with 2; auto with real.
apply Rmult_le_reg_l with x; auto with real.
apply Rle_trans with (ulp_flx x).
right; field; auto with real.
apply Rle_trans with x.
apply ulp_le_id; auto.
409
lra.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
410 411 412 413 414 415 416
assert (0 <= 1 - ulp_flx (x * x) / 2 / (x * x)).
apply Rplus_le_reg_l with (ulp_flx (x*x) / 2 / (x*x)); ring_simplify.
apply Rmult_le_reg_l with 2; auto with real.
apply Rmult_le_reg_l with (x*x); auto with real.
apply Rle_trans with (ulp_flx (x*x)).
right; field; auto with real.
apply Rle_trans with (x*x).
417
rewrite ulp_neq_0; try now apply Rgt_not_eq.
418
unfold cexp, FLX_exp.
419
destruct (mag beta (x*x)) as (e,He); simpl.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
420 421 422 423 424
apply Rle_trans with (bpow (e-1)).
apply bpow_le; auto with zarith.
rewrite <- (Rabs_right (x*x)).
apply He; auto with real.
apply Rle_ge; auto with real.
425
lra.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
426 427 428 429
assert (0 <= 1 + ulp_flx x / 2 / x).
rewrite <- (Rplus_0_l 0).
apply Rplus_le_compat; auto with real.
unfold Rdiv; apply Rmult_le_pos.
430
apply Rmult_le_pos; auto with real; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
431 432 433 434 435
left; auto with real.
assert (0 <= 1 + ulp_flx (x * x) / 2 / (x * x)).
rewrite <- (Rplus_0_l 0).
apply Rplus_le_compat; auto with real.
unfold Rdiv; apply Rmult_le_pos.
436
apply Rmult_le_pos; auto with real; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
437 438 439 440
left; auto with real.
(* *)
apply sym_eq, Rle_antisym.
(* . *)
441
apply round_N_ge_midp; try assumption...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
442
apply Rle_lt_trans with (x*(1-ulp_flx x / 2 / x)).
443
rewrite pred_eq_pos;[idtac|now left].
444
unfold pred_pos; rewrite Req_bool_false; trivial.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469
right; field; auto with real.
apply Rlt_le_trans with (sqrt ((x*x)*(1-ulp_flx (x*x)/2/(x*x)))).
rewrite sqrt_mult; auto with real.
rewrite sqrt_square; auto with real.
apply Rmult_lt_compat_l.
assumption.
rewrite <- (sqrt_square (1 - ulp_flx x / 2 / x)); try assumption.
apply sqrt_lt_1_alt.
split.
apply Rmult_le_pos; assumption.
apply Rmult_lt_reg_l with (2*x*x).
rewrite Rmult_assoc; apply Rmult_lt_0_compat; auto with real.
apply Rplus_lt_reg_r with (-2*x*x+2*x*ulp_flx x+ ulp_flx(x*x)).
apply Rle_lt_trans with (ulp_flx (x*x) + ulp_flx (x)*ulp_flx (x)/2).
right; field.
auto with real.
apply Rlt_le_trans with (2*x*ulp_flx x).
2: right; field; auto with real.
apply ulp_sqr_ulp_lt; auto.
apply sqrt_le_1_alt.
apply Rplus_le_reg_l with (-y+ulp_flx (x * x)/2).
apply Rle_trans with (-(y-x*x));[right; field; auto with real|idtac].
apply Rle_trans with (/2*ulp_flx (x*x));[idtac|right; field; auto with real].
apply Rle_trans with (1:=RRle_abs _).
rewrite Rabs_Ropp.
470
apply error_le_half_ulp...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
471
(* . *)
472
apply round_N_le_midp; try assumption...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
473
apply Rlt_le_trans with (x*(1+ulp_flx x / 2 / x)).
474
2: rewrite succ_eq_pos; try now left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
475 476 477 478 479 480 481
2: right; field; auto with real.
apply Rle_lt_trans with (sqrt ((x*x)*(1+ulp_flx (x*x)/2/(x*x)))).
apply sqrt_le_1_alt.
apply Rplus_le_reg_l with (-x*x).
apply Rle_trans with (y-x*x);[right; ring|idtac].
apply Rle_trans with (/2*ulp_flx (x*x));[idtac|right; field; auto with real].
apply Rle_trans with (1:=RRle_abs _).
482
apply error_le_half_ulp...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503
rewrite sqrt_mult; auto with real.
rewrite sqrt_square; auto with real.
apply Rmult_lt_compat_l.
assumption.
rewrite <- (sqrt_square (1 + ulp_flx x / 2 / x)); try assumption.
apply sqrt_lt_1_alt.
split; try assumption.
apply Rmult_lt_reg_l with (2*x*x).
rewrite Rmult_assoc; apply Rmult_lt_0_compat; auto with real.
apply Rplus_lt_reg_r with (-2*x*x - ulp_flx x*ulp_flx x/2).
apply Rle_lt_trans with (ulp_flx (x*x) - ulp_flx x*ulp_flx x/2).
right; field.
auto with real.
apply Rlt_le_trans with (2*x*ulp_flx x).
2: right; field; auto with real.
apply Rle_lt_trans with (ulp_flx (x * x) + ulp_flx x * ulp_flx x / 2).
2: apply ulp_sqr_ulp_lt; auto.
apply Rplus_le_compat_l.
apply Rplus_le_reg_r with (ulp_flx x * ulp_flx x / 2).
apply Rle_trans with 0;[right; ring|idtac].
apply Rle_trans with (ulp_flx x * ulp_flx x).
504
apply Rmult_le_pos; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
505 506 507 508 509
right; field.
Qed.

Theorem round_flx_sqr_sqrt_aux2: forall n,
 (0 <= n)%Z ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
510 511
 (2*IZR n * bpow (prec-1) * ulp_flx x * (1+bpow (1-prec)/2) < x) ->
  round_flx3(x/(x-IZR n*ulp_flx x)) <= 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
512 513
Proof with auto with typeclass_instances.
intros n Hnp Hn.
514
apply round_N_le_midp...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
515 516 517
replace 1 with (bpow 0) by reflexivity.
apply generic_format_bpow.
unfold FLX_exp; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
518 519
assert (0 < x - IZR n*ulp_flx x).
apply Rplus_lt_reg_r with (IZR n*ulp_flx x); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
520
apply Rle_lt_trans with (2:=Hn).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
521 522
apply Rle_trans with (IZR n*(ulp_flx x*((1*1)*(1+0))));[right; ring|idtac].
apply Rle_trans with (IZR n*(ulp_flx x*(2*bpow (prec - 1)* (1 + bpow (1 - prec) / 2))));[idtac|right; ring].
BOLDO Sylvie's avatar
BOLDO Sylvie committed
523
apply Rmult_le_compat_l.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
524
now apply IZR_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
525
apply Rmult_le_compat_l.
526
unfold ulp; apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
527 528 529 530 531 532 533 534 535 536 537 538
apply Rmult_le_compat.
rewrite Rmult_1_l; apply Rle_0_1.
rewrite Rplus_0_r; apply Rle_0_1.
apply Rmult_le_compat; try apply Rle_0_1.
auto with real.
replace 1 with (bpow 0) by reflexivity.
apply bpow_le.
omega.
apply Rplus_le_compat_l.
unfold Rdiv; apply Rmult_le_pos.
apply bpow_ge_0.
auto with real.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
539
apply Rmult_lt_reg_r with (x - IZR n*ulp_flx x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
540 541 542 543 544
assumption.
unfold Rdiv; rewrite Rmult_assoc.
rewrite Rinv_l.
rewrite Rmult_1_r.
2: auto with real.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
545
apply Rplus_lt_reg_r with (-x+IZR n*ulp_flx x + IZR n*ulp_flx x* ulp_flx 1 * / 2); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
546 547 548 549
apply Rmult_lt_reg_l with (2*/ulp_flx 1).
apply Rmult_lt_0_compat.
auto with real.
apply Rinv_0_lt_compat.
550 551
rewrite ulp_neq_0; try apply bpow_gt_0.
apply Rgt_not_eq, Rlt_0_1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
552 553 554 555 556 557 558 559 560
apply Rlt_le_trans with x.
apply Rle_lt_trans with (2:=Hn).
replace  (ulp_flx 1) with (bpow (1-prec)).
rewrite <- bpow_opp.
replace ((-(1-prec)))%Z with (prec -1)%Z by ring.
right; unfold Rdiv; ring.
replace 1 with (bpow 0) by reflexivity.
rewrite ulp_bpow.
apply f_equal; unfold FLX_exp; omega.
561
rewrite succ_eq_pos.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
562
right; field.
563 564 565 566
apply Rgt_not_eq.
rewrite ulp_neq_0; try apply bpow_gt_0.
apply Rgt_not_eq, Rlt_0_1.
apply Rle_0_1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
567 568 569
Qed.

End Sec1.
570

BOLDO Sylvie's avatar
BOLDO Sylvie committed
571
Section Sec2.
572

BOLDO Sylvie's avatar
BOLDO Sylvie committed
573 574 575 576 577 578 579
Open Scope R_scope.

Variable beta : radix.
Notation bpow e := (bpow beta e).

Variable prec : Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
580 581
Variable choice1 : Z -> bool.
Variable choice2 : Z -> bool.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
582 583

Notation format := (generic_format beta (FLX_exp prec)).
584 585
Notation round_flx1 :=(round beta (FLX_exp prec) (Znearest choice1)).
Notation round_flx2 :=(round beta (FLX_exp prec) (Znearest choice2)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
586 587 588 589 590 591 592 593
Notation ulp_flx :=(ulp beta (FLX_exp prec)).
Notation pred_flx := (pred beta (FLX_exp prec)).

Hypothesis pGt1: (1 < prec)%Z.

Variables x:R.
Hypothesis xPos: 0 < x.
Hypothesis Fx: format x.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
594
Hypothesis Hx: (sqrt (IZR (radix_val beta)) * bpow (mag beta x-1) < x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
595 596 597 598 599 600

Variable k:Z.
Hypothesis kpos: (0 <= k)%Z.
Hypothesis kLe: (k < radix_val beta)%Z.
Hypothesis kradix2 : (k = 0)%Z  \/  (2 < radix_val beta)%Z.

601 602
Let y:=round_flx1(x*x).
Let z:=round_flx2(sqrt y).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
603
Let xk :=  x-IZR k*ulp_flx x.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
604

605
Lemma xkgt:  bpow (mag beta x - 1) < xk.
606
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
607 608 609 610
unfold xk.
case kradix2.
intros V; rewrite V; simpl; ring_simplify.
apply Rle_lt_trans with (2:=Hx).
611
rewrite <- (Rmult_1_l (bpow (mag beta x - 1))) at 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
612 613 614 615 616
apply Rmult_le_compat_r.
apply bpow_ge_0.
assert (2 <= beta)%Z.
clear; destruct beta as (v, Hr); simpl.
now apply Zle_bool_imp_le.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
617
apply IZR_le in H.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
618 619
simpl in H; interval.
intros Hb.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
620 621
apply Rle_lt_trans with (sqrt (IZR beta) * bpow (mag beta x - 1)
    - IZR k * ulp_flx x).
622
rewrite ulp_neq_0; try now apply Rgt_not_eq.
623
unfold cexp, FLX_exp.
624
apply Rle_trans with (bpow (mag beta x - 1)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
625
   *(sqrt (IZR beta) -IZR k * bpow (1-prec))).
626
rewrite <- (Rmult_1_r (bpow (mag beta x - 1))) at 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
627 628
apply Rmult_le_compat_l.
apply bpow_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
629 630
apply Rplus_le_reg_l with (IZR k * bpow (1 - prec)).
apply Rle_trans with ((1-/IZR beta) * bpow (2 - prec)+1).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
631
apply Rplus_le_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
632
apply Rle_trans with (((1-/IZR beta)*IZR beta) * bpow (1 - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
633 634
apply Rmult_le_compat_r.
apply bpow_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
635 636
apply Rle_trans with (IZR (beta-1)).
apply IZR_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
637
omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
638
rewrite minus_IZR; right.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
639 640
simpl; field.
apply Rgt_not_eq, Rlt_gt.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
641
apply IZR_lt, radix_gt_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
642 643 644
right; rewrite Rmult_assoc; apply f_equal.
rewrite <- bpow_1, <- bpow_plus.
apply f_equal; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
645
ring_simplify (IZR k * bpow (1 - prec) + (sqrt (IZR beta) - IZR k * bpow (1 - prec))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
646
assert (H':(3 <= beta)%Z) by omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
647 648
apply IZR_le in H'; simpl in H'.
assert (H'':(1 <= IZR beta)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
649 650 651 652 653 654 655
apply Rle_trans with (2:=H').
apply Rplus_le_reg_l with (-1); ring_simplify.
auto with real.
(* because p=2 is possible *)
case (Zle_lt_or_eq 3 beta).
omega.
intros H1; assert (H1':(4 <= beta)%Z) by omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
656
apply IZR_le in H1'; simpl in H1'.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
657 658 659
apply Rle_trans with (1*1 +1).
apply Rplus_le_compat_r.
apply Rmult_le_compat.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
660
apply Rplus_le_reg_l with (/IZR beta); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
661 662 663 664 665 666
rewrite <- Rinv_1.
apply Rle_Rinv.
apply Rlt_0_1.
apply Rlt_le_trans with (2:=H''); auto with real.
exact H''.
apply bpow_ge_0.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
667
apply Rplus_le_reg_l with (/IZR beta-1); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
668 669 670 671 672 673 674 675
left; apply Rinv_0_lt_compat.
apply Rlt_le_trans with (2:=H''); auto with real.
apply Rle_trans with (bpow 0).
apply bpow_le.
omega.
right; reflexivity.
interval.
intros Hb'.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
676
apply Rle_trans with ((1 - / IZR beta) *1 +1).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
677 678
apply Rplus_le_compat_r.
apply Rmult_le_compat_l.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
679
apply Rplus_le_reg_l with (/IZR beta); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
680 681 682 683 684 685 686 687 688 689 690 691
rewrite <- Rinv_1.
apply Rle_Rinv.
apply Rlt_0_1.
apply Rlt_le_trans with (2:=H''); auto with real.
exact H''.
apply Rle_trans with (bpow 0).
apply bpow_le.
omega.
right; reflexivity.
rewrite <- Hb'; simpl; interval.
right; ring_simplify.
apply f_equal.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
692
apply trans_eq with (IZR k *(bpow (mag beta x - 1) * bpow (1 - prec))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
693 694 695 696 697 698 699 700 701
ring.
apply f_equal.
rewrite <- bpow_plus.
apply f_equal.
ring.
apply Rplus_lt_compat_r.
assumption.
Qed.

702
Lemma xklt: xk < bpow (mag beta x).
703
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
704 705 706
apply Rle_lt_trans with x.
apply Rplus_le_reg_l with (-xk); unfold xk; ring_simplify.
apply Rmult_le_pos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
707
apply IZR_le, kpos.
708
apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
709
apply Rle_lt_trans with (1:=RRle_abs _).
710
apply bpow_mag_gt.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
711 712 713
Qed.

Lemma xkpos: 0  < xk.
714
Proof.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
715 716 717 718 719 720 721 722 723
apply Rle_lt_trans with (2:=xkgt).
apply bpow_ge_0.
Qed.

Lemma format_xminusk: format xk.
Proof with auto with typeclass_instances.
apply generic_format_FLX.
unfold generic_format in Fx.
exists  (Float beta (Ztrunc (scaled_mantissa beta (FLX_exp prec) x) - k)
724
                    (cexp beta (FLX_exp prec) x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
725
unfold xk; rewrite Fx at 1; unfold F2R; simpl.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
726
rewrite minus_IZR; ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
727
apply f_equal.
728
rewrite ulp_neq_0; try now apply Rgt_not_eq.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
729 730 731 732 733 734 735
apply Rmult_comm; apply f_equal.
simpl.
rewrite Z.abs_eq.
apply Zle_lt_trans with (Ztrunc (scaled_mantissa beta (FLX_exp prec) x) - 0)%Z.
apply Zplus_le_compat_l.
omega.
rewrite Zminus_0_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
736
apply lt_IZR.
737
apply Rmult_lt_reg_r with (bpow (cexp beta (FLX_exp prec) x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
738 739 740 741
apply bpow_gt_0.
apply Rle_lt_trans with x.
rewrite Fx at 3.
right; unfold F2R; reflexivity.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
742
rewrite IZR_Zpower.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
743
rewrite <- bpow_plus.
744
unfold cexp, FLX_exp.
745
replace  (prec + (mag beta x - prec))%Z with (mag beta x+0)%Z by ring.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
746 747 748
rewrite Zplus_0_r.
apply Rle_lt_trans with (Rabs x).
apply RRle_abs.
749
apply bpow_mag_gt...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
750
omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
751
apply le_IZR.
752
apply Rmult_le_reg_r with (bpow (cexp beta (FLX_exp prec) x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
753 754 755 756
apply bpow_gt_0.
rewrite Rmult_0_l.
apply Rle_trans with xk.
left; apply xkpos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
757
right; unfold xk; rewrite Fx at 1; unfold F2R; simpl; rewrite minus_IZR; ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
758
apply f_equal.
759
rewrite ulp_neq_0; try now apply Rgt_not_eq.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
760 761 762
apply Rmult_comm; apply f_equal.
Qed.

Guillaume Melquiond's avatar
Guillaume Melquiond committed
763
Theorem round_flx_sqr_sqrt_aux1:
764
  (/ 2 * bpow (mag beta x) <
Guillaume Melquiond's avatar
Guillaume Melquiond committed
765 766
      (2 * IZR k + 1) * x -
           (IZR k + / 2) * (IZR k + / 2) * bpow (mag beta x - prec)) ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
767 768 769
  xk <= z.
Proof with auto with typeclass_instances.
intros V.
770
apply round_N_ge_midp...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
771
apply format_xminusk.
772 773
assert (Z:(mag_val beta xk (mag beta xk) = mag beta x)%Z).
apply mag_unique.
774
rewrite Rabs_right; try split.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
775 776 777
left; now apply xkgt.
now apply xklt.
apply Rle_ge; left; now apply xkpos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
778
apply Rle_lt_trans with (x - IZR k * ulp_flx x - ulp_flx x / 2).
779
rewrite pred_eq_pos.
780 781 782 783 784 785 786 787
2: left; apply xkpos.
unfold pred_pos; rewrite Req_bool_false.
2: apply Rgt_not_eq, Rlt_gt.
2: apply Rle_lt_trans with (2:=xkgt).
2: right; apply f_equal, f_equal2...
replace (ulp_flx xk) with (ulp_flx x).
unfold xk; right; field.
rewrite 2!ulp_neq_0; try apply Rgt_not_eq; try assumption.
788
unfold cexp; now rewrite Z.
789
apply xkpos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
790
apply Rle_lt_trans with (x-(IZR k+/2)*ulp_flx x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
791
right; unfold Rdiv; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
792
assert (0 <= x - (IZR k + / 2) * ulp_flx x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
793 794 795
apply Rplus_le_reg_l with (/2*ulp_flx x).
rewrite Rplus_0_r.
apply Rle_trans with xk.
796
apply Rle_trans with (1*bpow (mag beta x - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
797 798
apply Rmult_le_compat.
auto with real.
799
apply ulp_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
800
interval.
801
rewrite ulp_neq_0; try now apply Rgt_not_eq.
802
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
803 804 805 806
apply bpow_le; omega.
rewrite Rmult_1_l.
left; apply xkgt.
unfold xk; right; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
807
rewrite <- (sqrt_square (x - (IZR k + / 2) * ulp_flx x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
808 809 810 811
2: assumption.
apply sqrt_lt_1_alt.
split.
apply Rmult_le_pos; assumption.
812
apply Rlt_le_trans with (x*x - /2*bpow (2 * mag beta x  - prec)).
813
rewrite ulp_neq_0; try now apply Rgt_not_eq.
814
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
815 816
apply Rplus_lt_reg_r with (-x*x).
apply Rle_lt_trans with
Guillaume Melquiond's avatar
Guillaume Melquiond committed
817 818
  (- (bpow (mag beta x - prec)*((2*IZR k +1)*x -
            (IZR k + / 2)*  (IZR k + / 2) * bpow (mag beta x - prec)))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
819
right; field.
820 821
apply Rlt_le_trans with (- (bpow (mag beta x - prec)*
    (/2*bpow (mag beta x)))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
822 823 824 825
apply Ropp_lt_contravar.
apply Rmult_lt_compat_l.
apply bpow_gt_0.
exact V.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
826
right; apply trans_eq with
827
  ((-/2)*(bpow (mag beta x - prec)*bpow (mag beta x))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
828 829
ring.
rewrite <- bpow_plus.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
830
apply trans_eq with
831
  ((-/2)*bpow (2 * mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
832 833 834
apply f_equal.
apply f_equal; ring.
ring.
835
apply Rplus_le_reg_l with (-y+/2*bpow (2 * mag beta x  - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
836 837 838 839 840
ring_simplify.
apply Rle_trans with (-(y-x*x)).
right; ring.
apply Rle_trans with (1:=RRle_abs _).
rewrite Rabs_Ropp.
841
apply Rle_trans with (1:=error_le_half_ulp _ _ _ _)...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
842 843
apply Rmult_le_compat_l.
left; auto with real.
844 845
rewrite ulp_neq_0.
2: apply Rmult_integral_contrapositive_currified; now apply Rgt_not_eq.
846
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
847 848
apply bpow_le.
apply Zplus_le_compat_r.
849
apply mag_le_bpow.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
850 851
auto with real.
rewrite Rabs_mult.
852
replace (2*mag beta x)%Z with (mag beta x+mag beta x)%Z by ring.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
853
rewrite bpow_plus.
854
apply Rmult_le_0_lt_compat; try apply Rabs_pos; apply bpow_mag_gt.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
855 856
Qed.

857

BOLDO Sylvie's avatar
BOLDO Sylvie committed
858
Theorem round_flx_sqr_sqrt_aux1_simpl:
Guillaume Melquiond's avatar
Guillaume Melquiond committed
859
  (/ 2 * bpow (mag beta x) + bpow (2+mag beta x - prec) <= (2 * IZR k + 1) * x)
860
  -> xk <= z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
861 862
Proof.
intros H; apply round_flx_sqr_sqrt_aux1.
863
apply Rplus_lt_reg_r with (bpow (2 + mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
864
apply Rle_lt_trans with (1:=H).
Guillaume Melquiond's avatar
Guillaume Melquiond committed
865 866
apply Rplus_lt_reg_r with (-(2 * IZR k + 1) * x + (IZR k + / 2) * (IZR k + / 2) * bpow (mag beta x - prec)).
apply Rle_lt_trans with (((IZR k + / 2) * (IZR k + / 2) )* bpow (mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
867
right; ring.
868
apply Rlt_le_trans with (bpow (2 + mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
869 870 871 872 873 874
2: right; ring.
unfold Zminus; rewrite <- Zplus_assoc.
rewrite bpow_plus with (e1:=2%Z).
apply Rmult_lt_compat_r.
apply bpow_gt_0.
simpl; unfold Z.pow_pos; simpl.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
875 876 877
rewrite Zmult_1_r, mult_IZR.
assert (IZR k + /2 < IZR beta).
apply Rle_lt_trans with (IZR (beta -1) + /2).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
878
apply Rplus_le_compat_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
879
apply IZR_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
880
omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
881 882
rewrite minus_IZR; simpl.
apply Rplus_lt_reg_r with (-IZR beta + /2).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
883 884 885 886
field_simplify.
unfold Rdiv; apply Rmult_lt_compat_r.
apply Rinv_0_lt_compat, Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
exact Rlt_0_1.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
887
assert (0 <= IZR k + / 2).
888
replace 0 with (0+0) by (simpl; ring); apply Rplus_le_compat.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
889
apply IZR_le, kpos.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
890 891 892 893 894 895 896
apply Rlt_le, Rinv_0_lt_compat, Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
now apply Rmult_le_0_lt_compat.
Qed.




BOLDO Sylvie's avatar
BOLDO Sylvie committed
897 898 899
End Sec2.

Section Sec3.
900

BOLDO Sylvie's avatar
BOLDO Sylvie committed
901 902 903 904 905 906 907
Open Scope R_scope.

Variable beta : radix.
Notation bpow e := (bpow beta e).

Variable prec : Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
908 909 910
Variable choice1 : Z -> bool.
Variable choice2 : Z -> bool.
Variable choice3 : Z -> bool.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
911 912

Notation format := (generic_format beta (FLX_exp prec)).
913 914 915
Notation round_flx1 :=(round beta (FLX_exp prec) (Znearest choice1)).
Notation round_flx2 :=(round beta (FLX_exp prec) (Znearest choice2)).
Notation round_flx3 :=(round beta (FLX_exp prec) (Znearest choice3)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
916 917 918
Notation ulp_flx :=(ulp beta (FLX_exp prec)).
Notation pred_flx := (pred beta (FLX_exp prec)).

BOLDO Sylvie's avatar
BOLDO Sylvie committed
919
Hypothesis pGt2: (2 < prec)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
920 921 922 923 924 925

Hypothesis pGt1: (1 < prec)%Z.

Variables x:R.
Hypothesis xPos: 0 < x.
Hypothesis Fx: format x.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
926
Hypothesis Hx: (sqrt (IZR (radix_val beta)) * bpow (mag beta x-1) < x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
927

928 929
Let y:=round_flx1(x*x).
Let z:=round_flx2(sqrt y).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
930 931 932 933 934 935 936 937

Theorem round_flx_sqr_sqrt_exact_aux: (radix_val beta <= 4)%Z ->
    z=x.
Proof with auto with typeclass_instances.
intros Hb.
apply Rle_antisym.
(* . *)
apply round_flx_sqr_sqrt_le...
Guillaume Melquiond's avatar
Guillaume Melquiond committed
938
(* . *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
939 940 941 942
assert (((radix_val beta = 2)%Z \/ (radix_val beta=3)%Z) \/ (radix_val beta=4)%Z).
generalize (radix_gt_1 beta); omega.
destruct H.
(* .. radix is 2 or 3 *)
943
apply Rle_trans with (x - 0 * ulp_flx x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
944 945 946 947
right; simpl; ring.
apply round_flx_sqr_sqrt_aux1...
omega.
apply radix_gt_0.
948
apply Rlt_le_trans with (x-/4*bpow (mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
949
2: right; simpl; field.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
950
apply Rle_lt_trans with (sqrt (IZR beta) * bpow (mag beta x - 1)
951
   - / 4 * bpow (mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
952
2: apply Rplus_lt_compat_r; assumption.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
953
apply Rle_trans with ((IZR beta / 2)*bpow (mag beta x-1)).
954
replace (bpow (mag beta x)) with (bpow (1+(mag beta x-1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
955 956 957
rewrite bpow_plus, bpow_1.
right; unfold Rdiv; ring.
apply f_equal; ring.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
958
apply Rle_trans with ((sqrt (IZR beta) - /4* / IZR beta)
959
    * bpow (mag beta x-1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
960 961 962 963 964 965 966 967
apply Rmult_le_compat_r.
apply bpow_ge_0.
destruct H; rewrite H; simpl; interval.
rewrite Rmult_minus_distr_r.
apply Rplus_le_compat_l.
apply Ropp_le_contravar.
rewrite Rmult_assoc; apply Rmult_le_compat_l.
left; apply Rinv_0_lt_compat, Rmult_lt_0_compat; auto with real.
968
apply Rle_trans with (bpow (-1+(mag beta x - 1))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
969 970 971 972 973 974
apply bpow_le; omega.
rewrite bpow_plus.
right; apply f_equal2; try reflexivity.
rewrite <- bpow_1, <- bpow_opp.
apply f_equal; reflexivity.
(* .. radix is 4 *)
975
assert (2 * bpow (mag beta x - 1) < x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
976 977 978 979 980
apply Rle_lt_trans with (2:=Hx).
right; apply f_equal2; try reflexivity.
rewrite H; simpl.
apply sym_eq, sqrt_square; auto with real.
(* ... *)
981 982
assert ((2 * bpow (mag beta x - 1)+1*bpow (mag beta x - prec)) <= x).
assert (0 < 2 * bpow (mag beta x - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
983 984 985
apply Rmult_lt_0_compat.
auto with real.
apply bpow_gt_0.
986
assert (bpow (mag beta x - prec)=ulp_flx (2 * bpow (mag beta x - 1))).
987
rewrite ulp_neq_0; try now apply Rgt_not_eq.
988
unfold cexp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
989 990
apply f_equal.
apply f_equal2; try reflexivity.
991
apply sym_eq, mag_unique.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
992 993 994
rewrite Rabs_right.
2: apply Rle_ge; left; assumption.
split.
995
apply Rle_trans with (1*bpow (mag beta x - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
996 997 998 999
right; ring.
apply Rmult_le_compat_r.
apply bpow_ge_0.
auto with real.
1000
apply Rlt_le_trans with (4*bpow (mag beta x - 1)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1001 1002 1003 1004 1005 1006
apply Rmult_lt_compat_r.
apply bpow_gt_0.
interval.
rewrite <- H, <- bpow_1, <- bpow_plus.
right; apply f_equal; ring.
rewrite H2, Rmult_1_l.
1007
rewrite <- succ_eq_pos;[idtac|now left].
1008
apply succ_le_lt...
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1009
apply generic_format_FLX.
1010
exists (Float beta 2 (mag beta x -1)).
1011
easy.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1012 1013
rewrite H; apply Zlt_le_trans with (4^1)%Z.
simpl; unfold Z.pow_pos; simpl; omega.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1014 1015
apply (Zpower_le (Build_radix 4 eq_refl)).
now apply Zlt_le_weak.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1016
(* ... *)
1017
apply Rle_trans with (x - 0 * ulp_flx x).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1018 1019 1020 1021
right; simpl; ring.
apply round_flx_sqr_sqrt_aux1...
omega.
apply radix_gt_0.
1022
apply Rlt_le_trans with (x-/4*bpow (mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1023
2: right; simpl; field.
1024 1025
apply Rlt_le_trans with ( (2* bpow (mag beta x - 1) +  1 * bpow (mag beta x - prec))
   - / 4 * bpow (mag beta x - prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1026 1027
2: apply Rplus_le_compat_r; assumption.
apply Rlt_le_trans with ((/2 + (1-/4)*bpow (-prec))
1028
    * bpow (mag beta x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042
apply Rmult_lt_compat_r.
apply bpow_gt_0.
apply Rle_lt_trans with (/2+0);[right; ring|idtac].
apply Rplus_lt_compat_l.
apply Rmult_lt_0_compat.
interval.
apply bpow_gt_0.
unfold Zminus; repeat rewrite bpow_plus.
replace (bpow (- (1))) with (/4).
right; field.
rewrite bpow_opp, bpow_1, H; reflexivity.
Qed.


1043
Let k:=(Zceil (x*bpow (1-mag beta x)/(2+bpow (1-prec))) -1)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1044

BOLDO Sylvie's avatar
BOLDO Sylvie committed
1045
Lemma kpos: (0 <= k)%Z.
1046
Proof.
1047
assert (0 < x*bpow (1-mag beta x)/(2+bpow (1-prec))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1048
apply Fourier_util.Rlt_mult_inv_pos.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1049 1050
apply Rmult_lt_0_compat.
assumption.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1051
apply bpow_gt_0.
1052 1053 1054
apply Rplus_lt_0_compat.
now apply IZR_lt.
apply bpow_gt_0.
1055
assert (0 < Zceil (x * bpow (1 - mag beta x) / (2+bpow (1-prec))))%Z.
1056
apply lt_IZR.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1057 1058
apply Rlt_le_trans with (1:=H).
apply Zceil_ub.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1059
unfold k; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1060 1061 1062
Qed.

Lemma kLe: (k < radix_val beta)%Z.
1063
Proof.
1064
cut (Zceil (x * bpow (1 - mag beta x) / (2+bpow (1-prec))) <= beta)%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1065 1066 1067 1068 1069 1070
unfold k; omega.
apply Zceil_glb.
apply Rle_trans with (bpow 1 / 1).
unfold Rdiv; apply Rmult_le_compat.
apply Rmult_le_pos; try apply bpow_ge_0; now left.
apply Rlt_le, Rinv_0_lt_compat.
1071 1072 1073
apply Rplus_lt_0_compat.
now apply IZR_lt.
apply bpow_gt_0.
1074
apply Rle_trans with (bpow (mag beta x)*bpow (1 - mag beta x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1075 1076 1077
apply Rmult_le_compat_r.
apply bpow_ge_0.
apply Rle_trans with (1:=RRle_abs _).
1078
left; apply bpow_mag_gt.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1079
rewrite <- bpow_plus.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1080 1081 1082 1083
apply bpow_le; omega.
apply Rinv_le.
exact Rlt_0_1.
apply Rplus_le_reg_l with (-1); ring_simplify.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1084
apply Rlt_le, Rle_lt_0_plus_1, bpow_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1085
rewrite bpow_1; right; field.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1086 1087
Qed.

Guillaume Melquiond's avatar
Guillaume Melquiond committed
1088
Lemma kLe2: (k <= Zceil (IZR(radix_val beta)/ 2) -1)%Z.
1089
Proof.
1090
cut (Zceil (x * bpow (1 - mag beta x) / (2+bpow (1-prec)))
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1091
   <= Zceil (IZR(radix_val beta)/ 2))%Z.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1092 1093 1094
unfold k; omega.
apply Zceil_glb.
apply Rle_trans with (bpow 1 / 2).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1095 1096 1097 1098
unfold Rdiv; apply Rmult_le_compat.
apply Rmult_le_pos.
now left.
apply bpow_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1099
apply Rlt_le, Rinv_0_lt_compat.
1100 1101 1102
apply Rplus_lt_0_compat.
now apply IZR_lt.
apply bpow_gt_0.
1103
apply Rle_trans with (bpow (mag beta x)*bpow (1 - mag beta x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1104 1105 1106
apply Rmult_le_compat_r.
apply bpow_ge_0.
apply Rle_trans with (1:=RRle_abs _).
1107
left; apply bpow_mag_gt.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1108 1109
rewrite <- bpow_plus.
apply bpow_le; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1110 1111 1112 1113
apply Rinv_le.
apply Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
apply Rplus_le_reg_l with (-2); ring_simplify.
apply bpow_ge_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1114 1115 1116
rewrite bpow_1.
apply Zceil_ub.
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1117

BOLDO Sylvie's avatar
BOLDO Sylvie committed
1118 1119


1120
Lemma round_flx_sqr_sqrt_snd_deg:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1121
  (((radix_val beta=5)%Z /\ (3 < prec)%Z) \/ (5 < radix_val beta)%Z) ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1122 1123
    sqrt (IZR beta) + / 4 * bpow (3 - prec) <=
      IZR beta * (2 - bpow (1 - prec)) / (2 * (2 + bpow (1 - prec))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1124 1125 1126 1127
Proof.
intros H; destruct H.
(* radix=5 *)
destruct H as (H1,H2).
1128
apply Rle_trans with (sqrt (IZR beta) + / (4 *IZR beta)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1129 1130 1131 1132 1133 1134 1135 1136 1137 1138
apply Rplus_le_compat_l.
rewrite (Rinv_mult_distr 4).
apply Rmult_le_compat_l.
apply Rlt_le, Rinv_0_lt_compat, Rmult_lt_0_compat; apply Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
apply Rle_trans with (bpow (-(1))).
apply bpow_le; omega.
right; rewrite bpow_opp.
apply f_equal, bpow_1.
apply Rgt_not_eq, Rmult_lt_0_compat; apply Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
apply Rgt_not_eq, radix_pos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1139
apply Rle_trans with (IZR beta * (2-/(IZR beta*IZR beta)) / (2* (2 + /(IZR beta*IZR beta)))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152
rewrite H1;simpl; interval.
unfold Rdiv; rewrite 2!Rmult_assoc.
apply Rmult_le_compat_l.
left; apply radix_pos.
apply Rmult_le_compat.
rewrite H1; simpl; interval.
rewrite H1; simpl; interval.
apply Rplus_le_reg_l with (-2); ring_simplify.
apply Ropp_le_contravar.
apply Rle_trans with (bpow (-(2))).
apply bpow_le; omega.
right; rewrite bpow_opp.
apply f_equal; simpl; unfold Z.pow_pos; simpl.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1153
rewrite <- mult_IZR; apply f_equal; ring.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1154 1155 1156
apply Rinv_le.
apply Rmult_lt_0_compat.
apply Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
1157 1158 1159
apply Rplus_lt_0_compat.
now apply IZR_lt.
apply bpow_gt_0.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1160 1161 1162 1163 1164 1165 1166
apply Rmult_le_compat_l.
apply Rlt_le,Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
apply Rplus_le_compat_l.
apply Rle_trans with (bpow (-(2))).
apply bpow_le; omega.
right; rewrite bpow_opp.
apply f_equal; simpl; unfold Z.pow_pos; simpl.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1167
rewrite <- mult_IZR; apply f_equal; ring.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1168
(* radix > 5 *)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1169
apply Rle_trans with (sqrt (IZR beta) + /4*1).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1170 1171 1172 1173 1174 1175 1176
apply Rplus_le_compat_l.
apply Rmult_le_compat_l.
apply Rlt_le, Rinv_0_lt_compat,  Rmult_lt_0_compat; apply Rle_lt_0_plus_1, Rlt_le, Rlt_0_1.
apply Rle_trans with (bpow 0).
apply bpow_le; omega.
right; reflexivity.
rewrite Rmult_1_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1177 1178 1179 1180 1181 1182 1183
assert (6 <= IZR beta).
apply IZR_le; omega.
apply Rle_trans with (IZR beta*(12/25)).
apply Rplus_le_reg_l with (-sqrt (IZR beta)); ring_simplify.
apply Rle_trans with (IZR beta*(12/25-/sqrt(IZR beta))).
apply Rle_trans with (IZR beta*(71/1000)).
apply Rle_trans with (IZR 6*(71/1000)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1184 1185 1186
simpl; interval.
apply Rmult_le_compat_r.
interval.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1187
apply IZR_le; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1188 1189 1190
apply Rmult_le_compat_l.
left; apply radix_pos.
interval.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1191
assert (sqrt (IZR beta) <> 0).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1192 1193
apply Rgt_not_eq.
apply sqrt_lt_R0, radix_pos.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1194
apply Rplus_le_reg_l with (-12/25*IZR beta).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1195 1196
unfold Rdiv; ring_simplify.
right; rewrite Ropp_mult_distr_l_reverse; apply f_equal.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1197 1198
apply Rmult_eq_reg_l with (sqrt (IZR beta)); trivial.
apply trans_eq with (IZR beta).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213
field; trivial.
apply sym_eq, sqrt_sqrt.
left; apply radix_pos.
unfold Rdiv; rewrite (Rmult_assoc _ (2 - bpow (1 - prec))).
apply Rmult_le_compat_l.
left; apply radix_pos.
assert (0 <= bpow (1-prec) <= 1/32).
split.
apply bpow_ge_0.
apply Rle_trans with (bpow (-(2))).
apply bpow_le.
omega.
rewrite bpow_opp.
simpl; unfold Z.pow_pos; simpl.
rewrite Zmult_1_r.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1214
apply Rle_trans with (/IZR (6*6)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1215 1216
apply Rinv_le.
simpl; interval.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1217
apply IZR_le.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1218 1219 1220 1221 1222 1223 1224
apply Zmult_le_compat; omega.
simpl; interval.
interval.
Qed.



BOLDO Sylvie's avatar
BOLDO Sylvie committed
1225
Theorem round_flx_sqr_sqrt_aux: (4 < radix_val beta)%Z ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1226
 ((radix_val beta=5)%Z -> (3 < prec)%Z) ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1227
 (sqrt (IZR (radix_val beta)) * bpow (mag beta x-1) < x) ->
1228
  round_flx3(x/z) <= 1.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1229
Proof with auto with typeclass_instances.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1230
intros Hbeta Hprec5 H1.
Guillaume Melquiond's avatar
Guillaume Melquiond committed
1231
apply Rle_trans with (round_flx3 (x/(x-IZR k*ulp_flx x))).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243
apply round_le...
unfold Rdiv; apply Rmult_le_compat_l.
now left.
apply Rinv_le.
apply xkpos; try assumption.
apply kLe.
right; omega.
(* *)
apply round_flx_sqr_sqrt_aux1...
apply kpos.
apply kLe.
right; omega.