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.
Fixpoint Zsqrt_aux (p : positive) : Z × Z :=
match p with
| xH ⇒ (1, 0)%Z
| xO xH ⇒ (1, 1)%Z
| xI xH ⇒ (1, 2)%Z
| xO (xO p) ⇒
let (q, r) := Zsqrt_aux p in
let r' := (4 × r)%Z in
let d := (r' - (4 × q + 1))%Z in
if Zlt_bool d 0 then (2 × q, r')%Z else (2 × q + 1, d)%Z
| xO (xI p) ⇒
let (q, r) := Zsqrt_aux p in
let r' := (4 × r + 2)%Z in
let d := (r' - (4 × q + 1))%Z in
if Zlt_bool d 0 then (2 × q, r')%Z else (2 × q + 1, d)%Z
| xI (xO p) ⇒
let (q, r) := Zsqrt_aux p in
let r' := (4 × r + 1)%Z in
let d := (r' - (4 × q + 1))%Z in
if Zlt_bool d 0 then (2 × q, r')%Z else (2 × q + 1, d)%Z
| xI (xI p) ⇒
let (q, r) := Zsqrt_aux p in
let r' := (4 × r + 3)%Z in
let d := (r' - (4 × q + 1))%Z in
if Zlt_bool d 0 then (2 × q, r')%Z else (2 × q + 1, d)%Z
end.
Lemma Zsqrt_ind :
∀ P : positive → Prop,
P xH → P (xO xH) → P (xI xH) →
( ∀ p, P p → P (xO (xO p)) ∧ P (xO (xI p)) ∧ P (xI (xO p)) ∧ P (xI (xI p)) ) →
∀ p, P p.
Lemma Zsqrt_aux_correct :
∀ p,
let (q, r) := Zsqrt_aux p in
Zpos p = (q × q + r)%Z ∧ (0 ≤ r ≤ 2 × q)%Z.
Opaque Zmult. Transparent Zmult.
Computes the integer square root and its remainder, but
without carrying a proof, contrarily to the operation of
the standard libary.
Definition Zsqrt p :=
match p with
| Zpos p ⇒ Zsqrt_aux p
| _ ⇒ (0, 0)%Z
end.
Theorem Zsqrt_correct :
∀ x,
(0 ≤ x)%Z →
let (q, r) := Zsqrt x in
x = (q × q + r)%Z ∧ (0 ≤ r ≤ 2 × q)%Z.
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) := Zsqrt 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.