Library Interval.Poly.Taylor_model

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) 2013-2016, 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 ZArith Psatz.
From mathcomp.ssreflect Require Import ssreflect ssrbool ssrfun eqtype ssrnat seq bigop.

Require Import Interval.
Require Import Xreal.
Require Import Basic.
Require Import Interval_compl.
Require Import Datatypes.
Require Import Taylor_model_sharp.
Require Import Bound.
Require Import Bound_quad.
Require Import Univariate_sig.

Set Implicit Arguments.

Auxiliary lemmas on natural numbers


Local Open Scope nat_scope.

Lemma maxnS (m n : nat) : 0 < maxn m n.+1.

Lemma maxSn (m n : nat) : 0 < maxn m.+1 n.

Parameterized Module for Taylor Models


Module TM (I : IntervalOps) <: UnivariateApprox I.

Load the CoqApprox modules

Module Pol := SeqPolyInt I.
Module Bnd := PolyBoundHornerQuad I Pol.
Module Import TMI := TaylorModel I Pol Bnd.

Main type definitions


Inductive t_ := Const of I.type | Var | Tm of rpa.

Definition T := t_.

Definition U := (I.precision × nat)%type.

Definition i1 := I.fromZ_small 1.
Definition i2 := I.fromZ_small 2.
Definition im1 := I.fromZ_small (-1).

Auxiliary material


Definition tmsize (tm : rpa) := Pol.size (approx tm).

Definition tsize (t : T) : nat :=
  match t with
    | Const _ ⇒ 1
    | Var ⇒ 2
    | Tm tm ⇒ tmsize tm
  end.

Definition get_tm (u : U) X (t : T) : rpa :=
  match t with
    | Const c ⇒ TM_cst X c
    | Var ⇒ TM_var X (J.midpoint X)
    | Tm tm ⇒ tm
  end.

Lemma size_get_tm u X t :
  tmsize (get_tm u X t) = tsize t.

Define the main validity predicate


Definition approximates (X : I.type) (tf : T) (f : R → ExtendedR) : Prop :=
  not_empty (I.convert X) →
  match tf with
  | Const c ⇒ I.convert c = Inan ∨ is_const f (I.convert X) (I.convert c)
  | Var ⇒
    ∀ x : R, contains (I.convert X) (Xreal x) → f x = Xreal x
  | Tm tm ⇒
    let x0 := proj_val (I.F.convert (I.midpoint X)) in
    i_validTM x0 (I.convert X) tm f
  end.

Theorem approximates_ext f g xi t :
  (∀ x, f x = g x) →
  approximates xi t f → approximates xi t g.

Lemma contains_midpoint :
  ∀ X : I.type,
  not_empty (I.convert X) →
  contains (I.convert X) (Xreal (proj_val (I.F.convert (I.midpoint X)))).

Lemma get_tm_correct u Y tf f :
  approximates Y tf f →
  approximates Y (Tm (get_tm u Y tf)) f.

Main definitions and correctness claims


Definition const : I.type → T := Const.

Theorem const_correct (c : I.type) (r : R) :
  contains (I.convert c) (Xreal r) →
  ∀ (X : I.type),
  approximates X (const c) (fun _ ⇒ Xreal r).

Definition dummy : T := Const I.nai.

Theorem dummy_correct xi f :
  TM.approximates xi dummy f.

Definition var : T := Var.

Theorem var_correct (X : I.type) :
  approximates X var Xreal.

Definition eval (u : U) (t : T) (Y X : I.type) : I.type :=
  if I.subset X Y then
  match t with
  | Const c ⇒ I.mask c X
  | Var ⇒ X
  | Tm tm ⇒
    let X0 := J.midpoint Y in
    let tm := get_tm u Y t in
    I.add u.1
      (Bnd.ComputeBound u.1 (approx tm) (I.sub u.1 X X0))
      (error tm)
  end
  else I.nai.

Theorem eval_correct u (Y : I.type) tf f :
  approximates Y tf f → I.extension (Xbind f) (eval u tf Y).

Definition add_slow (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  let M1 := get_tm u X t1 in
  let M2 := get_tm u X t2 in
  Tm (TM_add u.1 M1 M2).

Definition add (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  match t1, t2 with
    | Const c1, Const c2 ⇒ Const (I.add u.1 c1 c2)
    | _, _ ⇒ add_slow u X t1 t2
  end.

Lemma add_slow_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (add_slow u Y tf tg) (fun x ⇒ Xadd (f x) (g x)).

Theorem add_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (add u Y tf tg) (fun x ⇒ Xadd (f x) (g x)).

Definition opp_slow (u : U) (X : I.type) (t : T) : T :=
  Tm (TM_opp (get_tm u X t)).

Definition opp (u : U) (X : I.type) (t : T) : T :=
  match t with
    | Const c ⇒ Const (I.neg c)
    | _ ⇒ opp_slow u X t
  end.

Lemma opp_slow_correct u (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (opp_slow u Y tf) (fun x ⇒ Xneg (f x)).

Theorem opp_correct u (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (opp u Y tf) (fun x ⇒ Xneg (f x)).

Definition sub_slow (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  let M1 := get_tm u X t1 in
  let M2 := get_tm u X t2 in
  Tm (TM_sub u.1 M1 M2).

Definition sub (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  match t1, t2 with
    | Const c1, Const c2 ⇒ Const (I.sub u.1 c1 c2)
  
    | _, _ ⇒ sub_slow u X t1 t2
  end.

Lemma sub_slow_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (sub_slow u Y tf tg) (fun x ⇒ Xsub (f x) (g x)).

Theorem sub_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (sub u Y tf tg) (fun x ⇒ Xsub (f x) (g x)).

Definition mul_slow (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  let M1 := get_tm u X t1 in
  let M2 := get_tm u X t2 in
  let X0 := J.midpoint X in
  Tm (TM_mul u.1 M1 M2 X0 X u.2).

Definition mul (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  match t1, t2 with
    | Const c1, Const c2 ⇒ Const (I.mul u.1 c1 c2)
    | Const c1, _ ⇒ Tm (TM_mul_mixed u.1 c1 (get_tm u X t2) )
    | _, Const c2 ⇒ Tm (TM_mul_mixed u.1 c2 (get_tm u X t1) )
    | _, _ ⇒ mul_slow u X t1 t2
  end.

Lemma mul_slow_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (mul_slow u Y tf tg) (fun x ⇒ Xmul (f x) (g x)).

Theorem mul_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (mul u Y tf tg) (fun x ⇒ Xmul (f x) (g x)).

Definition div_slow (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  let M1 := get_tm u X t1 in
  let M2 := get_tm u X t2 in
  let X0 := J.midpoint X in
  Tm (TM_div u.1 M1 M2 X0 X u.2).

Definition div (u : U) (X : I.type) (t1 : T) (t2 : T) : T :=
  match t1, t2 with
    | Const c1, Const c2 ⇒ Const (I.div u.1 c1 c2)
    | _, Const c2 ⇒ Tm (TM_div_mixed_r u.1 (get_tm u X t1) c2)
  
    | _, _ ⇒ div_slow u X t1 t2
  end.

Lemma div_slow_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (div_slow u Y tf tg) (fun x ⇒ Xdiv (f x) (g x)).

Theorem div_correct u (Y : I.type) tf tg f g :
  approximates Y tf f → approximates Y tg g →
  approximates Y (div u Y tf tg) (fun x ⇒ Xdiv (f x) (g x)).

Definition abs (u : U) (X : I.type) (t : T) : T :=
  let e := eval u t X X in
  match I.sign_large e with
  | Xeq | Xgt ⇒ t
  | Xlt ⇒ opp u X t
  | Xund ⇒ Tm (TM_any u.1 (I.abs e) X u.2)
  end.

Lemma Isign_large_Xabs (u : U) (tf : T) (Y X : I.type) f :
  approximates Y tf f →
  match I.sign_large (eval u tf Y X) with
    | Xeq ⇒
      ∀ x, contains (I.convert X) (Xreal x) →
      f x = Xabs (f x)
    | Xgt ⇒
      ∀ x, contains (I.convert X) (Xreal x) →
      f x = Xabs (f x)
    | Xlt ⇒
      ∀ x, contains (I.convert X) (Xreal x) →
      Xneg (f x) = Xabs (f x)
    | Xund ⇒ True
  end.

Local Ltac byp a b :=
  move⇒ x Hx; rewrite -a //; exact: b.
Local Ltac byp' a b := let Hne := fresh in
  move⇒ Hne x Hx; rewrite -a //; exact: (b Hne).
Local Ltac foo :=
  by move⇒ Hne; apply: TM_any_correct;
  [ exact: contains_midpoint
  | intros x Hx; apply: I.abs_correct; exact: eval_correct Hx].

Theorem abs_correct u (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (abs u Y tf) (fun x ⇒ Xabs (f x)).

Definition Iconst (i : I.type) :=
  I.bounded i && I.subset i (I.bnd (I.lower i) (I.lower i)).

Lemma Iconst_true_correct i x :
  Iconst i = true → contains (I.convert i) x → x = Xlower (I.convert i).

Definition nearbyint (m : rounding_mode)
     (u : U) (X : I.type) (t : T) : T :=
  let e := eval u t X X in
  let e1 := I.nearbyint m e in
  if Iconst e1 then Const (I.bnd (I.lower e1) (I.lower e1))
  else
  let (p, i) :=
     match m with
     | rnd_UP ⇒
         let vm1 := I.lower im1 in
         let v1 := I.upper i1 in
         let i := I.bnd vm1 v1 in
         (I.div u.1 i1 i2, I.div u.1 i i2)
     | rnd_DN ⇒
         let vm1 := I.lower im1 in
         let v1 := I.upper i1 in
         let i := I.bnd vm1 v1 in
         (I.div u.1 im1 i2, I.div u.1 i i2)
     | rnd_ZR ⇒
          match I.sign_large e with
         | Xlt ⇒
             let vm1 := I.lower im1 in
             let v1 := I.upper i1 in
             let i := I.bnd vm1 v1 in
             (I.div u.1 i1 i2, I.div u.1 i i2)
         | Xund ⇒
            let vm1 := I.lower im1 in
            let v1 := I.upper i1 in
            (I.zero, I.bnd vm1 v1)
         | _ ⇒
             let vm1 := I.lower im1 in
             let v1 := I.upper i1 in
             let i := I.bnd vm1 v1 in
             (I.div u.1 im1 i2, I.div u.1 i i2)
         end
     | rnd_NE ⇒
         let vm1 := I.lower im1 in
         let v1 := I.upper i1 in
         let i := I.bnd vm1 v1 in
         (I.zero, I.div u.1 i i2)
     end in
  add u X t (Tm {|approx := Pol.polyC p;
                  error := I.mask i e1 |}).

Lemma contains_fromZ_lower_upper z1 z2 i :
  (-256 ≤ z1 ≤ 0)%Z → (0 ≤ z2 ≤ 256)%Z →
  contains
  (I.convert
     (I.mask (I.bnd (I.lower (I.fromZ_small z1)) (I.upper (I.fromZ_small z2))) i))
  (Xreal 0).

Lemma contains_fromZ_lower_upper_div prec z1 z2 z3 i :
  (-256 ≤ z1 ≤ 0)%Z → (0 ≤ z2 ≤ 256)%Z → (0 < z3 ≤ 256)%Z →
  contains
  (I.convert
     (I.mask
       (I.div prec
          (I.bnd (I.lower (I.fromZ_small z1)) (I.upper (I.fromZ_small z2)))
           (I.fromZ_small z3)) i))
    (Xreal 0).

Theorem nearbyint_correct m u (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (nearbyint m u Y tf) (fun x ⇒ Xnearbyint m (f x)).

Definition error_flt (u : U) (m : rounding_mode) (emin : Z) (prec : positive)
     (X : I.type) (t : T) : T :=
  let e := eval u t X X in
  let err := I.error_flt u.1 m emin prec e in
  Tm (TM_any u.1 err X 0).

Theorem error_flt_correct u m e p (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (error_flt u m e p Y tf) (fun x ⇒ Xerror_flt m e p (f x)).

Definition round_flt (u : U) (m : rounding_mode) (emin : Z) (prec : positive)
     (X : I.type) (t : T) : T :=
  add u X t (error_flt u m emin prec X t).

Theorem round_flt_correct u m e p (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (round_flt u m e p Y tf) (fun x ⇒ Xround_flt m e p (f x)).

Definition error_fix (u : U) (m : rounding_mode) (emin : Z) (X : I.type) (t : T) : T :=
  let e := eval u t X X in
  let err := I.error_fix u.1 m emin e in
  Tm (TM_any u.1 err X 0).

Theorem error_fix_correct u m e (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (error_fix u m e Y tf) (fun x ⇒ Xerror_fix m e (f x)).

Definition round_fix (u : U) (m : rounding_mode) (emin : Z) (X : I.type) (t : T) : T :=
  add u X t (error_fix u m emin X t).

Theorem round_fix_correct u m e (Y : I.type) tf f :
  approximates Y tf f →
  approximates Y (round_fix u m e Y tf) (fun x ⇒ Xround_fix m e (f x)).

Generic implementation of basic functions


Definition fun_gen
  (fi : I.precision → I.type → I.type)
  (ftm : I.precision → TM_type)
  (u : U) (X : I.type) (t : T) : T :=
  match t with
    | Const c ⇒ Const (fi u.1 c)
    | Var ⇒ let X0 := J.midpoint X in Tm (ftm u.1 X0 X u.2)
    | Tm tm ⇒ let X0 := J.midpoint X in
      Tm (TM_comp u.1 (ftm u.1) tm X0 X u.2)
  end.

Lemma fun_gen_correct
  (fi : I.precision → I.type → I.type)
  (ftm : I.precision → TM_type)
  (fx : R → ExtendedR)
  (ft := fun_gen fi ftm) :
  (∀ prec : I.precision, I.extension (Xbind fx) (fi prec)) →
  (∀ prec X0 X n, tmsize (ftm prec X0 X n) = n.+1) →
  (∀ prec x0 X0 X n,
    subset' (I.convert X0) (I.convert X) →
    contains (I.convert X0) (Xreal x0) →
    i_validTM x0 (I.convert X) (ftm prec X0 X n) fx) →
  ∀ (u : U) (X : I.type) (tf : T) (f : R → ExtendedR),
  approximates X tf f →
  approximates X (ft u X tf) (fun x ⇒ Xbind fx (f x)).



Definition inv := Eval hnf in fun_gen I.inv TM_inv.

Theorem inv_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (inv u Y tf) (fun x ⇒ Xinv (f x)).

Definition sqrt := Eval hnf in fun_gen I.sqrt TM_sqrt.

Theorem sqrt_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (sqrt u Y tf) (fun x ⇒ Xsqrt (f x)).

Definition exp := Eval hnf in fun_gen I.exp TM_exp.

Theorem exp_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (exp u Y tf) (fun x ⇒ Xexp (f x)).

Definition ln := Eval hnf in fun_gen I.ln TM_ln.

Theorem ln_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (ln u Y tf) (fun x ⇒ Xln (f x)).

Definition cos := Eval hnf in fun_gen I.cos TM_cos.

Theorem cos_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (cos u Y tf) (fun x ⇒ Xcos (f x)).

Definition sin := Eval hnf in fun_gen I.sin TM_sin.

Theorem sin_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (sin u Y tf) (fun x ⇒ Xsin (f x)).

Definition tan := Eval hnf in fun_gen I.tan TM_tan.

Theorem tan_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (tan u Y tf) (fun x ⇒ Xtan (f x)).

Definition atan := Eval hnf in fun_gen I.atan TM_atan.

Theorem atan_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (atan u Y tf) (fun x ⇒ Xatan (f x)).

Definition power_int p := Eval cbv delta[fun_gen] beta in
  match p with

  | 1%Z ⇒ fun u xi t ⇒ t
  | _ ⇒ fun_gen
           (fun prec x ⇒ I.power_int prec x p)
           (fun prec ⇒ TM_power_int prec p)
  end.

Theorem power_int_correct :
  ∀ (p : Z) u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (power_int p u Y tf) (fun x ⇒ Xbind (fun y ⇒ Xpower_int (Xreal y) p) (f x)).

Definition sqr := power_int 2.

Theorem sqr_correct :
  ∀ u (Y : I.type) tf f,
  approximates Y tf f →
  approximates Y (sqr u Y tf) (fun x ⇒ Xsqr (f x)).

End TM.