Library Flocq.Calc.Sqrt

This file is part of the Flocq formalization of floating-point arithmetic in Coq: https://flocq.gitlabpages.inria.fr/
Copyright (C) 2010-2018 Sylvie Boldo
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.


From Coq Require Import ZArith Reals Lia.

Require Import Zaux Raux Defs Digits Generic_fmt Float_prop Bracket.

Set Implicit Arguments.

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