Library Flocq.Pff.Pff2Flocq

This file is part of the Flocq formalization of floating-point arithmetic in Coq: https://flocq.gitlabpages.inria.fr/
Copyright (C) 2016-2018 Sylvie Boldo
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 Psatz.

Require Import Core Plus_error Mult_error Operations Sterbenz.
Require Import Pff Pff2FlocqAux.
Require Import Nat2Z_compat.
Open Scope R_scope.

Section FTS.
Variable emin prec : Z.
Variable choice : Z → bool.
Hypothesis precisionNotZero : (1 < prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.
Hypothesis choice_sym: ∀ x, choice x = negb (choice (- (x + 1))).

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) (Znearest choice)).
Notation bpow e := (bpow radix2 e).

Lemma round_N_opp_sym: ∀ x, round_flt (- x) =
       - round_flt x.

inputs
Variable x y:R.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

algorithm
Let a := round_flt (x+y).
Let b := round_flt (y+round_flt (x-a)).

Theorem
Theorem Fast2Sum_correct: Rabs y ≤ Rabs x → a+b=x+y.

End FTS.

Section TS.

Variable emin prec : Z.
Variable choice : Z → bool.
Hypothesis precisionNotZero : (1 < prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.
Hypothesis choice_sym: ∀ x, choice x = negb (choice (- (x + 1))).

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) (Znearest choice)).
Notation bpow e := (bpow radix2 e).

inputs
Variable x y:R.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

algorithm
Let a := round_flt (x+y).
Let x' := round_flt (a-x).
Let dx := round_flt (x- round_flt (a-x')).
Let dy := round_flt (y - x').
Let b := round_flt (dx + dy).

Theorem
Theorem TwoSum_correct: a+b=x+y.

End TS.

Section Veltkamp.

Variable beta : radix.
Variable emin prec : Z.
Variable choice : Z → bool.
Variable s:Z.
Hypothesis precisionGe3 : (3 ≤ prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format beta (FLT_exp emin prec)).
Notation round_flt :=(round beta (FLT_exp emin prec) (Znearest choice)).
Notation round_flt_s :=(round beta (FLT_exp emin (prec-s)) (Znearest choice)).
Notation ulp_flt :=(ulp beta (FLT_exp emin prec)).
Notation bpow e := (bpow beta e).

inputs
Hypothesis SLe: (2 ≤ s)%Z.
Hypothesis SGe: (s ≤ prec-2)%Z.
Variable x:R.
Hypothesis Fx: format x.

algorithm
Let p := round_flt (x*(bpow s+1)).
Let q:= round_flt (x-p).
Let hx:=round_flt (q+p).
Let tx:=round_flt (x-hx).

Theorems

Lemma C_format: format (bpow s +1).

Theorem Veltkamp_Even:
  (choice = fun z ⇒ negb (Z.even z)) →
   hx = round_flt_s x.

Theorem Veltkamp: ∃ choice': Z→bool ,
   hx = round beta (FLT_exp emin (prec-s)) (Znearest choice') x.

Theorem Veltkamp_tail:
 x = hx+tx ∧ generic_format beta (FLT_exp emin s) tx.

End Veltkamp.

Section Underf_mult_aux.

Variable beta : radix.
Variable b: Fbound.
Variable prec: Z.

Hypothesis pGivesBound : Zpos (vNum b) = Zpower_nat beta (Z.abs_nat prec).
Hypothesis precisionGt1 : (1 < prec)%Z.

Lemma underf_mult_aux:
 ∀ e x y,
  Fbounded b x →
  Fbounded b y →
   (bpow beta (e + 2 × prec - 1)%Z ≤ Rabs (FtoR beta x × FtoR beta y)) →
     (e ≤ Pff.Fexp x + Pff.Fexp y)%Z.

Lemma underf_mult_aux':
 ∀ x y,
  Fbounded b x →
  Fbounded b y →
   (bpow beta (-dExp b + 2 × prec - 1)%Z ≤ Rabs (FtoR beta x × FtoR beta y)) →
     (-dExp b ≤ Pff.Fexp x + Pff.Fexp y)%Z.

End Underf_mult_aux.

Section Dekker.

Variable beta : radix.
Variable emin prec: Z.
Variable choice : Z → bool.
Let s:= (prec- Z.div2 prec)%Z.

Hypothesis precisionGe4 : (4 ≤ prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin < 0)%Z.

Notation format := (generic_format beta (FLT_exp emin prec)).
Notation round_flt :=(round beta (FLT_exp emin prec) (Znearest choice)).
Notation round_flt_s :=(round beta (FLT_exp emin (prec-s)) (Znearest choice)).
Notation ulp_flt :=(ulp beta (FLT_exp emin prec)).
Notation bpow e := (bpow beta e).

inputs
Variable x y:R.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

first Veltkamp
Let px := round_flt (x*(bpow s+1)).
Let qx:= round_flt (x-px).
Let hx:=round_flt (qx+px).
Let tx:=round_flt (x-hx).

second Veltkamp
Let py := round_flt (y*(bpow s+1)).
Let qy:= round_flt (y-py).
Let hy:=round_flt (qy+py).
Let ty:=round_flt (y-hy).

all products
Let x1y1:=round_flt (hx×hy).
Let x1y2:=round_flt (hx×ty).
Let x2y1:=round_flt (tx×hy).
Let x2y2:=round_flt (tx×ty).

final summation
Let r :=round_flt(x×y).
Let t1 :=round_flt (-r+x1y1).
Let t2 :=round_flt (t1+x1y2).
Let t3 :=round_flt (t2+x2y1).
Let t4 :=round_flt (t3+x2y2).

Theorem Dekker: (radix_val beta=2)%Z ∨ (Z.Even prec) →
  (x×y=0 ∨ bpow (emin + 2 × prec - 1) ≤ Rabs (x × y) → (x×y=r+t4)%R) ∧
    (Rabs (x×y-(r+t4)) ≤ (7/2)*bpow emin)%R.

End Dekker.

Section ErrFMA_V1.

Variable beta : radix.
Variable emin prec : Z.
Hypothesis Even_radix: Even beta.
Hypothesis precisionGe3 : (3 ≤ prec)%Z.
Variable choice : Z → bool.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format beta (FLT_exp emin prec)).
Notation round_flt :=(round beta (FLT_exp emin prec) (Znearest choice)).
Notation ulp_flt :=(ulp beta (FLT_exp emin prec)).

inputs
Variable a x y:R.
Hypothesis Fa: format a.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

algorithm
Let r1 := round_flt (a×x+y).
Let u1 := round_flt (a×x).
Let u2 := a×x-u1.
Let alpha1 := round_flt (y+u2).
Let alpha2 := (y+u2)-alpha1.
Let beta1 := round_flt (u1+alpha1).
Let beta2 := (u1+alpha1) - beta1.
Let gamma := round_flt (round_flt (beta1-r1)+beta2).
Let r2 := round_flt (gamma+alpha2).
Let r3 := (gamma+alpha2) -r2.

Non-underflow hypotheses
Hypothesis V1_Und1: a × x = 0 ∨ bpow beta (emin + 2 × prec - 1) ≤ Rabs (a × x).
Hypothesis V1_Und2: alpha1 = 0 ∨ bpow beta (emin + prec) ≤ Rabs alpha1.
Hypothesis V1_Und4: beta1 = 0 ∨ bpow beta (emin + prec+1) ≤ Rabs beta1.
Hypothesis V1_Und5: r1 = 0 ∨ bpow beta (emin + prec-1) ≤ Rabs r1.

Deduced from non-underflow hypotheses
Lemma V1_Und3': u1 = 0 ∨ bpow beta (emin + 2×prec-1) ≤ Rabs u1.

Lemma V1_Und3: u1 = 0 ∨ bpow beta (emin + prec) ≤ Rabs u1.

Theorems
Lemma ErrFMA_bounded: format r1 ∧ format r2 ∧ format r3.

Lemma ErrFMA_correct:
   a×x+y = r1+r2+r3.

End ErrFMA_V1.

Section ErrFMA_V2.

Variable beta : radix.
Variable emin prec : Z.
Hypothesis Even_radix: Even beta.
Hypothesis precisionGe3 : (3 ≤ prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format beta (FLT_exp emin prec)).
Notation round_flt :=(round beta (FLT_exp emin prec) ZnearestE).
Notation ulp_flt :=(ulp beta (FLT_exp emin prec)).

Lemma mult_error_FLT_ge_bpow': ∀ a b e, format a → format b →
   a×b = 0 ∨ bpow beta e ≤ Rabs (a×b) →
   a×b-round_flt (a×b) = 0 ∨
     bpow beta (e+1-2×prec) ≤ Rabs (a×b-round_flt (a×b)).

inputs
Variable a x y:R.
Hypothesis Fa: format a.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

algorithm
Let r1 := round_flt (a×x+y).
Let u1 := round_flt (a×x).
Let u2 := a×x-u1.
Let alpha1 := round_flt (y+u2).
Let alpha2 := (y+u2)-alpha1.
Let beta1 := round_flt (u1+alpha1).
Let beta2 := (u1+alpha1) - beta1.
Let gamma := round_flt (round_flt (beta1-r1)+beta2).
Let r2 := round_flt (gamma+alpha2).
Let r3 := (gamma+alpha2) -r2.

Non-underflow hypotheses
Hypothesis U1: a × x = 0 ∨ bpow beta (emin + 4×prec - 3) ≤ Rabs (a × x).
Hypothesis U2: y = 0 ∨ bpow beta (emin + 2×prec) ≤ Rabs y.

Lemma ErrFMA_bounded_simpl: format r1 ∧ format r2 ∧ format r3.

Lemma V2_Und4: a×x ≠ 0 → beta1 = 0 ∨ bpow beta (emin + prec+1) ≤ Rabs beta1.

Lemma V2_Und2: y ≠ 0 → alpha1 = 0 ∨ bpow beta (emin + prec) ≤ Rabs alpha1.

Lemma V2_Und5: a×x ≠ 0 → r1 = 0 ∨ bpow beta (emin + prec-1) ≤ Rabs r1.

Lemma ErrFMA_correct_simpl:
   a×x+y = r1+r2+r3.

End ErrFMA_V2.

Section ErrFmaApprox.

Variable beta : radix.
Variable emin prec : Z.
Hypothesis precisionGe3 : (4 ≤ prec)%Z.
Variable choice : Z → bool.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format beta (FLT_exp emin prec)).
Notation round_flt :=(round beta (FLT_exp emin prec) (Znearest choice)).
Notation ulp_flt :=(ulp beta (FLT_exp emin prec)).

inputs
Variable a x y:R.
Hypothesis Fa: format a.
Hypothesis Fx: format x.
Hypothesis Fy: format y.

algorithm
Let r1 := round_flt (a×x+y).
Let u1 := round_flt (a×x).
Let u2 := a×x-u1.
Let v1 := round_flt (y+u1).
Let v2 := (y+u1)-v1.
Let t1 := round_flt (v1-r1).
Let t2 := round_flt (u2+v2).
Let r2 := round_flt (t1+t2).

Non-underflow hypotheses
Hypothesis Und1: a × x = 0 ∨ bpow beta (emin + 2 × prec - 1) ≤ Rabs (a × x).
Hypothesis Und2: v1 = 0 ∨ bpow beta (emin + prec - 1) ≤ Rabs v1.
Hypothesis Und3: r1 = 0 ∨ bpow beta (emin + prec - 1) ≤ Rabs r1.
Hypothesis Und4: r2 = 0 ∨ bpow beta (emin + prec - 1) ≤ Rabs r2.
Hypothesis Und5: t2 = 0 ∨ bpow beta (emin + prec - 1) ≤ Rabs t2.


Lemma ErrFmaAppr_correct:
   Rabs (r1+r2 -(a×x+y)) ≤ (3×beta/2+/2) × bpow beta (2-2×prec)%Z × Rabs (r1).

End ErrFmaApprox.

Section Axpy.
Variable emin prec : Z.
Variable choice : Z → bool.
Hypothesis precisionGt : (1 < prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) (Znearest choice)).
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
Notation bpow e := (bpow radix2 e).

inputs
Variable a x y:R.
Variable ta tx ty:R.
Hypothesis Fta: format ta.
Hypothesis Ftx: format tx.
Hypothesis Fty: format ty.

algorithm
Let tv := round_flt (ty+ round_flt (ta×tx)).

Hypotheses
Hypothesis H1 : (5+4×bpow (-prec))/(1-bpow (-prec))×
   (Rabs (ta×tx)+ bpow (emin-1)) ≤ Rabs ty.

Hypothesis H2 : Rabs (y-ty) + Rabs (a×x-ta×tx) ≤
    bpow (-prec-2)*(1-bpow (1-prec))*Rabs ty -
    bpow (-prec-2)×Rabs (ta×tx) - bpow (emin-2).

Theorem Axpy :
  tv = round radix2 (FLT_exp emin prec) Zfloor (y+a×x) ∨
     tv = round radix2 (FLT_exp emin prec) Zceil (y+a×x).

End Axpy.

Section Discri1.
Variable emin prec : Z.
Hypothesis precisionGt : (1 < prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) ZnearestE).
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
Notation bpow e := (bpow radix2 e).

inputs
Variable a b c:R.
Hypothesis Fa: format a.
Hypothesis Fb: format b.
Hypothesis Fc: format c.

algorithm
Let p:=round_flt (b×b).
Let q:=round_flt (a×c).
Let dp:= b×b-p.
Let dq:= a×c-q.
Let d:= if (Rle_bool (p+q) (3×Rabs (p-q)))
            then round_flt (p-q)
            else round_flt (round_flt (p-q) + round_flt(dp-dq)).

No-underflow hypotheses
Hypothesis U1 : b × b ≠ 0 →
    bpow (emin + 3 × prec) ≤ Rabs (b × b).
Hypothesis U2 : a × c ≠ 0 →
    bpow (emin + 3 × prec) ≤ Rabs (a × c).
Hypothesis Zd : d ≠ 0.

Lemma U3_discri1 : b × b ≠ 0 → a × c ≠ 0 → p - q ≠ 0 →
   bpow (emin + 2×prec) ≤ Rabs (round_flt (p-q)).

Lemma U4_discri1 : b × b ≠ 0 → a × c ≠ 0 → p - q ≠ 0 →
    bpow (emin + prec) ≤ Rabs d.

Lemma format_dp : format dp.

Lemma format_dq : format dq.

Lemma format_d_discri1 : format d.

Lemma U5_discri1_aux : ∀ x y e, format x → format y
   → (emin ≤ e)%Z → bpow e ≤ Rabs x → bpow e ≤ Rabs y
   → round_flt (x+y) ≠ x+y
   → bpow e ≤ Rabs (round_flt (x+y)).

Lemma U5_discri1 : b × b ≠ 0 → a×c ≠ 0 →
  round_flt (dp - dq) ≠ dp - dq →
  bpow (emin + prec - 1) ≤ Rabs (round_flt (dp - dq)).

Correctness of d
Theorem discri_correct_test :
 Rabs (d-(b×b-a×c)) ≤ 2× ulp_flt d.

End Discri1.

Section Discri2.
Variable emin prec : Z.
Hypothesis precisionGt : (4 ≤ prec)%Z.
Context { prec_gt_0_ : Prec_gt_0 prec }.
Hypothesis emin_neg: (emin ≤ 0)%Z.

Notation format := (generic_format radix2 (FLT_exp emin prec)).
Notation round_flt :=(round radix2 (FLT_exp emin prec) ZnearestE).
Notation ulp_flt :=(ulp radix2 (FLT_exp emin prec)).
Notation bpow e := (bpow radix2 e).

inputs
Variable a b c:R.
Hypothesis Fa: format a.
Hypothesis Fb: format b.
Hypothesis Fc: format c.

algorithm
Let p:=round_flt (b×b).
Let q:=round_flt (a×c).
Let dp:= b×b-p.
Let dq:= a×c-q.
Let d:= if (Rle_bool (round_flt (p+q)) (round_flt (3×Rabs (round_flt (p-q)))))
            then round_flt (p-q)
            else round_flt (round_flt (p-q) + round_flt(dp-dq)).

No-underflow hypotheses
Hypothesis U1 : b × b ≠ 0 →
    bpow (emin + 3 × prec) ≤ Rabs (b × b).
Hypothesis U2 : a × c ≠ 0 →
    bpow (emin + 3 × prec) ≤ Rabs (a × c).
Hypothesis Zd : d ≠ 0.

Lemma format_d_discri2 : format d.

Correctness of d
Theorem discri_fp_test :
 Rabs (d-(b×b-a×c)) ≤ 2× ulp_flt d.

End Discri2.