Library Interval.Interval.Float_full_primfloat
This file is part of the Coq.Interval library for proving bounds of
real-valued expressions in Coq: https://coqinterval.gitlabpages.inria.fr/
Copyright (C) 2023-2023, Inria
This library is governed by the CeCILL-C license under French law and
abiding by the rules of distribution of free software. You can use,
modify and/or redistribute the library under the terms of the CeCILL-C
license as circulated by CEA, CNRS and Inria at the following URL:
http://www.cecill.info/
As a counterpart to the access to the source code and rights to copy,
modify and redistribute granted by the license, users are provided
only with a limited warranty and the library's author, the holder of
the economic rights, and the successive licensors have only limited
liability. See the COPYING file for more details.
From Coq Require Import Reals Psatz PrimInt63 Uint63 Sint63 PArray Floats.
From Flocq Require Import Core PrimFloat BinarySingleNaN.
Require Import Missing.Flocq.
Require Import Xreal.
Require Import Basic.
Require Import Sig.
Require Import Interval.
Require Import Float.
Require Import Float_full.
Require Import Specific_bigint Specific_ops.
Require Import Primitive_ops.
Require Import Lang_expr Lang_tac.
Require Import Interval_helper.
Local Open Scope R_scope.
Lemma sub_sub_assoc : ∀ x y z, (x - (y - z) = x - y + z)%uint63.
Lemma asr_land : ∀ x y, (Uint63.to_Z y ≤ 62)%Z →
to_Z x = (to_Z (asr x y) × (2 ^ Uint63.to_Z y) + Uint63.to_Z (x land (lsl 1 y - 1)))%Z.
Lemma ulp_FLX_le_FLT : ∀ (beta : radix) (emin prec : Z), Prec_gt_0 prec →
∀ x : R, ulp beta (FLX_exp prec) x ≤ ulp beta (FLT_exp emin prec) x.
Lemma ulp_FLX_eq_FLT_large : ∀ (beta : radix) (emin prec : Z), Prec_gt_0 prec →
∀ x : R, bpow beta (emin + prec - 1) ≤ Rabs x →
ulp beta (FLX_exp prec) x = ulp beta (FLT_exp emin prec) x.
Lemma succ_FLX_le_FLT : ∀ beta emin prec, Prec_gt_0 prec → ∀ x,
succ beta (FLX_exp prec) x ≤ succ beta (FLT_exp emin prec) x.
Lemma pred_FLT_le_FLX : ∀ beta emin prec, Prec_gt_0 prec → ∀ x,
pred beta (FLT_exp emin prec) x ≤ pred beta (FLX_exp prec) x.
Lemma succ_FLX_eq_FLT_large : ∀ beta emin prec, Prec_gt_0 prec → ∀ x,
bpow beta (emin + prec - 1) ≤ x →
succ beta (FLX_exp prec) x = succ beta (FLT_exp emin prec) x.
Lemma FLX_FLT_aux : ∀ (emin prec e : Z),
(emin + prec - 1 < e)%Z → FLX_exp prec e = FLT_exp emin prec e.
Lemma pred_FLT_eq_FLX_large : ∀ beta emin prec, Prec_gt_0 prec → ∀ x,
bpow beta (emin + prec - 1) < x →
pred beta (FLT_exp emin prec) x = pred beta (FLX_exp prec) x.
Lemma succ_round_gt_id : ∀ (beta : radix) (fexp : Z → Z), Valid_exp fexp →
∀ rnd : R → Z, Valid_rnd rnd → ∀ x : R, ulp beta fexp x ≠ 0 →
x < succ beta fexp (Generic_fmt.round beta fexp rnd x).
Lemma succ_round_FLX_le_FLT : ∀ beta emin prec, Prec_gt_0 prec →
∀ rnd, Valid_rnd rnd → ∀ x, generic_format beta (FLX_exp prec) x →
succ beta (FLX_exp prec) x ≤
succ beta (FLT_exp emin prec) (Generic_fmt.round beta (FLT_exp emin prec) rnd x).
Lemma pred_round_FLT_le_FLX : ∀ beta emin prec, Prec_gt_0 prec →
∀ rnd, Valid_rnd rnd → ∀ x, generic_format beta (FLX_exp prec) x →
pred beta (FLT_exp emin prec) (Generic_fmt.round beta (FLT_exp emin prec) rnd x)
≤ pred beta (FLX_exp prec) x.
Lemma pred_FLX_exact_shift : ∀ beta prec, Prec_gt_0 prec → ∀ x e,
pred beta (FLX_exp prec) (x × bpow beta e) = pred beta (FLX_exp prec) x × bpow beta e.
Lemma succ_FLT_shift_ge : ∀ beta emin prec, Prec_gt_0 prec →
∀ rnd, Valid_rnd rnd → ∀ x, generic_format beta (FLT_exp emin prec) x →
bpow beta (emin + prec - 1) ≤ x → ∀ e,
succ beta (FLT_exp emin prec) x × bpow beta e ≤
succ beta (FLT_exp emin prec) (Generic_fmt.round beta (FLT_exp emin prec) rnd (x × bpow beta e)).
Lemma pred_FLT_shift_le : ∀ beta emin prec, Prec_gt_0 prec →
∀ rnd, Valid_rnd rnd → ∀ x, generic_format beta (FLT_exp emin prec) x →
bpow beta (emin + prec - 1) < x → ∀ e,
pred beta (FLT_exp emin prec) (Generic_fmt.round beta (FLT_exp emin prec) rnd (x × bpow beta e))
≤ pred beta (FLT_exp emin prec) x × bpow beta e.
Module PrimFloatIntervalFull <: IntervalOps.
Module Faux := SpecificFloat BigIntRadix2.
Module Iaux := FloatIntervalFull Faux.
Module IT := IntervalTacticAux Iaux.
Import IT.SimpleTactic.
Module I := FloatIntervalFull PrimitiveFloat.
Import I.
Definition pi (prec : F.precision) : type :=
(Ibnd 0x1.921fb54442d18p+1 0x1.921fb54442d19p+1)%float.
Theorem pi_correct : ∀ prec, contains (convert (pi prec)) (Xreal PI).
Include FloatInterval PrimitiveFloat.
Definition cos := cos.
Definition sin := sin.
Definition tan := tan.
Definition atan := atan.
Definition ln := ln.
Definition cos_correct := cos_correct.
Definition sin_correct := sin_correct.
Definition tan_correct := tan_correct.
Definition atan_correct := atan_correct.
Definition ln_correct := ln_correct.
Module ExpImpl.
Definition consts := [|
0x1.0000000000000p0%float;
0x1.02c9a3e778061p0%float;
0x1.059b0d3158574p0%float;
0x1.0874518759bc8p0%float;
0x1.0b5586cf9890fp0%float;
0x1.0e3ec32d3d1a2p0%float;
0x1.11301d0125b51p0%float;
0x1.1429aaea92de0p0%float;
0x1.172b83c7d517bp0%float;
0x1.1a35beb6fcb75p0%float;
0x1.1d4873168b9aap0%float;
0x1.2063b88628cd6p0%float;
0x1.2387a6e756238p0%float;
0x1.26b4565e27cddp0%float;
0x1.29e9df51fdee1p0%float;
0x1.2d285a6e4030bp0%float;
0x1.306fe0a31b715p0%float;
0x1.33c08b26416ffp0%float;
0x1.371a7373aa9cbp0%float;
0x1.3a7db34e59ff7p0%float;
0x1.3dea64c123422p0%float;
0x1.4160a21f72e2ap0%float;
0x1.44e086061892dp0%float;
0x1.486a2b5c13cd0p0%float;
0x1.4bfdad5362a27p0%float;
0x1.4f9b2769d2ca7p0%float;
0x1.5342b569d4f82p0%float;
0x1.56f4736b527dap0%float;
0x1.5ab07dd485429p0%float;
0x1.5e76f15ad2148p0%float;
0x1.6247eb03a5585p0%float;
0x1.6623882552225p0%float;
0x1.6a09e667f3bcdp0%float;
0x1.6dfb23c651a2fp0%float;
0x1.71f75e8ec5f74p0%float;
0x1.75feb564267c9p0%float;
0x1.7a11473eb0187p0%float;
0x1.7e2f336cf4e62p0%float;
0x1.82589994cce13p0%float;
0x1.868d99b4492edp0%float;
0x1.8ace5422aa0dbp0%float;
0x1.8f1ae99157736p0%float;
0x1.93737b0cdc5e5p0%float;
0x1.97d829fde4e50p0%float;
0x1.9c49182a3f090p0%float;
0x1.a0c667b5de565p0%float;
0x1.a5503b23e255dp0%float;
0x1.a9e6b5579fdbfp0%float;
0x1.ae89f995ad3adp0%float;
0x1.b33a2b84f15fbp0%float;
0x1.b7f76f2fb5e47p0%float;
0x1.bcc1e904bc1d2p0%float;
0x1.c199bdd85529cp0%float;
0x1.c67f12e57d14bp0%float;
0x1.cb720dcef9069p0%float;
0x1.d072d4a07897cp0%float;
0x1.d5818dcfba487p0%float;
0x1.da9e603db3285p0%float;
0x1.dfc97337b9b5fp0%float;
0x1.e502ee78b3ff6p0%float;
0x1.ea4afa2a490dap0%float;
0x1.efa1bee615a27p0%float;
0x1.f50765b6e4540p0%float;
0x1.fa7c1819e90d8p0%float|
0%float|].
Definition InvLog2_64 := 0x1.71547652b82fep6%float.
Definition Log2div64h := 0x1.62e42fefap-7%float.
Definition Log2div64l := 0x1.cf79abc9e3b3ap-46%float.
Definition q1 := 0x1p0%float.
Definition q2 := 0x1.fffffffffdb3bp-2%float.
Definition q3 := 0x1.555555555653ep-3%float.
Definition q4 := 0x1.555573f218b93p-5%float.
Definition q5 := 0x1.111112d9f54c8p-7%float.
Definition g0 : ArithExpr (BinFloat :: nil) BinFloat :=
Op MUL (Var 0) (Op ADD (BinFl q1) (Op MUL (Var 0) (Op ADD (BinFl q2) (Op MUL (Var 0) (Op ADD (BinFl q3) (Op MUL (Var 0) (Op ADD (BinFl q4) (Op MUL (Var 0) (BinFl q5))))))))).
Definition Papprox (t : PrimFloat.float) := Eval cbv in evalPrim g0 (t, tt).
Definition exp_aux (x : F.type) :=
if PrimFloat.ltb x (-0x1.74385446d71c4p9)%float then (0%float, 0x1p-1074%float) else
if PrimFloat.ltb 0x1.62e42fefa39efp9%float x then (0x1.fffffffffff2ap1023%float, infinity) else
let k0 := (x × InvLog2_64 + 0x1.8p52)%float in let kf := (k0 - 0x1.8p52)%float in
let tf := (x - kf × Log2div64h - kf × Log2div64l)%float in
let ki := (normfr_mantissa (fst (frshiftexp k0)) - 6755399440921280)%uint63 in
let C := consts.[PrimInt63.land ki 63] in
let kq := PrimInt63.asr ki 6 in let y := (C × Papprox tf)%float in
let lb := (C + (y + -0x1.25p-57))%float in let ub := (C + (y + 0x1.25p-57))%float in
(next_down (ldshiftexp lb kq), next_up (ldshiftexp ub kq)).
Lemma exp_aux_correct :
∀ x, is_finite (Prim2B x) = true →
(let lb := fst (exp_aux x) in
F.valid_lb lb = true ∧
match F.toX lb with
| Xreal.Xnan ⇒ True
| Xreal.Xreal r ⇒ r ≤ Rtrigo_def.exp (B2R (Prim2B x))
end) ∧
(let ub := snd (exp_aux x) in
F.valid_ub ub = true ∧
match F.toX ub with
| Xreal.Xnan ⇒ True
| Xreal.Xreal r ⇒ Rtrigo_def.exp (B2R (Prim2B x)) ≤ r
end).
End ExpImpl.
Import ExpImpl.
Definition exp (prec : F.precision) xi :=
let aux x :=
let k0 := (x × InvLog2_64 + 0x1.8p52)%float in let kf := (k0 - 0x1.8p52)%float in
let tf := (x - kf × Log2div64h - kf × Log2div64l)%float in
let ki := (normfr_mantissa (fst (frshiftexp k0)) - 6755399440921280)%uint63 in
let C := consts.[PrimInt63.land ki 63] in
let kq := PrimInt63.asr ki 6 in let y := (C × Papprox tf)%float in
(C, y, kq) in
match xi with
| Ibnd xl xu ⇒
Ibnd
(if F.real xl then
if PrimFloat.ltb xl (-0x1.74385446d71c4p9)%float then 0%float else
if PrimFloat.ltb 0x1.62e42fefa39efp9%float xl then 0x1.fffffffffff2ap1023%float else
let '(C, y, kq) := aux xl in
next_down (ldshiftexp (C + (y + -0x1.25p-57))%float kq)
else 0%float)
(if F.real xu then
if PrimFloat.ltb xu (-0x1.74385446d71c4p9)%float then 0x1p-1074%float else
if PrimFloat.ltb 0x1.62e42fefa39efp9%float xu then infinity else
let '(C, y, kq) := aux xu in
next_up (ldshiftexp (C + (y + 0x1.25p-57))%float kq)
else nan)
| Inan ⇒ Inan
end.
Theorem exp_correct :
∀ prec, extension Xexp (exp prec).
End PrimFloatIntervalFull.