Commit 7d4d7827 authored by Guillaume Melquiond's avatar Guillaume Melquiond
Browse files

Added comments on the Calc part of the libary.

parent 1965172b
......@@ -17,7 +17,9 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Locations: where is a real positioned compared to floats *)
(** * Locations: where a real number is positioned with respect to
its rounded-down value in an arbitrary format. *)
Require Import Fcore.
Section Fcalc_bracket.
......@@ -35,6 +37,8 @@ Definition inbetween_loc :=
| _ => loc_Exact
end.
(** Locates a real number with respect to the middle of two other numbers. *)
Inductive inbetween : location -> Prop :=
| inbetween_Exact : x = d -> inbetween loc_Exact
| inbetween_Inexact l : (d < x < u)%R -> Rcompare x ((d + u) / 2)%R = l -> inbetween (loc_Inexact l).
......@@ -326,6 +330,10 @@ simpl (Z2R 2).
field.
Qed.
(** Computes a new location when the interval containing a real
number is split into nb_steps subintervals and the real is
in the k-th one. (Even radix.) *)
Definition new_location_even k l :=
if Zeq_bool k 0 then
match l with loc_Exact => l | _ => loc_Inexact Lt end
......@@ -379,6 +387,10 @@ exact Hk1.
apply Hk.
Qed.
(** Computes a new location when the interval containing a real
number is split into nb_steps subintervals and the real is
in the k-th one. (Odd radix.) *)
Definition new_location_odd k l :=
if Zeq_bool k 0 then
match l with loc_Exact => l | _ => loc_Inexact Lt end
......@@ -501,6 +513,8 @@ Variable fexp : Z -> Z.
Hypothesis prop_exp : valid_exp fexp.
Notation format := (generic_format beta fexp).
(** Specialization of inbetween for two consecutive floating-point numbers. *)
Definition inbetween_float m e x l :=
inbetween (F2R (Float beta m e)) (F2R (Float beta (m + 1) e)) x l.
......@@ -520,6 +534,8 @@ now apply Rlt_le.
apply Hx.
Qed.
(** Specialization of inbetween for two consecutive integers. *)
Definition inbetween_int m x l :=
inbetween (Z2R m) (Z2R (m + 1)) x l.
......@@ -565,6 +581,8 @@ now apply inbetween_float_new_location.
apply Zmult_1_r.
Qed.
(** Relates location and rounding. *)
Theorem inbetween_float_round :
forall rnd choice,
( forall x m l, inbetween_int m x l -> Zrnd rnd x = choice m l ) ->
......@@ -582,6 +600,8 @@ apply bpow_gt_0.
now rewrite scaled_mantissa_bpow.
Qed.
(** Gives a computable way to round a real number down. *)
Theorem inbetween_float_DN :
forall x m l,
let e := canonic_exponent beta fexp x in
......@@ -604,6 +624,8 @@ Definition round_UP l :=
| _ => true
end.
(** Gives a computable way to round a real number up. *)
Theorem inbetween_float_UP :
forall x m l,
let e := canonic_exponent beta fexp x in
......@@ -637,6 +659,8 @@ Definition round_NE (p : bool) l :=
| loc_Inexact Gt => true
end.
(** Gives a computable way to round a real number to nearest even. *)
Theorem inbetween_float_NE :
forall x m l,
let e := canonic_exponent beta fexp x in
......
......@@ -17,14 +17,23 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Number of digits for "multi-precision" floats *)
Require Import Fcore.
(** * Function of computing the number of digits of integers
and related theorems. *)
Require Import Fcore_Raux.
Require Import Fcore_defs.
Require Import Fcore_float_prop.
Section Fcalc_digits.
Variable beta : radix.
Notation bpow e := (bpow beta e).
(** Computes the number of bits (radix 2) of a positive integer.
It serves as an upper bound on the number of digits to ensure termination.
*)
Fixpoint digits2_Pnat (n : positive) : nat :=
match n with
| xH => O
......@@ -205,19 +214,27 @@ apply He.
now apply (Z2R_neq _ 0).
Qed.
Theorem ln_beta_F2R_digits :
forall m e, m <> Z0 ->
(ln_beta beta (F2R (Float beta m e)) = digits m + e :> Z)%Z.
Proof.
intros m e Hm.
rewrite ln_beta_F2R with (1 := Hm).
apply (f_equal (fun v => Zplus v e)).
apply sym_eq.
now apply digits_ln_beta.
Qed.
Theorem digits_shift :
forall m e,
m <> Z0 -> (0 <= e)%Z ->
digits (m * Zpower beta e) = (digits m + e)%Z.
Proof.
intros m e Hm He.
rewrite 2!digits_ln_beta.
rewrite <- ln_beta_F2R_digits with (1 := Hm).
rewrite digits_ln_beta.
rewrite Z2R_mult.
rewrite Z2R_Zpower with (1 := He).
change (Z2R m * bpow e)%R with (F2R (Float beta m e)).
apply ln_beta_F2R.
exact Hm.
exact Hm.
now rewrite Z2R_Zpower with (1 := He).
contradict Hm.
apply Zmult_integral_l with (2 := Hm).
apply neq_Z2R.
......@@ -272,6 +289,12 @@ cut (y <= x -> digits y <= digits x)%Z. omega.
now apply digits_le.
Qed.
(** Characterizes the number digits of a product.
This strong version is needed for proofs of division and square root
algorithms, since they involve operation remainders.
*)
Theorem digits_mult_strong :
forall x y,
(0 <= x)%Z -> (0 <= y)%Z ->
......@@ -336,4 +359,21 @@ rewrite Zmult_0_l, Zplus_0_r, 2!Zplus_0_l.
apply Zle_refl.
Qed.
Theorem digits_mult :
forall x y,
(digits (x * y) <= digits x + digits y)%Z.
Proof.
intros x y.
rewrite <- digits_abs.
rewrite <- (digits_abs x).
rewrite <- (digits_abs y).
apply Zle_trans with (digits (Zabs x + Zabs y + Zabs x * Zabs y)).
apply digits_le.
apply Zabs_pos.
rewrite Zabs_Zmult.
generalize (Zabs_pos x) (Zabs_pos y).
omega.
apply digits_mult_strong ; apply Zabs_pos.
Qed.
End Fcalc_digits.
......@@ -17,8 +17,12 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Computing the rounded division *)
Require Import Fcore.
(** * Helper function and theorem for computing the rounded quotient
of two floating-point numbers. *)
Require Import Fcore_Raux.
Require Import Fcore_defs.
Require Import Fcore_float_prop.
Require Import Fcalc_bracket.
Require Import Fcalc_digits.
......@@ -27,13 +31,17 @@ Section Fcalc_div.
Variable beta : radix.
Notation bpow e := (bpow beta e).
(**
- Shift dividend mantissa so that it has at least p2 + p digits.
- Perform the euclidean division.
- Compute position with remainder.
(** Computes a mantissa of precision p, the corresponding exponent,
and the position with respect to the real quotient of the
input floating-point numbers.
The algorithm performs the following steps:
- Shift dividend mantissa so that it has at least p2 + p digits.
- Perform the Euclidean division.
- Compute the position according to the division remainder.
Complexity is fine as long as p1 <= 2p and p2 <= p.
*)
Complexity is fine as long as p1 <= 2p and p2 <= p.
*)
Definition Fdiv_aux prec m1 e1 m2 e2 :=
let d1 := digits beta m1 in
......
......@@ -17,7 +17,8 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * How to round a multi-precision float into a usable float *)
(** * Helper function for computing the rounded value of a real number. *)
Require Import Fcore.
Require Import Fcalc_bracket.
Require Import Fcalc_digits.
......@@ -31,6 +32,10 @@ Variable fexp : Z -> Z.
Hypothesis prop_exp : valid_exp fexp.
Notation format := (generic_format beta fexp).
(** Given a triple (mantissa, exponent, position), this function
computes a triple with a canonic exponent, assuming the
original triple had enough precision. *)
Definition truncate t :=
let '(m, e, l) := t in
let k := (fexp (digits beta m + e) - e)%Z in
......@@ -195,6 +200,8 @@ rewrite <- H.
now apply round_generic.
Qed.
(** Truncating a triple is sufficient to round a real number. *)
Theorem round_trunc_any_correct :
forall x m e l,
(0 <= x)%R ->
......@@ -211,6 +218,8 @@ Qed.
End round_dir.
(** * Instances of the theorems above, for the usual rounding modes. *)
Definition round_DN_correct :=
round_any_correct _ (fun m _ => m)
(fun _ => refl_equal _) (inbetween_float_DN beta fexp).
......
......@@ -17,8 +17,12 @@ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Computing the rounded square root *)
Require Import Fcore.
(** * 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_float_prop.
Require Import Fcalc_bracket.
Require Import Fcalc_digits.
......@@ -103,6 +107,10 @@ generalize (Zlt_cases (4 * r + 3 - (4 * q + 1)) 0).
case Zlt_bool ; ( split ; [ ring | omega ] ).
Qed.
(** 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
......@@ -125,17 +133,22 @@ Qed.
Variable beta : radix.
Notation bpow e := (bpow beta e).
(*
* 1. Shift the mantissa so that it has at least 2p-1 digits;
* shift it one digit more if the new exponent is not even.
* 2. Compute the square root s (at least p digits) of the new
* mantissa, and its remainder r.
* 3. Current position: r == 0 => Eq,
* r <= s => Lo,
* r >= s => Up.
*
* Complexity is fine as long as p1 <= 2p-1.
*)
(** 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_aux prec m e :=
let d := digits beta m in
......@@ -330,4 +343,4 @@ apply (Z2R_le 0).
omega.
Qed.
End Fcalc_sqrt.
\ No newline at end of file
End Fcalc_sqrt.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment