...
 
Commits (142)
......@@ -4,6 +4,9 @@ ChangeLog
config.log
config.status
configure
deps.dot
deps.png
deps.map
install-sh
lia.cache
Remakefile
......@@ -11,7 +14,7 @@ Remakefile
remake
remake.exe
html/
src/Flocq_version.v
src/Version.v
*.vo
*.glob
.*.aux
......
Sylvie Boldo <sylvie.boldo@inria.fr> Sylvie Boldo <sboldo@quetsche.inria.fr>
Pierre Roux <pierre@roux01.fr> Pierre Roux <proux@lri.fr>
Sylvie Boldo <sylvie.boldo@inria.fr> <sboldo@quetsche.inria.fr>
Pierre Roux <pierre.roux@onera.fr> <pierre@roux01.fr>
Pierre Roux <pierre.roux@onera.fr> <proux@lri.fr>
Prerequisites
-------------
You will need the Coq proof assistant (>= 8.4) with a Reals theory
compiled in.
The .tar.gz file is distributed with a working set of configure files. They
are not in the git repository though. Consequently, if you are building from
git, you will need autoconf (>= 2.59).
Configuring, compiling and installing
-------------------------------------
Ideally, you should just have to type:
$ ./configure && ./remake && ./remake install
The environment variable COQC can be passed to the configure script in order
to set the Coq compiler command. The configure script defaults to "coqc".
Similarly, COQDEP can be used to specify the location of "coqdep". The
COQBIN environment variable can be used to set both variables at once.
The option "--libdir=DIR" can be set to the directory where the compiled
library files should be installed by "make install". By default, the
target directory is "`$COQC -where`/user-contrib/Flocq".
The files are compiled at a logical location starting with "Flocq".
Installation instructions
=========================
Prerequisites
-------------
You will need the [Coq proof assistant](https://coq.inria.fr/) (>= 8.7)
with the `Reals` theory compiled in.
The `.tar.gz` file is distributed with a working set of configure files. They
are not in the git repository though. Consequently, if you are building from
git, you will need `autoconf` (>= 2.59).
Configuring, compiling and installing
-------------------------------------
Ideally, you should just have to type:
./configure && ./remake --jobs=2 && ./remake install
The environment variable `COQC` can be passed to the configure script in order
to set the Coq compiler command. The configure script defaults to `coqc`.
Similarly, `COQDEP` can be used to specify the location of `coqdep`. The
`COQBIN` environment variable can be used to set both variables at once.
The option `--libdir=DIR` can be set to the directory where the compiled
library files should be installed by `./remake install`. By default, the
target directory is `` `$COQC -where`/user-contrib/Flocq ``.
The files are compiled at a logical location starting with `Flocq`.
Version 2.5.2:
- ensured compatibility from Coq 8.4 to 8.6
Version 2.5.1:
- ensured compatibility with both Coq 8.4 and 8.5
Version 2.5.0:
- ensured compatibility with both Coq 8.4 and 8.5
(Flocq now provides its own version of iter_pos)
- redefined ulp, so that ulp(0) is meaningful
- renamed, generalized, and added lemmas in Fcore_ulp
- extended predecessor and successor to nonpositive values
(the previous definition of pred has been renamed pred_pos)
- removed some hypotheses on lemmas of Fprop_relative
- added more examples
. Average: proof on Sterbenz's average and correctly-rounded average
. Cody_Waite: Cody & Waite's approximation of exponential
. Compute: effective FP computations with an example of sqrt(sqr(x))
in radix 5 and precision 3
. Division_u16: integer division using floating-point FMA
. Triangle: Kahan's algorithm for the area of a triangle
Version 2.4.0:
- moved some lemmas from Fcalc_digits to Fcore_digits and made them axiom-free
- added theorems about double rounding being innocuous (Fappli_double_round.v)
- added example about double rounding in odd radix
- improved a bit the efficiency of IEEE-754 arithmetic
Version 2.3.0:
- moved some lemmas from Fcalc_digits to Fcore_digits and made them axiom-free
- used the square root from the standard library instead of a custom one
- added an example about sqrt(sqr(x))
Version 2.2.2:
- fixed install target for case-insensitive filesystems
Version 2.2.1:
- fixed regeneration of Flocq_version.v
Version 2.2.0:
- added theorems about rounding to odd and double rounding
- improved handling of special values of IEEE-754 arithmetic
Version 2.1.0:
- ensured compatibility with both Coq 8.3 and 8.4
- improved support for rounding toward and away from zero
- proved that formats are stable by arbitrary rounding or ulp addition
- generalized usage of ln_beta
Version 2.0.0:
- changed rounding modes from records to R -> Z functions as arguments
to the round function
. rndDN -> Zfloor
. rndUP -> Zceil
. rndZE -> Ztrunc
. rndNE -> ZnearestE
. rndNA -> ZnearestA
- added typeclasses for automatic inference of properties:
. Valid_exp for Z -> Z functions describing formats
. Valid_rnd for R -> Z functions describing rounding modes
. Exp_not_FTZ for Z -> Z functions describing formats with subnormal
. Monotone_exp for Z -> Z functions describing regular formats
. Exists_NE for radix/formats supporting rounding to nearest even
- removed theorems superseded by typeclasses:
. FIX_exp_correct, FLX_exp_correct, FLT_exp_correct, FTZ_exp_correct
. FIX_exp_monotone, FLX_exp_monotone, FLT_monotone
. Rnd_NE_pt_FIX, Rnd_NE_pt_FLX, Rnd_NE_pt_FLT
. round_NE_pt_FLX, round_NE_pt_FLT
. Zrnd_opp_le, Zrnd_Z2R_opp
- removed example file Fappli_sqrt_FLT_ne
- split theorems on equivalence between specific and generic formats
into both directions:
. FIX_format_generic and generic_format_FIX
. FLX_format_generic and generic_format_FLX
. FLT_format_generic and generic_format_FLT
. FTZ_format_generic and generic_format_FTZ
- modified correctness theorems for IEEE-754 operators so that they also
mention finiteness of the results
- added a Flocq_version.Flocq_version Coq variable for Flocq detection
and version testing in configure scripts
- moved parts of some files into other files:
. Fcore_Raux -> Fcore_Zaux
. Fcalc_digits -> Fcore_digits
. Fappli_IEEE -> Fappli_IEEE_bits
- renamed functions:
. canonic_exponent -> canonic_exp
. digits -> Zdigits
. bounded_prec -> canonic_mantissa
. Bsign_FF -> sign_FF
. shl -> shl_align
. shl_fexp -> shl_align_fexp
. binary_round_sign -> binary_round_aux
. binary_round_sign_shl -> binary_round
- renamed theorems more uniformly:
. Rabs_Rminus_pos -> Rabs_minus_le
. exp_increasing_weak -> exp_le
. ln_beta_monotone -> ln_beta_le
. ln_beta_monotone_abs -> ln_beta_le_abs
. abs_F2R -> F2R_Zabs (changed direction)
. opp_F2R -> F2R_Zopp (changed direction)
. scaled_mantissa_bpow -> scaled_mantissa_mult_bpow
. round_monotone -> round_le
. round_monotone_l -> round_ge_generic
. round_monotone_r -> round_le_generic
. round_monotone_abs_l -> abs_round_ge_generic
. round_monotone_abs_r -> abs_round_le_generic
. generic_format_canonic_exponent -> generic_format_F2R (modified hypothesis)
. canonic_exponent_round -> canonic_exp_round_ge
. generic_N_pt -> round_N_pt
. round_pred_pos_imp_rnd -> round_pred_ge_0
. round_pred_rnd_imp_pos -> round_pred_gt_0
. round_pred_neg_imp_rnd -> round_pred_le_0
. round_pred_rnd_imp_neg -> round_pred_lt_0
. ulp_le_pos -> ulp_le_id
. succ_lt_le -> succ_le_lt
. ulp_monotone -> ulp_le
. ulp_DN_pt_eq -> ulp_DN
. format_pred -> generic_format_pred
. pred_ulp -> pred_plus_ulp
. pred_lt -> pred_lt_id
. FLT_generic_format_FLX -> generic_format_FLT_FLX
. FLX_generic_format_FLT -> generic_format_FLX_FLT
. FIX_generic_format_FLT -> generic_format_FIX_FLT
. FLT_generic_format_FIX -> generic_format_FLT_FIX
. FLT_round_FLX -> round_FLT_FLX
. FTZ_round -> round_FTZ_FLX
. FTZ_round_small -> round_FTZ_small
. FLT_canonic_FLX -> canonic_exp_FLT_FLX
. FLT_canonic_FIX -> canonic_exp_FLT_FIX
. canonic_exponent_opp -> canonic_exp_opp
. canonic_exponent_abs -> canonic_exp_abs
. canonic_exponent_fexp -> canonic_exp_fexp
. canonic_exponent_fexp_pos -> canonic_exp_fexp_pos
. canonic_exponent_DN -> canonic_exp_DN
. canonic_exp_ge -> abs_lt_bpow_prec (modified hypotheses)
. Fopp_F2R -> F2R_opp
. Fabs_F2R -> F2R_abs
. plus_F2R -> F2R_plus
. minus_F2R -> F2R_minus
. mult_F2R -> F2R_mult
. digits_abs -> Zdigits_abs
. digits_ge_0 -> Zdigits_ge_0
. digits_gt_0 -> Zdigits_gt_0
. ln_beta_F2R_Zdigits -> ln_beta_F2R_Zdigits
. digits_shift -> Zdigits_mult_Zpower
. digits_Zpower -> Zdigits_Zpower
. digits_le -> Zdigits_le
. lt_digits -> lt_Zdigits
. Zpower_le_digits -> Zpower_le_Zdigits
. digits_le_Zpower -> Zdigits_le_Zpower
. Zpower_gt_digits -> Zpower_gt_Zdigits
. digits_gt_Zpower -> Zdigits_gt_Zpower
. digits_mult_strong -> Zdigits_mult_strong
. digits_mult -> Zdigits_mult
. digits_mult_ge -> Zdigits_mult_ge
. digits_shr -> Zdigits_div_Zpower
. format_add -> generic_format_plus_prec (modified hypothesis)
. format_nx -> ex_Fexp_canonic
. generic_relative_error -> relative_error
. generic_relative_error_ex -> relative_error_ex
. generic_relative_error_F2R -> relative_error_F2R_emin
. generic_relative_error_F2R_ex -> relative_error_F2R_emin_ex
. generic_relative_error_2 -> relative_error_round
. generic_relative_error_F2R_2 -> relative_error_round_F2R_emin
. generic_relative_error_N -> relative_error_N
. generic_relative_error_N_ex -> relative_error_N_ex
. generic_relative_error_N_F2R -> relative_error_N_F2R_emin
. generic_relative_error_N_F2R_ex -> relative_error_N_F2R_emin_ex
. generic_relative_error_N_2 -> relative_error_N_round
. generic_relative_error_N_F2R_2 -> relative_error_N_round_F2R_emin
. relative_error_FLT_F2R -> relative_error_FLT_F2R_emin
. relative_error_FLT_F2R_ex -> relative_error_FLT_F2R_emin_ex
. relative_error_N_FLT_2 -> relative_error_N_FLT_round
. relative_error_N_FLT_F2R -> relative_error_N_FLT_F2R_emin
. relative_error_N_FLT_F2R_ex -> relative_error_N_FLT_F2R_emin_ex
. relative_error_N_FLT_F2R_2 -> relative_error_N_FLT_round_F2R_emin
. relative_error_FLX_2 -> relative_error_FLX_round
. relative_error_N_FLX_2 -> relative_error_N_FLX_round
. canonic_bounded_prec -> canonic_canonic_mantissa
. B2R_lt_emax -> abs_B2R_lt_emax
. binary_unicity -> B2FF_inj
. finite_binary_unicity -> B2R_inj
. binary_round_sign_correct -> binary_round_aux_correct
. shl_correct -> shl_align_correct
. snd_shl -> snd_shl_align
. shl_fexp_correct -> shl_align_fexp_correct
. binary_round_sign_shl_correct -> binary_round_correct
Version 1.4.0:
- improved efficiency of IEEE-754 addition
- fixed compilation with Coq 8.3
Version 1.3:
- fixed overflow handling in IEEE-754 formats with directed rounding
Version 1.2:
- added IEEE-754 binary32 and 64 formats, including infinities and NaN
Version 1.1:
- simplified effective rounding for negative reals
- proved monotonicity of exponent functions for common formats
Version 1.0:
- initial release
This diff is collapsed.
The Flocq library provides vernacular files formalizing multi-radix
multi-precision fixed- and floating-point arithmetic for the Coq proof
assistant.
This package is free software; you can redistribute it and/or modify it
under the terms of GNU Lesser General Public License (see the COPYING
file). Authors are Sylvie Boldo <sylvie.boldo@inria.fr> and Guillaume
Melquiond <guillaume.melquiond@inria.fr>.
FLOCQ
=====
The Flocq library provides vernacular files formalizing multi-radix
multi-precision fixed- and floating-point arithmetic for the
[Coq proof assistant](https://coq.inria.fr/).
PROJECT HOME
------------
Homepage: http://flocg.gforge.inria.fr/
Repository: https://gitlab.inria.fr/flocq/flocq
Bug tracker: https://gitlab.inria.fr/flocq/flocq/issues
COPYRIGHT
---------
This package is free software; you can redistribute it and/or modify it
under the terms of GNU Lesser General Public License (see the
[COPYING](COPYING) file). Authors are Sylvie Boldo <sylvie.boldo@inria.fr>
and Guillaume Melquiond <guillaume.melquiond@inria.fr>.
INSTALLATION
------------
See the file [INSTALL.md](INSTALL.md).
FILES = \
Flocq_version.v \
Core/Fcore_Raux.v \
Core/Fcore_Zaux.v \
Core/Fcore_defs.v \
Core/Fcore_digits.v \
Core/Fcore_float_prop.v \
Core/Fcore_FIX.v \
Core/Fcore_FLT.v \
Core/Fcore_FLX.v \
Core/Fcore_FTZ.v \
Core/Fcore_generic_fmt.v \
Core/Fcore_rnd.v \
Core/Fcore_rnd_ne.v \
Core/Fcore_ulp.v \
Core/Fcore.v \
Calc/Fcalc_bracket.v \
Calc/Fcalc_digits.v \
Calc/Fcalc_div.v \
Calc/Fcalc_ops.v \
Calc/Fcalc_round.v \
Calc/Fcalc_sqrt.v \
Prop/Fprop_div_sqrt_error.v \
Prop/Fprop_mult_error.v \
Prop/Fprop_plus_error.v \
Prop/Fprop_relative.v \
Prop/Fprop_Sterbenz.v \
Appli/Fappli_rnd_odd.v \
Appli/Fappli_IEEE.v \
Appli/Fappli_IEEE_bits.v \
Appli/Fappli_double_round.v
Version.v \
Core/Raux.v \
Core/Zaux.v \
Core/Defs.v \
Core/Digits.v \
Core/Float_prop.v \
Core/FIX.v \
Core/FLT.v \
Core/FLX.v \
Core/FTZ.v \
Core/Generic_fmt.v \
Core/Round_pred.v \
Core/Round_NE.v \
Core/Ulp.v \
Core/Core.v \
Calc/Bracket.v \
Calc/Div.v \
Calc/Operations.v \
Calc/Round.v \
Calc/Sqrt.v \
Prop/Div_sqrt_error.v \
Prop/Mult_error.v \
Prop/Plus_error.v \
Prop/Relative.v \
Prop/Sterbenz.v \
Prop/Round_odd.v \
Prop/Double_rounding.v \
IEEE754/Binary.v \
IEEE754/Bits.v \
Pff/Pff.v \
Pff/Pff2FlocqAux.v \
Pff/Pff2Flocq.v
OBJS = $(addprefix src/,$(addsuffix o,$(FILES)))
EXAMPLES = \
Average.v \
Compute.v \
Double_round_beta_odd.v \
Double_rounding_odd_radix.v \
Homogen.v
MORE_EXAMPLES = \
......@@ -58,8 +60,8 @@ check-more: $(MOBJS)
Remakefile: Remakefile.in config.status
./config.status Remakefile
src/Flocq_version.v: src/Flocq_version.v.in config.status
./config.status src/Flocq_version.v
src/Version.v: src/Version.v.in config.status
./config.status src/Version.v
configure config.status: configure.in
autoconf
......@@ -69,9 +71,13 @@ configure config.status: configure.in
@COQDEP@ -R src Flocq $< | @REMAKE@ -r $@
@COQC@ -R src Flocq $<
examples/%.vo: examples/%.v
@COQDEP@ -R src Flocq -R examples FlocqEx $< | @REMAKE@ -r $@
@COQC@ -R src Flocq -R examples FlocqEx $<
clean:
rm -f $(OBJS) $(EOBJS) $(MOBJS) src/*.glob examples/*.glob
for d in src src/Core src/Calc src/Prop src/Appli examples; do \
for d in src src/Core src/Calc src/Prop src/Pff src/IEEE754 examples; do \
rm -f $d/.coq-native/*.o $d/.coq-native/*.cm*; done
find . -type d -name ".coq-native" -empty -prune -exec rmdir "{}" \;
......@@ -79,8 +85,39 @@ html/index.html: $(OBJS)
rm -rf html
mkdir -p html
@COQDOC@ -toc -interpolate -utf8 -html -g -R src Flocq -d html \
--coqlib http://coq.inria.fr/distrib/current/stdlib/ \
--coqlib http://coq.inria.fr/distrib/current/stdlib \
$(addprefix src/,$(FILES))
for f in html/*.html; do
sed -e 's;<a href="index.html">Index</a>;Go back to the <a href="../index.html">Main page</a> or <a href="index.html">Index</a>.;' -i $f
done
deps.dot: $(addprefix src/,$(FILES)) Remakefile.in
(echo "digraph flocq_deps { pack=true; rank=max;"
echo "node [shape=ellipse, style=filled, URL=\"html/Flocq.\N.html\", color=black];"
echo '{ rank=same; "Core.Zaux"; "Core.Raux"; }'
echo '{ rank=same; "Core.FLX"; "Core.FIX"; "Core.Round_NE"; "Calc.Operations"; }'
echo '{ rank=same; "Core.FLT"; "Core.FTZ"; }'
echo '{ rank=same; "Core.Generic_fmt"; "Core.Ulp"; }'
echo '{ rank=same; "IEEE754.Binary"; "IEEE754.Bits"; }'
echo '{ rank=same; "Pff.Pff2FlocqAux"; "Pff.Pff2Flocq"; }'
(cd src ; @COQDEP@ -R . Flocq $(FILES)) |
sed -n -e 's,/,.,g;s/[.]vo.*: [^ ]*[.]v//p' |
while read src dst; do
color=$$(echo "$src" | sed -e 's,Core.*,turquoise,;s,Calc.*,plum,;s,Prop.*,lightcoral,;s,Pff.*,yellow,;s,IEEE754.*,cornflowerblue,;s,Version.*,white,')
echo "\"$src\" [fillcolor=$color];"
for d in $dst; do
echo "\"$src\" -> \"${d%.vo}\" ;"
done
done;
echo "}") | tred > $@
deps.png: deps.dot
dot -T png deps.dot > deps.png
deps.map: deps.dot
dot -T cmap deps.dot | sed -e 's,>$,/>,' > deps.map
doc: html/index.html
......@@ -88,16 +125,15 @@ install:
prefix=@prefix@
exec_prefix=@exec_prefix@
mkdir -p @libdir@
for d in Core Calc Prop Appli; do mkdir -p @libdir@/$d; done
for d in Core Calc Prop IEEE754 Pff; do mkdir -p @libdir@/$d; done
for f in $(OBJS); do cp $f @libdir@/${f#src/}; done
for f in $(FILES); do cp src/$f @libdir@/$f; done
( cd src && find . -type d -name ".coq-native" -exec cp -RT "{}" "@libdir@/{}" \; )
EXTRA_DIST = \
configure
REMOVE_FROM_DIST = \
src/Appli/Fappli_Axpy.v \
src/Translate/
REMOVE_FROM_DIST =
dist: $(EXTRA_DIST)
PACK=@PACKAGE_TARNAME@-@PACKAGE_VERSION@
......
-R src Flocq
AC_INIT([Flocq], [2.5.2],
AC_INIT([Flocq], [3.1.0],
[Sylvie Boldo <sylvie.boldo@inria.fr>, Guillaume Melquiond <guillaume.melquiond@inria.fr>],
[flocq])
......@@ -68,5 +68,5 @@ MINGW*)
;;
esac
AC_CONFIG_FILES([Remakefile src/Flocq_version.v])
AC_CONFIG_FILES([Remakefile src/Version.v])
AC_OUTPUT
This diff is collapsed.
Require Import Reals Fcore.
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2014-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.
*)
Require Import Reals Flocq.Core.Core.
Require Import Gappa.Gappa_tactic Interval.Interval_tactic.
Open Scope R_scope.
Lemma rel_helper :
forall xa xe b : R,
xe <> 0 ->
(Rabs ((xa - xe) / xe) <= b <-> Rabs (xa - xe) <= b * Rabs xe).
Proof.
intros xa xe b Zx.
unfold Rdiv.
rewrite Rabs_mult, Rabs_Rinv by exact Zx.
split ; intros H.
- apply Rmult_le_reg_r with (/ Rabs xe).
apply Rinv_0_lt_compat.
now apply Rabs_pos_lt.
rewrite Rmult_assoc, Rinv_r, Rmult_1_r.
exact H.
now apply Rabs_no_R0.
- apply Rmult_le_reg_r with (Rabs xe).
now apply Rabs_pos_lt.
rewrite Rmult_assoc, Rinv_l, Rmult_1_r.
exact H.
now apply Rabs_no_R0.
Qed.
Lemma Rdiv_compat_r :
forall a b c : R, a <> 0 -> c <> 0 ->
(a*b) / (a*c) = b/c.
Proof.
intros a b c Ha Hc.
field.
apply conj.
apply Hc.
apply Ha.
Qed.
Notation pow2 := (bpow radix2).
Definition mul x y := (round radix2 (FLT_exp (-1074) 53) ZnearestE (x * y)).
Definition sub x y := (round radix2 (FLT_exp (-1074) 53) ZnearestE (x - y)).
Notation rnd := (round radix2 (FLT_exp (-1074) 53) ZnearestE).
Definition add x y := rnd (x + y).
Definition sub x y := rnd (x - y).
Definition mul x y := rnd (x * y).
Definition div x y := rnd (x / y).
Definition nearbyint x := round radix2 (FIX_exp 0) ZnearestE x.
Definition Log2h := 3048493539143 * pow2 (-42).
Definition Log2l := 544487923021427 * pow2 (-93).
Definition InvLog2 := 3248660424278399 * pow2 (-51).
Definition p0 := 1 * pow2 (-2).
Definition p1 := 4002712888408905 * pow2 (-59).
Definition p2 := 1218985200072455 * pow2 (-66).
Definition q0 := 1 * pow2 (-1).
Definition q1 := 8006155947364787 * pow2 (-57).
Definition q2 := 4573527866750985 * pow2 (-63).
Definition cw_exp (x : R) :=
let k := nearbyint (mul x InvLog2) in
let t := sub (sub x (mul k Log2h)) (mul k Log2l) in
let t2 := mul t t in
let p := add p0 (mul t2 (add p1 (mul t2 p2))) in
let q := add q0 (mul t2 (add q1 (mul t2 q2))) in
pow2 (Zfloor k + 1) * (add (div (mul t p) (sub q (mul t p))) (1/2)).
Lemma method_error :
forall t : R,
let t2 := t * t in
let p := p0 + t2 * (p1 + t2 * p2) in
let q := q0 + t2 * (q1 + t2 * q2) in
let f := 2 * ((t * p) / (q - t * p) + 1/2) in
Rabs t <= 355 / 1024 ->
Rabs ((f - exp t) / exp t) <= 23 * pow2 (-62).
Proof.
intros t t2 p q f Ht.
unfold f, q, p, t2, p0, p1, p2, q0, q1, q2 ;
interval with (i_bisect_taylor t 9, i_prec 70).
Qed.
Lemma argument_reduction :
forall x : R,
generic_format radix2 (FLT_exp (-1074) 53) x ->
Rabs x <= 710 ->
let k := round radix2 (FIX_exp 0) ZnearestE (mul x InvLog2) in
-746 <= x <= 710 ->
let k := nearbyint (mul x InvLog2) in
let tf := sub (sub x (mul k Log2h)) (mul k Log2l) in
let te := x - k * ln 2 in
Rabs tf <= 355 / 1024 /\
Rabs (tf - te) <= 65537 * pow2 (-71).
Proof with auto with typeclass_instances.
intros x Fx Bx k tf te.
assert (Rabs x <= 5/16 \/ 5/16 <= Rabs x <= 710) as [Bx'|Bx'] by gappa.
assert (Rabs x <= 5/16 \/ 5/16 <= Rabs x <= 746) as [Bx'|Bx'] by gappa.
- assert (k = 0).
clear -Bx'.
refine (let H := _ in Rle_antisym _ _ (proj2 H) (proj1 H)).
unfold k, mul, InvLog2, Log2h ; gappa.
unfold k, mul, nearbyint, InvLog2, Log2h ; gappa.
unfold tf, te, mul, sub.
rewrite H.
rewrite 2!Rmult_0_l.
......@@ -36,13 +120,58 @@ assert (Rabs x <= 5/16 \/ 5/16 <= Rabs x <= 710) as [Bx'|Bx'] by gappa.
rewrite 2!round_generic with (2 := Fx)...
gappa.
- assert (Hl: - 1 * pow2 (-102) <= Log2l - (ln 2 - Log2h) <= 0).
unfold Log2l, Log2h ; simpl bpow ;
unfold Log2l, Log2h ;
interval with (i_prec 110).
assert (Ax: x = x * InvLog2 * (1 / InvLog2)).
field.
unfold InvLog2 ; simpl ; interval.
unfold InvLog2 ; interval.
unfold te.
replace (x - k * ln 2) with (x - k * Log2h - k * (ln 2 - Log2h)) by ring.
revert Hl Ax.
unfold tf, te, k, sub, mul, InvLog2, Log2h, Log2l ; gappa.
unfold tf, te, k, sub, mul, nearbyint, InvLog2, Log2h, Log2l ; gappa.
Qed.
Theorem exp_correct :
forall x : R,
generic_format radix2 (FLT_exp (-1074) 53) x ->
-746 <= x <= 710 ->
Rabs ((cw_exp x - exp x) / exp x) <= 1 * pow2 (-51).
Proof.
intros x Fx Bx.
generalize (argument_reduction x Fx Bx).
unfold cw_exp.
set (k := nearbyint (mul x InvLog2)).
set (t := sub (sub x (mul k Log2h)) (mul k Log2l)).
set (t2 := mul t t).
intros [Bt Et].
clearbody t.
generalize (method_error t Bt).
intros Ef.
rewrite bpow_plus, Rmult_assoc.
assert (exp x = pow2 (Zfloor k) * exp (x - k * ln 2)) as ->.
assert (exists k', k = IZR k') as [k' ->] by (eexists ; apply Rmult_1_r).
rewrite Zfloor_IZR, bpow_exp, <- exp_plus.
apply f_equal.
simpl ; ring.
rewrite <- Rmult_minus_distr_l.
rewrite Rdiv_compat_r.
2: apply Rgt_not_eq, bpow_gt_0.
2: apply Rgt_not_eq, exp_pos.
clearbody k.
revert Et.
generalize (x - k * ln 2).
clear x Fx Bx.
intros x Ex.
assert (Rabs ((exp t - exp x) / exp x) <= 33 * pow2 (-60)).
unfold Rdiv.
rewrite Rmult_minus_distr_r.
rewrite Rinv_r by apply Rgt_not_eq, exp_pos.
rewrite <- exp_Ropp, <- exp_plus.
revert Ex.
unfold Rminus ; generalize (t + - x) ; clear.
intros r Hr ; interval with (i_prec 60).
apply rel_helper. apply Rgt_not_eq, exp_pos.
apply rel_helper in H. 2: apply Rgt_not_eq, exp_pos.
apply rel_helper in Ef. 2: apply Rgt_not_eq, exp_pos.
unfold t2, add, mul, sub, div, p0, p1, p2, q0, q1, q2 in * ; gappa.
Qed.
Require Import Flocq.Core.Fcore.
Require Import Flocq.Calc.Fcalc_bracket Flocq.Calc.Fcalc_round Flocq.Calc.Fcalc_ops Flocq.Calc.Fcalc_div Flocq.Calc.Fcalc_sqrt.
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2015-2018 Sylvie Boldo
#<br />#
Copyright (C) 2015-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.
*)
From Flocq Require Import Core Bracket Round Operations Div Sqrt.
Section Compute.
......@@ -29,15 +47,15 @@ Proof.
intros m e.
case Zlt_bool_spec ; intros H.
apply Rlt_bool_true.
now apply F2R_lt_0_compat.
now apply F2R_lt_0.
apply Rlt_bool_false.
now apply F2R_ge_0_compat.
now apply F2R_ge_0.
Qed.
Definition plus (x y : float beta) :=
let (m, e) := Fplus beta x y in
let (m, e) := Fplus x y in
let s := Zlt_bool m 0 in
let '(m', e', l) := truncate beta fexp (Zabs m, e, loc_Exact) in
let '(m', e', l) := truncate beta fexp (Z.abs m, e, loc_Exact) in
Float beta (cond_Zopp s (choice s m' l)) e'.
Theorem plus_correct :
......@@ -47,8 +65,8 @@ Proof.
intros x y.
unfold plus.
rewrite <- F2R_plus.
destruct (Fplus beta x y) as [m e].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ (Zabs m) e loc_Exact).
destruct (Fplus x y) as [m e].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ (Z.abs m) e loc_Exact).
3: now right.
destruct truncate as [[m' e'] l'].
apply (f_equal (fun s => F2R (Float beta (cond_Zopp s (choice s _ _)) _))).
......@@ -58,9 +76,9 @@ apply sym_eq, F2R_Zabs.
Qed.
Definition mult (x y : float beta) :=
let (m, e) := Fmult beta x y in
let (m, e) := Fmult x y in
let s := Zlt_bool m 0 in
let '(m', e', l) := truncate beta fexp (Zabs m, e, loc_Exact) in
let '(m', e', l) := truncate beta fexp (Z.abs m, e, loc_Exact) in
Float beta (cond_Zopp s (choice s m' l)) e'.
Theorem mult_correct :
......@@ -70,8 +88,8 @@ Proof.
intros x y.
unfold mult.
rewrite <- F2R_mult.
destruct (Fmult beta x y) as [m e].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ (Zabs m) e loc_Exact).
destruct (Fmult x y) as [m e].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ (Z.abs m) e loc_Exact).
3: now right.
destruct truncate as [[m' e'] l'].
apply (f_equal (fun s => F2R (Float beta (cond_Zopp s (choice s _ _)) _))).
......@@ -81,9 +99,8 @@ apply sym_eq, F2R_Zabs.
Qed.
Definition sqrt (x : float beta) :=
let (m, e) := x in
if Zlt_bool 0 m then
let '(m', e', l) := truncate beta fexp (Fsqrt_core beta prec m e) in
if Zlt_bool 0 (Fnum x) then
let '(m', e', l) := truncate beta fexp (Fsqrt fexp x) in
Float beta (choice false m' l) e'
else Float beta 0 0.
......@@ -91,20 +108,18 @@ Theorem sqrt_correct :
forall x : float beta,
round beta fexp rnd (R_sqrt.sqrt (F2R x)) = F2R (sqrt x).
Proof.
intros [m e].
intros x.
unfold sqrt.
case Zlt_bool_spec ; intros Hm.
generalize (Fsqrt_core_correct beta prec m e Hm).
destruct Fsqrt_core as [[m' e'] l].
generalize (Fsqrt_correct fexp x (F2R_gt_0 _ _ Hm)).
destruct Fsqrt as [[m' e'] l].
intros [Hs1 Hs2].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ m' e' l).
rewrite (round_trunc_sign_any_correct' beta fexp rnd choice rnd_choice _ m' e' l).
destruct truncate as [[m'' e''] l'].
now rewrite Rlt_bool_false by apply sqrt_ge_0.
now rewrite Rabs_pos_eq by apply sqrt_ge_0.
left.
apply Zle_trans with (2 := fexp_prec _).
clear -Hs1 ; omega.
destruct (Req_dec (F2R (Float beta m e)) 0) as [Hz|Hz].
now left.
destruct (Req_dec (F2R x) 0) as [Hz|Hz].
rewrite Hz, sqrt_0, F2R_0.
now apply round_0.
unfold R_sqrt.sqrt.
......@@ -113,13 +128,13 @@ rewrite F2R_0.
now apply round_0.
destruct H as [H|H].
elim Rgt_not_le with (1 := H).
now apply F2R_le_0_compat.
now apply F2R_le_0.
now elim Hz.
Qed.
Lemma Rsgn_div :
forall x y : R,
x <> R0 -> y <> R0 ->
x <> 0%R -> y <> 0%R ->
Rlt_bool (x / y) 0 = xorb (Rlt_bool x 0) (Rlt_bool y 0).
Proof.
intros x y Hx0 Hy0.
......@@ -162,51 +177,54 @@ now elim Hy0.
Qed.
Definition div (x y : float beta) :=
let (mx, ex) := x in
let (my, ey) := y in
if Zeq_bool mx 0 then Float beta 0 0
if Zeq_bool (Fnum x) 0 then Float beta 0 0
else
let '(m, e, l) := truncate beta fexp (Fdiv_core beta prec (Zabs mx) ex (Zabs my) ey) in
let s := xorb (Zlt_bool mx 0) (Zlt_bool my 0) in
let '(m, e, l) := truncate beta fexp (Fdiv fexp (Fabs x) (Fabs y)) in
let s := xorb (Zlt_bool (Fnum x) 0) (Zlt_bool (Fnum y) 0) in
Float beta (cond_Zopp s (choice s m l)) e.
Theorem div_correct :
forall x y : float beta,
F2R y <> R0 ->
F2R y <> 0%R ->
round beta fexp rnd (F2R x / F2R y) = F2R (div x y).
Proof.
intros [mx ex] [my ey] Hy.
intros x y Hy.
unfold div.
case Zeq_bool_spec ; intros Hm.
rewrite Hm, 2!F2R_0.
unfold Rdiv.
rewrite Rmult_0_l.
now apply round_0.
assert (Hx: (0 < Zabs mx)%Z) by now apply Z.abs_pos.
assert (Hy': (0 < Zabs my)%Z).
apply Z.abs_pos.
destruct x as [mx ex].
simpl in Hm.
rewrite Hm, 2!F2R_0.
unfold Rdiv.
rewrite Rmult_0_l.
now apply round_0.
assert (0 < F2R (Fabs x))%R as Hx.
destruct x as [mx ex].
now apply F2R_gt_0, Z.abs_pos.
assert (0 < F2R (Fabs y))%R as Hy'.
destruct y as [my ey].
apply F2R_gt_0, Z.abs_pos.
contradict Hy.
rewrite Hy.
apply F2R_0.
generalize (Fdiv_core_correct beta prec (Zabs mx) ex (Zabs my) ey Hprec Hx Hy').
destruct Fdiv_core as [[m e] l].
generalize (Fdiv_correct fexp _ _ Hx Hy').
destruct Fdiv as [[m e] l].
intros [Hs1 Hs2].
rewrite (round_trunc_sign_any_correct beta fexp rnd choice rnd_choice _ m e l).
destruct truncate as [[m' e'] l'].
apply (f_equal (fun s => F2R (Float beta (cond_Zopp s (choice s _ _)) _))).
rewrite Rsgn_div.
apply f_equal2 ; apply Rsgn_F2R.
contradict Hm.
apply F2R_eq_0_reg with (1 := Hm).
exact Hy.
unfold Rdiv.
rewrite Rabs_mult, Rabs_Rinv.
rewrite <- 2!F2R_Zabs.
exact Hs2.
exact Hy.
left.
apply Zle_trans with (2 := fexp_prec _).
clear -Hs1 ; omega.
rewrite (round_trunc_sign_any_correct' beta fexp rnd choice rnd_choice _ m e l).
- destruct truncate as [[m' e'] l'].
apply (f_equal (fun s => F2R (Float beta (cond_Zopp s (choice s _ _)) _))).
rewrite Rsgn_div.
apply f_equal2 ; apply Rsgn_F2R.
contradict Hm.
apply eq_0_F2R with (1 := Hm).
exact Hy.
- unfold Rdiv.
rewrite Rabs_mult, Rabs_Rinv by exact Hy.
now rewrite <- 2!F2R_abs.
- left.
rewrite <- cexp_abs.
unfold Rdiv.
rewrite Rabs_mult, Rabs_Rinv by exact Hy.
now rewrite <- 2!F2R_abs.
Qed.
End Compute.
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2014-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.
*)
Require Import Reals Psatz.
Require Import Fcore Gappa.Gappa_tactic.
Require Import Flocq.Core.Core Gappa.Gappa_tactic.
Open Scope R_scope.
......@@ -16,12 +33,11 @@ field.
now apply Rgt_not_eq.
Qed.
Lemma FIX_format_Z2R :
forall beta x, FIX_format beta 0 (Z2R x).
Lemma FIX_format_IZR :
forall beta x, FIX_format beta 0 (IZR x).
Proof.
intros beta x.
exists (Float beta x 0).
split.
apply sym_eq, Rmult_1_r.
apply eq_refl.
Qed.
......@@ -48,15 +64,15 @@ now apply Rplus_le_reg_r with a.
now apply Rplus_lt_reg_r with a.
Qed.
Coercion Z2R : Z >-> R.
Coercion IZR : Z >-> R.
Lemma Z2R_le_le :
Lemma IZR_le_le :
forall a b c,
(a <= b <= c)%Z ->
a <= b <= c.
Proof.
intros a b c.
now split ; apply Z2R_le.
now split ; apply IZR_le.
Qed.
Notation pow2 := (bpow radix2).
......@@ -89,29 +105,26 @@ set (q0 := fma a y0 0).
set (e0 := fnma b y0 (1 + pow2 (-17))).
set (q1 := fma e0 q0 q0).
apply Zfloor_imp.
rewrite Z2R_plus.
rewrite plus_IZR.
simpl.
apply Rmult_le_lt_reg_l with b.
apply (Z2R_lt 0) ; lia.
apply IZR_lt ; lia.
apply Rplus_le_lt_reg_r with (-a).
replace (b * (a / b)%Z + - a) with (-(a - (a / b)%Z * b)) by ring.
replace (b * ((a / b)%Z + 1) + - a) with (b - (a - (a / b)%Z * b)) by ring.
rewrite <- Z2R_mult, <- Z2R_minus.
rewrite <- mult_IZR, <- minus_IZR.
rewrite <- Zmod_eq_full by lia.
cut (0 <= b * q1 - a < 1).
cut (0 <= Zmod a b <= b - 1).
lra.
change 1 with (Z2R 1).
rewrite <- Z2R_minus.
apply (Z2R_le_le 0).
rewrite <- minus_IZR.
apply IZR_le_le.
generalize (Z_mod_lt a b).
lia.
assert (Ba': 1 <= a <= 65535).
change (1%Z <= a <= 65535%Z).
now split ; apply Z2R_le.
now split ; apply IZR_le.
assert (Bb': 1 <= b <= 65535).
change (1%Z <= b <= 65535%Z).
now split ; apply Z2R_le.
now split ; apply IZR_le.
refine (let '(conj H1 H2) := _ in conj H1 (Rnot_le_lt _ _ H2)).
set (err := (q1 - a / b) / (a / b)).
replace (b * q1 - a) with (a * err) by abstract (unfold err ; field ; lra).
......@@ -123,6 +136,6 @@ assert (H: (Mq1 - a / b) / (a / b) = -(eps0 * eps0) + (1 + eps0) * pow2 (-17))
by abstract (unfold Mq1, Me0, Mq0, eps0 ; field ; lra).
revert H.
unfold Mq1, Me0, Mq0, eps0, err, q1, e0, q0, y0.
generalize (frcpa_spec b) (FIX_format_Z2R radix2 a) (FIX_format_Z2R radix2 b).
generalize (frcpa_spec b) (FIX_format_IZR radix2 a) (FIX_format_IZR radix2 b).
gappa.
Qed.
Require Import Reals.
Require Import Fcore.
Require Import Fcalc_ops.
Require Import Psatz.
Require Import Fprop_relative.
Require Import Fprop_plus_error.
(**
This example is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2016-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.
*)
Require Import Reals Psatz.
From Flocq Require Import Core Operations Relative Plus_error.
Section Theory.
......@@ -67,7 +80,7 @@ Lemma hombnd_plus :
forall m M u1 v1 b1 B1 u2 v2 b2 B2,
hombnd m M u1 v1 b1 B1 ->
hombnd m M u2 v2 b2 B2 ->
hombnd m M (u1 + u2) (v1 + v2) (Fplus radix2 b1 b2) (Fplus radix2 B1 B2).
hombnd m M (u1 + u2) (v1 + v2) (Fplus b1 b2) (Fplus B1 B2).
Proof.
intros m M u1 v1 b1 B1 u2 v2 b2 B2 [H11 [H12 H1]] [H21 [H22 H2]].
unfold hombnd.
......@@ -91,7 +104,7 @@ Lemma hombnd_minus :
forall m M u1 v1 b1 B1 u2 v2 b2 B2,
hombnd m M u1 v1 b1 B1 ->
hombnd m M u2 v2 b2 B2 ->
hombnd m M (u1 - u2) (v1 - v2) (Fplus radix2 b1 b2) (Fplus radix2 B1 B2).
hombnd m M (u1 - u2) (v1 - v2) (Fplus b1 b2) (Fplus B1 B2).
Proof.
intros m M u1 v1 b1 B1 u2 v2 b2 B2 H1 [H21 [H22 H2]].
apply hombnd_plus with (1 := H1).
......@@ -104,13 +117,13 @@ now split.
Qed.
Definition mult_err b1 B1 b2 B2 :=
Fplus radix2 (Fplus radix2 (Fmult radix2 b1 B2) (Fmult radix2 B1 b2)) (Fmult radix2 b1 b2).
Fplus (Fplus (Fmult b1 B2) (Fmult B1 b2)) (@Fmult radix2 b1 b2).
Lemma hombnd_mult :
forall m M1 u1 v1 b1 B1 M2 u2 v2 b2 B2,
hombnd m M1 u1 v1 b1 B1 ->
hombnd m M2 u2 v2 b2 B2 ->
hombnd m (M1 * M2) (u1 * u2) (v1 * v2) (mult_err b1 B1 b2 B2) (Fmult radix2 B1 B2).
hombnd m (M1 * M2) (u1 * u2) (v1 * v2) (mult_err b1 B1 b2 B2) (Fmult B1 B2).
Proof.
intros m M1 u1 v1 b1 B1 M2 u2 v2 b2 B2 [H11 [H12 H1]] [H21 [H22 H2]].
unfold hombnd, mult_err.
......@@ -151,7 +164,7 @@ now apply Rmult_le_compat ; try apply Rabs_pos.
Qed.
Definition round_err b B :=
Fplus radix2 (Fmult radix2 (Fplus radix2 b B) (Float radix2 1 (- prec))) b.
Fplus (Fmult (Fplus b B) (Float radix2 1 (- prec))) b.
Lemma hombnd_rnd :
forall m M u v b B,
......@@ -166,7 +179,7 @@ apply Rplus_le_le_0_compat with (2 := Ho1).
apply Rmult_le_pos.
apply Rplus_le_le_0_compat with (1 := Ho1).
now apply Rle_trans with (1 := Rle_0_1).
now apply F2R_ge_0_compat.
now apply F2R_ge_0.
apply (conj Ho2).
intros H.
specialize (Ho (Rle_trans _ _ _ H (Rmin_l _ _))).
......@@ -233,7 +246,7 @@ Lemma hombnd_sub_init :
Proof.
intros u v Fu Fv.
split.
now apply F2R_ge_0_compat.
now apply F2R_ge_0.
unfold F2R at 1 3 ; simpl.
rewrite 2!Rmult_1_l.
repeat split ; try apply Rle_refl.
......@@ -243,7 +256,7 @@ rewrite round_generic.
unfold Rminus at 1.
rewrite Rplus_opp_r, Rabs_R0.
apply Rmult_le_pos.
now apply F2R_ge_0_compat.
now apply F2R_ge_0.
apply Rabs_pos.
apply valid_rnd_N.
apply FLT_format_plus_small ; try easy.
......@@ -282,7 +295,7 @@ Lemma hombnd_add :
forall {m M u1 v1 b1 B1 u2 v2 b2 B2},
hombnd' m M u1 v1 b1 B1 ->
hombnd' m M u2 v2 b2 B2 ->
hombnd' m M (u1 + u2) (v1 + v2) (Fplus radix2 b1 b2) (Fplus radix2 B1 B2).
hombnd' m M (u1 + u2) (v1 + v2) (Fplus b1 b2) (Fplus B1 B2).
Proof.
apply hombnd_plus.
Qed.
......@@ -291,7 +304,7 @@ Lemma hombnd_sub :
forall {m M u1 v1 b1 B1 u2 v2 b2 B2},
hombnd' m M u1 v1 b1 B1 ->
hombnd' m M u2 v2 b2 B2 ->
hombnd' m M (u1 - u2) (v1 - v2) (Fplus radix2 b1 b2) (Fplus radix2 B1 B2).
hombnd' m M (u1 - u2) (v1 - v2) (Fplus b1 b2) (Fplus B1 B2).
Proof.
apply hombnd_minus.
Qed.
......@@ -300,7 +313,7 @@ Lemma hombnd_mul :
forall {m M1 u1 v1 b1 B1 M2 u2 v2 b2 B2},
hombnd' m M1 u1 v1 b1 B1 ->
hombnd' m M2 u2 v2 b2 B2 ->
hombnd' m (M1 * M2) (u1 * u2) (v1 * v2) (mult_err b1 B1 b2 B2) (Fmult radix2 B1 B2).
hombnd' m (M1 * M2) (u1 * u2) (v1 * v2) (mult_err b1 B1 b2 B2) (Fmult B1 B2).
Proof.
apply hombnd_mult.
Qed.
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
(**
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
#<br />#
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.
*)
(** * Necessary conditions to compute a*x+y faithfully *)
Require Import Fcore.
Section Axpy.
Variable beta : radix.
Notation bpow e := (bpow beta e).
Variable prec : Z.
Variable Hp : Zlt 0 prec.
(* FLX first - probably correct in FLTand maybe FTZ? *)
Notation format := (generic_format beta (FLX_exp prec)).
Notation cexp := (canonic_exp beta (FLX_exp prec)).
Notation ulp := (ulp beta (FLX_exp prec)).
Theorem pred_gt_0: forall f,
format f -> (0 < f)%R -> (0 < pred beta (FLX_exp prec) f)%R.
intros f Hf Zf.
unfold pred, Fcore_ulp.ulp, canonic_exp, FLX_exp.
destruct (ln_beta beta f) as (ef,Hef).
simpl.
assert (Zf2: (f <>0)%R) by now apply Rgt_not_eq.
specialize (Hef Zf2).
rewrite Rabs_right in Hef.
2: now apply Rgt_ge.
case (Zle_lt_or_eq 1 prec).
omega.
intros Hp1.
case Req_bool_spec; intros H; apply Rlt_Rminus;
apply Rlt_le_trans with (2:=proj1 Hef);
apply bpow_lt; omega.
(* special case for p = 1 *)
intros Hp1.
case Req_bool_spec; intros H; apply Rlt_Rminus.
apply Rlt_le_trans with (2:=proj1 Hef).
apply bpow_lt; omega.
rewrite <- Hp1.
case (proj1 Hef); trivial.
intros H'.
now contradict H.
Qed.
Definition MinOrMax x f :=
((f = round beta (FLX_exp prec) Zfloor x)
\/ (f = round beta (FLX_exp prec) Zceil x)).
Theorem MinOrMax_opp: forall x f,
MinOrMax x f <-> MinOrMax (-x) (-f).
assert (forall x f, MinOrMax x f -> MinOrMax (- x) (- f)).
unfold MinOrMax; intros x f [H|H].
right.
now rewrite H, round_UP_opp.
left.
now rewrite H, round_DN_opp.
intros x f; split; intros H1.
now apply H.
rewrite <- (Ropp_involutive x), <- (Ropp_involutive f).
now apply H.
Qed.
Theorem implies_DN_lt_ulp:
forall x f, format f ->
(0 < f <= x)%R ->
(Rabs (f-x) < ulp f)%R ->
(f = round beta (FLX_exp prec) Zfloor x)%R.
intros x f Hf Hxf1 Hxf2.
apply sym_eq.
replace x with (f+-(f-x))%R by ring.
apply round_DN_succ; trivial.
apply Hxf1.
replace (- (f - x))%R with (Rabs (f-x)).
split; trivial; apply Rabs_pos.
rewrite Rabs_left1; trivial.
now apply Rle_minus.
Qed.
Theorem MinOrMax_ulp_pred:
forall x f, format f ->
(0 < f)%R ->
(Rabs (f-x) < ulp (pred beta (FLX_exp prec) f))%R ->
MinOrMax x f.
intros x f Ff Zf Hf.
case (Rlt_or_le x f); intros Hxf2.
(* *)
(* assert (ulp (f-ulp f) = ulp f)%R.
admit.
right; apply sym_eq.
replace f with ((f-ulp f) + (ulp (f-ulp f)))%R.
2: rewrite H; ring.
replace x with ((f-ulp f)+-(f-ulp f-x))%R by ring.
apply round_UP_succ; trivial.
apply Hxf1.
replace (- (f - x))%R with (Rabs (f-x)).
split; trivial; apply Rabs_pos.
rewrite Rabs_left1; trivial.
now apply Rle_minus.
*)
admit.
(* *)
left.
apply implies_DN_lt_ulp; trivial.
now split.
apply Rlt_le_trans with (1:=Hf).
apply ulp_monotone.
clear; intros m n H.
unfold FLX_exp.
omega.
now apply pred_gt_0.
left.
apply pred_lt.
Qed.
Theorem implies_MinOrMax_bpow:
forall x f, format f ->
f = bpow (ln_beta beta f) ->
(Rabs (f-x) < /2* ulp f)%R ->
MinOrMax x f.
intros x f Hf1 Hf2 Hxf.
Admitted.
Variable choice : Z -> bool.
Variable a1 x1 y1 a x y:R.
Hypothesis Ha: format a.
Hypothesis Hx: format x.
Hypothesis Hy: format y.
Notation t := (round beta (FLX_exp prec) (Znearest choice) (a*x)).
Notation u := (round beta (FLX_exp prec) (Znearest choice) (t+y)).
(*
Axpy_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> abs(t-a*x) <= Fulp(b)(Fpred(b)(u))/4
=> abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
=> MinOrMax?(y1+a1*x1,u)
Axpy_aux1_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Fnormal?(b)(t) => radix*abs(t) <= Fpred(b)(u)
=> abs(t-a*x) <= Fulp(b)(Fpred(b)(u))/4
Axpy_aux1_aux2 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Fsubnormal?(b)(t) => 1-dExp(b) <= Fexp(Fpred(b)(u))
=> abs(t-a*x) <= Fulp(b)(Fpred(b)(u))/4
Axpy_aux2 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Fsubnormal?(b)(t) => u=t+y
=> abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
=> MinOrMax?(y1+a1*x1,u)
Axpy_aux3 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Fsubnormal?(b)(t)
=> -dExp(b) = Fexp(Fpred(b)(u)) => 1-dExp(b) <= Fexp(u)
=> abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
=> MinOrMax?(y1+a1*x1,u)
AxpyPos : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> (Fnormal?(b)(t) => radix*abs(t) <= Fpred(b)(u))
=> abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
=> MinOrMax?(y1+a1*x1,u)
Axpy_opt_aux1_aux1 : lemma Fnormal?(b)(t) => Fnormal?(b)(u) => 0 < u
=> Prec(b) >= 3 =>
(1+radix*(1+1/(2*abs(Fnum(u))))*(1+1/abs(Fnum(Fpred(b)(u)))))/(1-1/(2*abs(Fnum(t))))
<= 1+radix+radix^(4-Prec(b))
Axpy_opt_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Prec(b) >= 3 => Fnormal?(b)(t)
=> (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
=> radix*abs(t) <= Fpred(b)(u)
Axpy_opt_aux2 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Prec(b) >= 6 => Fnormal?(b)(t)
=> (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
=> abs(y)*radix^(1-Prec(b))/(radix+1) < Fulp(b)(Fpred(b)(u))
Axpy_opt_aux3 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Prec(b) >= 3 => Fsubnormal?(b)(t)
=> (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
=> abs(y)*radix^(1-Prec(b))/(radix+radix/2) <= Fulp(b)(Fpred(b)(u))
Axpy_optPos : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
=> Prec(b) >= 6
=> (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
=> abs(y1-y+a1*x1-a*x) < abs(y)*radix^(1-Prec(b))/(6*radix)
=> MinOrMax?(y1+a1*x1,u)