Library Flocq.Calc.Fcalc_sqrt
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2010-2013 Sylvie Boldo
Copyright (C) 2010-2013 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.
Copyright (C) 2010-2013 Guillaume Melquiond
Helper functions and theorems for computing the rounded square root of a floating-point number.
Require Import Fcore_Raux.
Require Import Fcore_defs.
Require Import Fcore_digits.
Require Import Fcore_float_prop.
Require Import Fcalc_bracket.
Require Import Fcalc_digits.
Section Fcalc_sqrt.
Variable beta : radix.
Notation bpow e := (bpow beta e).
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.
Definition Fsqrt_core prec m e :=
let d := Zdigits beta m in
let s := Zmax (2 × prec - d) 0 in
let e´ := (e - s)%Z in
let (s´, e´´) := if Zeven e´ then (s, e´) else (s + 1, e´ - 1)%Z in
let m´ :=
match s´ with
| Zpos p ⇒ (m × Zpower_pos beta p)%Z
| _ ⇒ m
end in
let (q, r) := Z.sqrtrem m´ 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, Zdiv2 e´´, l).
Theorem Fsqrt_core_correct :
∀ prec m e,
(0 < m)%Z →
let ´(m´, e´, l) := Fsqrt_core prec m e in
(prec ≤ Zdigits beta m´)%Z ∧
inbetween_float beta m´ e´ (sqrt (F2R (Float beta m e))) l.
End Fcalc_sqrt.