# Library Flocq.Calc.Sqrt

This file is part of the Flocq formalization of floating-point arithmetic in Coq: http://flocq.gforge.inria.fr/
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.

Section Fcalc_sqrt.

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 :
m1 e1,
(0 < m1)%Z
mag beta (sqrt (F2R (Float beta m1 e1))) = Z.div2 (Zdigits beta m1 + e1 + 1) :> Z.

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 :
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.

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 :
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.

End Fcalc_sqrt.