Library Interval.Interval.Transcend

This file is part of the Coq.Interval library for proving bounds of real-valued expressions in Coq: https://coqinterval.gitlabpages.inria.fr/

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.
From Flocq Require Import Zaux Raux Digits.

Require Import Stdlib.
Require Import Xreal.
Require Import Basic.
Require Import Sig.
Require Import Interval.
Require Import Float.

Module TranscendentalFloatFast (F : FloatOps with Definition sensible_format := true).

Module I := FloatInterval F.
Module J := IntervalBasicExt I.
Module F' := FloatExt F.

CoInductive hidden_constant : Type :=
  | HConst : I.type → hidden_constant.

CoInductive constants : Type :=
  | Consts: hidden_constant → constants → constants.

Fixpoint constant_getter_aux n cst :=
  match n, cst with
  | O, Consts (HConst xi) _ ⇒ xi
  | S p, Consts _ c ⇒ constant_getter_aux p c
  end.

Definition constant_getter cst prec :=
  let nb := Z.abs_nat (Z.pred (fst (Zdiv.Zdiv_eucl_POS (F.prec prec + 30)%positive 31%Z))) in
  constant_getter_aux nb cst.

CoFixpoint hidden_constant_generator gen n :=
  HConst (gen (F.PtoP (n × 31)%positive)).

CoFixpoint constant_generator_aux gen n :=
  Consts (hidden_constant_generator gen n) (constant_generator_aux gen (Pos.succ n)).

Definition constant_generator gen :=
  constant_generator_aux gen 1.

Definition Z2nat x :=
  match x with
  | Zpos p ⇒ nat_of_P p
  | _ ⇒ O
  end.

Definition Z2P x :=
  match x with
  | Zpos p ⇒ p
  | _ ⇒ xH
  end.

Definition c1 := F.fromZ 1.
Definition c2 := F.fromZ 2.
Definition c3 := F.fromZ 3.
Definition cm1 := F.fromZ (-1).
Definition i1 := I.fromZ_small 1.
Definition i2 := I.fromZ_small 2.
Definition i3 := I.fromZ_small 3.
Definition i4 := I.fromZ_small 4.
Definition i5 := I.fromZ_small 5.
Definition i6 := I.fromZ_small 6.
Definition i239 := I.fromZ_small 239.
Definition c1_2 := F.div2 c1.
Definition c1_8 := iter_pos F.div2 8 c1.
Definition c1_p_c1_8 := F.add_DN (F.PtoP 52) c1 c1_8.

Lemma c1_n_correct :
  ∀ x n,
  (/ 256 × 2^Pos.to_nat n ≤ Rabs (F.toR x))%R →
  F.toX (iter_pos F.div2 n x) = Xmul (F.toX x) (Xreal ((/2)^Pos.to_nat n)).

Lemma c1_2_correct : F.toX c1_2 = Xreal (/ 2).

Lemma c1_8_correct : F.toX c1_8 = Xreal (/ 256).

Ltac bound_tac :=
  unfold Xround, Xbind ;
  match goal with
  | |- (round ?r_DN ?p ?v ≤ ?v)%R ⇒
    apply (proj1 (proj2 (Generic_fmt.round_DN_pt F.radix (FLX.FLX_exp (Zpos p)) v)))
  | |- (?v ≤ round ?r rnd_UP ?p ?v)%R ⇒
    apply (proj1 (proj2 (Generic_fmt.round_UP_pt F.radix (FLX.FLX_exp (Zpos p)) v)))
  | |- (round ?r_DN ?p ?v ≤ ?w)%R ⇒
    apply Rle_trans with (1 := proj1 (proj2 (Generic_fmt.round_DN_pt F.radix (FLX.FLX_exp (Zpos p)) v)))
  | |- (?w ≤ round ?r rnd_UP ?p ?v)%R ⇒
    apply Rle_trans with (2 := proj1 (proj2 (Generic_fmt.round_UP_pt F.radix (FLX.FLX_exp (Zpos p)) v)))
  end.

Notation toR := F.toR (only parsing).

Fixpoint atan_fast0_aux prec thre powi sqri divi (nb : nat) { struct nb } :=
  let npwi := I.mul prec powi sqri in
  let vali := I.div prec npwi divi in
  match F.cmp (I.upper vali) thre, nb with
  | Xlt, _
  | _, O ⇒ I.bnd F.zero (I.upper vali)
  | _, S n ⇒
    I.sub prec vali
      (atan_fast0_aux prec thre npwi sqri (I.add prec divi i2) n)
  end.

Definition atan_fast0 prec xi :=
  let x2i := I.sqr prec xi in
  let p := F.prec prec in
  let thre := F.scale c1 (F.ZtoS (Zneg p)) in
  let rem := atan_fast0_aux prec thre i1 x2i i3 (nat_of_P p) in
  I.mul prec (I.sub prec i1 rem) xi.

Definition pi4_gen prec :=
  I.sub prec
   (I.mul2 prec (I.mul2 prec (atan_fast0 prec (I.inv prec i5))))
   (atan_fast0 prec (I.inv prec i239)).

Definition pi4_seq := constant_generator pi4_gen.
Definition pi4 := constant_getter pi4_seq.

Definition atan_fastP prec x :=
  let xi := I.bnd x x in
  if F'.lt c1_2 x then
    let prec := F.incr_prec prec 2 in
    let pi4i := pi4 prec in
    if F'.lt c2 x then
      I.sub prec
       (I.mul2 prec pi4i)
       (atan_fast0 prec (I.inv prec xi))
    else
      let xm1i := I.sub prec xi i1 in
      let xp1i := I.add prec xi i1 in
      I.add prec pi4i
       (atan_fast0 prec (I.div prec xm1i xp1i))
  else
    atan_fast0 (F.incr_prec prec 1) xi.

Definition atan_fast prec x :=
  match F'.cmp x F.zero with
  | Xeq ⇒ I.zero
  | Xlt ⇒ I.neg (atan_fastP prec (F.neg x))
  | Xgt ⇒ atan_fastP prec x
  | Xund ⇒ I.nai
  end.

Lemma atan_fast0_correct :
  ∀ prec xi x,
  (Rabs x ≤ /2)%R →
  contains (I.convert xi) (Xreal x) →
  contains (I.convert (atan_fast0 prec xi)) (Xreal (atan x)).

Lemma pi4_correct :
  ∀ prec, contains (I.convert (pi4 prec)) (Xreal (PI /4)).

Lemma atan_fastP_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x)%R →
  contains (I.convert (atan_fastP prec x)) (Xreal (atan (toR x))).

Lemma atan_fast_correct :
  ∀ prec x,
  contains (I.convert (atan_fast prec x)) (Xatan (F.toX x)).

Fixpoint ln1p_fast0_aux prec thre powi xi divi (nb : nat) { struct nb } :=
  let npwi := I.mul prec powi xi in
  let vali := I.div prec npwi divi in
  match F'.cmp (I.upper vali) thre, nb with
  | Xlt, _
  | _, O ⇒ I.bnd F.zero (I.upper vali)
  | _, S n ⇒
    I.sub prec vali (ln1p_fast0_aux prec thre npwi xi (I.add prec divi i1) n)
  end.

Definition ln1p_fast0 prec xi :=
  let p := F.prec prec in
  let thre := F.scale c1 (F.ZtoS (Zneg p)) in
  let rem := ln1p_fast0_aux prec thre i1 xi i2 (nat_of_P p) in
  I.mul prec (I.sub prec i1 rem) xi.

Definition ln_fast1P prec xi :=
  let th := c1_p_c1_8 in
  match F'.le' (I.upper xi) th with
  | true ⇒
    ln1p_fast0 prec (I.sub prec xi i1)
  | false ⇒
    let m := Digits.Zdigits2 (F.StoZ (F.mag (I.upper xi))) in
    let prec := F.incr_prec prec 10 in
    let fix reduce xi (nb : nat) {struct nb} :=
      match F'.le' (I.upper xi) th, nb with
      | true, _ ⇒ ln1p_fast0 prec (I.sub prec xi i1)
      | _, O ⇒ I.bnd F.zero F.nan
      | _, S n ⇒ I.mul2 prec (reduce (I.sqrt prec xi) n)
      end in
    reduce xi (8 + Z2nat m)
  end.

Definition ln_fast prec x :=
  match F.cmp x F.zero with
  | Xgt ⇒
    let xi := I.bnd x x in
    match F.cmp x c1 with
    | Xeq ⇒ I.zero
    | Xlt ⇒
      let m := Z.opp (F.StoZ (F.mag (F.sub_UP prec c1 x))) in
      let prec := F.incr_prec prec (Z2P m) in
      I.neg (ln_fast1P prec (I.inv prec xi))
    | Xgt ⇒ if F.real x then ln_fast1P prec xi else I.nai
    | Xund ⇒ I.nai
    end
  | _ ⇒ I.nai
  end.

Lemma ln1p_fast0_correct :
  ∀ prec xi x,
  (0 ≤ x ≤ /2)%R →
  contains (I.convert xi) (Xreal x) →
  contains (I.convert (ln1p_fast0 prec xi)) (Xreal (ln (1 + x))).

Lemma ln_fast1P_correct :
  ∀ prec xi x,
  (1 ≤ x)%R →
  contains (I.convert xi) (Xreal x) →
  contains (I.convert (ln_fast1P prec xi)) (Xreal (ln x)).

Theorem ln_fast_correct :
  ∀ prec x,
  contains (I.convert (ln_fast prec x)) (Xln (F.toX x)).

Fixpoint cos_fast0_aux prec thre powi sqri facti divi (nb : nat) { struct nb } :=
  let npwi := I.mul prec powi sqri in
  let vali := I.div prec npwi divi in
  match F'.cmp (I.upper vali) thre, nb with
  | Xlt, _
  | _, O ⇒ I.bnd F.zero (I.upper vali)
  | _, S n ⇒
    let nfacti := I.add prec facti i2 in
    let ndivi := I.mul prec divi (I.mul prec facti (I.add prec facti i1)) in
    I.sub prec vali (cos_fast0_aux prec thre npwi sqri nfacti ndivi n)
  end.

Definition cos_fast0 prec x :=
  let xi := I.bnd x x in
  let x2i := I.sqr prec xi in
  let p := F.prec prec in
  let thre := F.scale c1 (F.ZtoS (Zneg p)) in
  let rem := cos_fast0_aux prec thre i1 x2i i3 i2 (nat_of_P p) in
  I.sub prec i1 rem.

Definition sin_cos_reduce prec x :=
  let th := c1_2 in
  let fix reduce x (nb : nat) {struct nb} :=
    match F'.le x th, nb with
    | true, _ ⇒ (Gt , cos_fast0 prec x )
    | _, O ⇒ (Eq , I.bnd cm1 c1 )
    | _, S n ⇒
      match reduce (F.div2 x) n with
      | (s, c) ⇒
       (match s, I.sign_large c with
        | Lt, Xgt ⇒ Lt
        | Lt, Xlt ⇒ Gt
        | Lt, _ ⇒ Eq
        | Gt, Xlt ⇒ Lt
        | Gt, Xgt ⇒ Gt
        | Gt, _ ⇒ Eq
        | _, _ ⇒ s
        end,
        I.sub prec (I.mul2 prec (I.sqr prec c)) i1 )
      end
    end in
  reduce x.

Lemma cos_fast0_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (Rabs (toR x) ≤ /2)%R →
  contains (I.convert (cos_fast0 prec x)) (Xreal (cos (toR x))).

Lemma sin_cos_reduce_correct :
  ∀ prec nb x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x)%R →
  match sin_cos_reduce prec x nb with
  | (ss, ci) ⇒
    contains (I.convert ci) (Xreal (cos (toR x))) ∧
    match ss with
    | Lt ⇒ (sin (toR x) ≤ 0)%R
    | Gt ⇒ (0 ≤ sin (toR x))%R
    | _ ⇒ True
    end
  end.

Definition cos_fastP prec x :=
  let th := c1_2 in
  match F'.le' x th with
  | true ⇒ cos_fast0 prec x
  | _ ⇒
    let m := F.StoZ (F.mag x) in
    let prec := F.incr_prec prec (Z2P (m + 6)) in
    snd (sin_cos_reduce prec x (S (Z2nat m)))
  end.

Lemma cos_fastP_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x)%R →
  contains (I.convert (cos_fastP prec x)) (Xreal (cos (toR x))).

Definition cos_fast prec x :=
  match F'.cmp x F.zero with
  | Xeq ⇒ i1
  | Xlt ⇒ cos_fastP prec (F.neg x)
  | Xgt ⇒ cos_fastP prec x
  | Xund ⇒ I.nai
  end.

Theorem cos_fast_correct :
  ∀ prec x,
  contains (I.convert (cos_fast prec x)) (Xcos (F.toX x)).

Definition sin_fast0 prec x :=
  let xi := I.bnd x x in
  let x2i := I.sqr prec xi in
  let p := F.prec prec in
  let thre := F.scale c1 (F.ZtoS (Zneg p)) in
  let rem := cos_fast0_aux prec thre i1 x2i i4 i6 (nat_of_P p) in
  I.mul prec (I.sub prec i1 rem) xi.

Lemma sin_fast0_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (Rabs (toR x) ≤ /2)%R →
  contains (I.convert (sin_fast0 prec x)) (Xreal (sin (toR x))).

Definition sin_fastP prec x :=
  let th := c1_2 in
  match F'.le' x th with
  | true ⇒ sin_fast0 (F.incr_prec prec 1) x
  | _ ⇒
    let m := F.StoZ (F.mag x) in
    let prec := F.incr_prec prec (Z2P (m + 6)) in
    match sin_cos_reduce prec x (S (Z2nat m)) with
    | (s, c) ⇒
      let v := I.sqrt prec (I.sub prec i1 (I.sqr prec c)) in
      match s with
      | Lt ⇒ I.neg v
      | Gt ⇒ v
      | _ ⇒ I.bnd cm1 c1
      end
    end
  end.

Lemma sin_fastP_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x)%R →
  contains (I.convert (sin_fastP prec x)) (Xreal (sin (toR x))).

Definition sin_fast prec x :=
  match F'.cmp x F.zero with
  | Xeq ⇒ I.zero
  | Xlt ⇒ I.neg (sin_fastP prec (F.neg x))
  | Xgt ⇒ sin_fastP prec x
  | Xund ⇒ I.nai
  end.

Theorem sin_fast_correct :
  ∀ prec x,
  contains (I.convert (sin_fast prec x)) (Xsin (F.toX x)).

Definition tan_fastP prec x :=
  let th := c1_2 in
  match F'.le' x th with
  | true ⇒
    let prec := F.incr_prec prec 2 in
    let s := sin_fast0 prec x in
    I.div prec s (I.sqrt prec (I.sub prec i1 (I.sqr prec s)))
  | _ ⇒
    let m := F.StoZ (F.mag x) in
    let prec := F.incr_prec prec (Z2P (m + 7)) in
    match sin_cos_reduce prec x (S (Z2nat m)) with
    | (s, c) ⇒
      let v := I.sqrt prec (I.sub prec (I.inv prec (I.sqr prec c)) i1) in
      match s, I.sign_large c with
      | Lt, Xgt ⇒ I.neg v
      | Gt, Xlt ⇒ I.neg v
      | Lt, Xlt ⇒ v
      | Gt, Xgt ⇒ v
      | _, _ ⇒ I.nai
      end
    end
  end.

Definition tan_fast prec x :=
  match F'.cmp x F.zero with
  | Xeq ⇒ I.zero
  | Xlt ⇒ I.neg (tan_fastP prec (F.neg x))
  | Xgt ⇒ tan_fastP prec x
  | Xund ⇒ I.nai
  end.

Lemma tan_fastP_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x)%R →
  contains (I.convert (tan_fastP prec x)) (Xtan (Xreal (toR x))).

Theorem tan_fast_correct :
  ∀ prec x,
  contains (I.convert (tan_fast prec x)) (Xtan (F.toX x)).

Definition semi_extension f fi :=
  ∀ x, contains (I.convert (fi x)) (f (F.toX x)).

Definition cos_correct : ∀ prec, semi_extension Xcos (cos_fast prec) := cos_fast_correct.
Definition sin_correct : ∀ prec, semi_extension Xsin (sin_fast prec) := sin_fast_correct.
Definition tan_correct : ∀ prec, semi_extension Xtan (tan_fast prec) := tan_fast_correct.
Definition atan_correct : ∀ prec, semi_extension Xatan (atan_fast prec) := atan_fast_correct.

Fixpoint expn_fast0_aux prec thre powi xi facti divi (nb : nat) { struct nb } :=
  let npwi := I.mul prec powi xi in
  let vali := I.div prec npwi divi in
  match F'.cmp (I.upper vali) thre, nb with
  | Xlt, _
  | _, O ⇒ I.bnd F.zero (I.upper vali)
  | _, S n ⇒
    let nfacti := I.add prec facti i1 in
    let ndivi := I.mul prec divi facti in
    I.sub prec vali (expn_fast0_aux prec thre npwi xi nfacti ndivi n)
  end.

Definition expn_fast0 prec x :=
  let p := F.prec prec in
  let thre := F.scale c1 (F.ZtoS (Zneg p)) in
  let xi := I.bnd x x in
  let rem := expn_fast0_aux prec thre xi xi i3 i2 (nat_of_P p) in
  I.sub prec i1 (I.sub prec xi rem).

Definition expn_reduce prec x :=
  let th := c1_8 in
  match F'.le' x th with
  | true ⇒ expn_fast0 (F.incr_prec prec 1) x
  | false ⇒
    let m := F.StoZ (F.mag x) in
    let prec := F.incr_prec prec (Z2P (9 + m)) in
    let fix reduce x (nb : nat) {struct nb} :=
      match F'.le' x th, nb with
      | true, _ ⇒ expn_fast0 prec x
      | _, O ⇒ I.bnd F.zero c1
      | _, S n ⇒ I.sqr prec (reduce (F.div2 x) n)
      end in
    reduce x (8 + Z2nat m)
  end.

Definition exp_fast prec x :=
  match F'.cmp x F.zero with
  | Xeq ⇒ i1
  | Xlt ⇒ expn_reduce prec (F.neg x)
  | Xgt ⇒
    let prec := F.incr_prec prec 1 in
    match I.invnz prec (expn_reduce prec x) with
    | Ibnd _ _ as b ⇒ b
    | Inai ⇒ I.bnd c1 F.nan
    end
  | Xund ⇒ I.nai
  end.

Lemma expn_fast0_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 ≤ toR x ≤ /2)%R →
  contains (I.convert (expn_fast0 prec x)) (Xreal (exp (- toR x))).

Lemma expn_reduce_correct :
  ∀ prec x,
  F.toX x = Xreal (toR x) →
  (0 < toR x)%R →
  contains (I.convert (expn_reduce prec x)) (Xreal (exp (- toR x))).

Theorem exp_fast_correct :
  ∀ prec x,
  contains (I.convert (exp_fast prec x)) (Xexp (F.toX x)).

End TranscendentalFloatFast.