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-2011 Sylvie Boldo
Copyright (C) 2010-2011 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.


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 :
  forall P : positive -> Prop,
  P xH -> P (xO xH) -> P (xI xH) ->
  ( forall p, P p -> P (xO (xO p)) /\ P (xO (xI p)) /\ P (xI (xO p)) /\ P (xI (xI p)) ) ->
  forall p, P p.

Lemma Zsqrt_aux_correct :
  forall 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 :
  forall 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.
Complexity is fine as long as p1 <= 2p-1.

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