Library CoqApprox.taylor_model_int_sharp

This file is part of the CoqApprox formalization of rigorous polynomial approximation in Coq: http://tamadi.gforge.inria.fr/CoqApprox/
Copyright (c) 2010-2013, ENS de Lyon and 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.

Require Import ZArith Psatz.
Require Import Fcore_Raux.
Require Import Interval_xreal.
Require Import Interval_generic Interval_interval.
Require Import Interval_definitions.
Require Import Interval_specific_ops.
Require Import Interval_float_sig.
Require Import Interval_interval_float.
Require Import Interval_interval_float_full.
Require Import Interval_xreal.
Require Import Interval_xreal_derive.
Require Import Interval_missing.
Require Import Interval_generic_proof.
Require Import Rstruct Classical.
Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq fintype bigop ssralg.
Require Import xreal_ssr_compat.
Require Import seq_compl.
Require Import interval_compl.
Require Import poly_datatypes.
Require Import taylor_poly.
Require Import coeff_inst.
Require Import rpa_inst.
Require Import derive_compl.

This theory implements Taylor Models, with sharp remainders for univariate base functions, using the so-called Zumkeller's technique which is partly based on the standard enclosure of the Taylor/Lagrange remainder.

Set Implicit Arguments.

Import GRing.Theory.

Local Open Scope nat_scope.

Module TaylorModel
  (I : IntervalOps)
  (Import Pol : IntMonomPolyOps I)
  (PolX : ExactMonomPolyOps FullXR)
  (Link : LinkIntX I Pol PolX).
Module Import RPA := RigPolyApproxInt I Pol.
Module Import TI := TaylorPoly Pol.Int Pol.
Import Int.
Module TX := TaylorPoly FullXR PolX.




Section TaylorModel.
The presence of this section variable does not hinder any sophisticated handling of the precision inside the functions that are defined or called below.
Variable prec : I.precision.
Variable M : rpa.

Variable Tcoeffs : T -> nat -> Pol.T.
For complexity reasons, Tcoeffs must return more than one coefficient.
The generic functions TLrem/Ztech are intended to ease the computation of the interval remainder for basic functions.
Definition TLrem X0 X n :=
  let N := S n in
  let NthCoeff := Pol.tnth (Tcoeffs X N) N in
  let NthPower :=
    I.power_int prec (I.sub prec X X0) (Z_of_nat N) in
  I.mul prec NthCoeff NthPower.

The following predicate, involved in Ztech below, will contribute to the conciseness of the Coq goals, notably when issuing tactic rewrite /Ztech.
Definition isNNegOrNPos (X : I.type) : bool :=
  if I.sign_large X is Xund then false else true.

The first argument of Ztech will be instantiated with Tcoeffs X0 n.
Definition Ztech (P : Pol.T) F X0 X n :=
  let N := S n in
  let NthCoeff := Pol.tnth (Tcoeffs X N) N in
  if isNNegOrNPos NthCoeff && I.bounded X then
    let a := I.lower X in let b := I.upper X in
    let A := I.bnd a a in let B := I.bnd b b in
    let Da := I.sub prec (F A) (teval prec P (I.sub prec A X0)) in
    let Db := I.sub prec (F B) (teval prec P (I.sub prec B X0)) in
    let Dx0 := I.sub prec (F X0) (teval prec P (I.sub prec X0 X0)) in
    I.join (I.join Da Db) Dx0
  else
    let NthPower :=
      I.power_int prec (I.sub prec X X0) (Z_of_nat N) in
    I.mul prec NthCoeff NthPower.

Lemma ZtechE1 P F X0 X n :
  isNNegOrNPos (Pol.tnth (Tcoeffs X n.+1) n.+1) = false ->
  Ztech P F X0 X n = TLrem X0 X n.

Lemma ZtechE2 P F X0 X n :
  I.bounded X = false ->
  Ztech P F X0 X n = TLrem X0 X n.


End TaylorModel.

Section PrecArgument.
The presence of this section variable does not hinder any sophisticated handling of the precision inside the functions that are defined or called below.
Variable prec : I.precision.

Note that Zumkeller's technique is not necessary for TM_cst & TM_var.
Definition TM_cst c (X0 X : I.type) (n : nat) : rpa :=
  RPA (T_cst c X0 n) (TLrem prec (T_cst c) X0 X n).

Definition TM_var X0 X (n : nat) :=
  RPA (T_var X0 n) (TLrem prec T_var X0 X n).

Definition TM_inv X0 X (n : nat) :=
  let P := (T_inv prec X0 n) in
  
Note that this let-in is essential in call-by-value context.
  RPA P (Ztech prec (T_inv prec) P (tinv prec) X0 X n).

Definition TM_exp X0 X (n : nat) : rpa :=
  let P := (T_exp prec X0 n) in
  RPA P (Ztech prec (T_exp prec) P (texp prec) X0 X n).

Definition TM_sqrt X0 X (n : nat) : rpa :=
  let P := (T_sqrt prec X0 n) in
  RPA P (Ztech prec (T_sqrt prec) P (tsqrt prec) X0 X n).

Definition TM_invsqrt X0 X (n : nat) : rpa :=
  let P := (T_invsqrt prec X0 n) in
  RPA P (Ztech prec (T_invsqrt prec) P (tinvsqrt prec) X0 X n).

Definition TM_sin X0 X (n : nat) : rpa :=
  let P := (T_sin prec X0 n) in
  RPA P (Ztech prec (T_sin prec) P (tsin prec) X0 X n).

Definition TM_cos X0 X (n : nat) : rpa :=
  let P := (T_cos prec X0 n) in
  RPA P (Ztech prec (T_cos prec) P (tcos prec) X0 X n).

Definition TM_add (Mf Mg : rpa) : rpa :=
  RPA (Pol.tadd prec (approx Mf) (approx Mg))
    (I.add prec (error Mf) (error Mg)).


Definition i_validTM (X0 X : interval)
  (M : rpa) (f : ExtendedR -> ExtendedR) :=
  [/\ contains (I.convert (error M)) (Xreal R0),
    I.subset_ X0 X &
    let N := Pol.tsize (approx M) in
    forall fi0, contains X0 fi0 ->
    exists alf,
    [/\ PolX.tsize alf = N,
      forall k, (k < N)%nat ->
        contains (I.convert (Pol.tnth (approx M) k)) (PolX.tnth alf k) &
      forall x, contains X x -> contains (I.convert (error M))
        (FullXR.tsub tt (f x) (PolX.teval tt alf (FullXR.tsub tt x fi0)))]].

Local Notation Ibnd2 x := (I.bnd x x) (only parsing).

Some extra support lemmas about intervals.

Lemma mul_0_contains_0_l y Y X :
  contains (I.convert Y) y ->
  contains (I.convert X) (Xreal 0) ->
  contains (I.convert (I.mul prec X Y)) (Xreal 0).

Lemma mul_0_contains_0_r y Y X :
  contains (I.convert Y) y ->
  contains (I.convert X) (Xreal 0) ->
  contains (I.convert (I.mul prec Y X)) (Xreal 0).

Lemma pow_contains_0 (X : I.type) (n : Z) :
  (n > 0)%Z ->
  contains (I.convert X) (Xreal 0) ->
  contains (I.convert (I.power_int prec X n)) (Xreal 0).

Lemma subset_sub_contains_0 x0 (X0 X : I.type) :
  contains (I.convert X0) x0 ->
  I.subset_ (I.convert X0) (I.convert X) ->
  contains (I.convert (I.sub prec X X0)) (Xreal 0).

Lemma subset_real_contains X rx ry c :
  contains (I.convert X) (Xreal rx) ->
  contains (I.convert X) (Xreal ry) -> (rx <= c <= ry)%Re ->
  contains (I.convert X) (Xreal c).

Lemma Isub_Inan_propagate_l y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.sub prec X Y) = Interval_interval.Inan.

Lemma Isub_Inan_propagate_r y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.sub prec Y X) = Interval_interval.Inan.

Lemma Imul_Inan_propagate_l y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.mul prec X Y) = Interval_interval.Inan.

Lemma Imul_Inan_propagate_r y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.mul prec Y X) = Interval_interval.Inan.

Lemma Ipow_Inan_propagate (X : I.type) (n : Z) :
  I.convert X = Interval_interval.Inan ->
  I.convert (I.power_int prec X n) = Interval_interval.Inan.

Lemma Iadd_Inan_propagate_l y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.add prec X Y) = Interval_interval.Inan.

Lemma Iadd_Inan_propagate_r y Y X :
  contains (I.convert Y) y ->
  I.convert X = Interval_interval.Inan ->
  I.convert (I.add prec Y X) = Interval_interval.Inan.

Lemma Iadd_zero_subset_l (a b : I.type) :
  (exists t, contains (I.convert a) t) ->
  contains (I.convert b) (Xreal 0) ->
  I.subset_ (I.convert a) (I.convert (I.add prec b a)).

Lemma Iadd_zero_subset_r a b :
  (exists t, contains (I.convert a) t) ->
  contains (I.convert b) (Xreal 0) ->
  I.subset_ (I.convert a) (I.convert (I.add prec a b)).

Lemma Iadd_Isub_aux b a B D :
  contains (I.convert B) b ->
  contains (I.convert D) (a - b)%XR ->
  contains (I.convert (I.add prec B D)) a.

Some auxiliary lemmas related to polynomials.

Lemma is_horner_pos (p : PolX.T) (x : ExtendedR) :
  0 < PolX.tsize p ->
  PolX.teval tt p x =
  \big[Xadd/FullXR.tzero]_(i < PolX.tsize p)
    FullXR.tmul tt (PolX.tnth p i)(FullXR.tpow tt x i).

Lemma tpolyNil_size0 p :
  PolX.tsize p = 0 -> p = PolX.tpolyNil.

Lemma is_horner_mask (p : PolX.T) (x : FullXR.T) :
  PolX.teval tt p x =
  Xmask (\big[Xadd/Xreal 0]_(i < PolX.tsize p)
    FullXR.tmul tt (PolX.tnth p i)(FullXR.tpow tt x i)) x.

Another lemma that will be useful to prove the correctness of TM_comp.

Lemma i_validTM_subset_X0 (X0 : I.type) (SX0 : interval) (X : I.type) f Mf :
  I.subset_ SX0 (I.convert X0) ->
  i_validTM (I.convert X0) (I.convert X) Mf f ->
  i_validTM (SX0) (I.convert X) Mf f.

Section ProofOfRec.

Variable XF0 : FullXR.T -> FullXR.T.
Variable XDn : nat -> FullXR.T -> FullXR.T.
Class nth_Xderive := Nth_Xderive : forall x:FullXR.T, nth_Xderive_pt XF0 XDn x.

Section Corollaries.
Context { XDn_ : nth_Xderive}.

Lemma Xder n x : Xderive_pt (XDn n) x (XDn (S n) x).

Lemma XDn_0 x : XDn 0 x = XF0 x.

Lemma XDn_Xnan k i x : XDn i x = Xnan -> XDn (k+i) x = Xnan.
End Corollaries.

Hypothesis XF0_Xnan : XF0 Xnan = Xnan.
Context { XDn_ : nth_Xderive}.
Variable F0 : T -> T.
Hypothesis F0_contains : I.extension XF0 F0.

Lemma XDn_Xnan' (n : nat) : XDn n Xnan = Xnan.

Section GenericProof.
Generic proof for TLrem/Ztech.

Variable XP : FullXR.T -> nat -> PolX.T.
Class validXPoly : Prop := ValidXPoly {
  XPoly_size : forall (xi0 : FullXR.T) n, PolX.tsize (XP xi0 n) = n.+1;
  XPoly_nth : forall (xi0 : FullXR.T) n k, k < n.+1 ->
    PolX.tnth (XP xi0 n) k = Xdiv (XDn k xi0) (Xreal (INR (fact k))) }.

Variable IP : T -> nat -> Pol.T.
Class validPoly : Prop := ValidPoly {
  Poly_size : forall (X0 : I.type) xi0 n,
    PolX.tsize (XP xi0 n) = tsize (IP X0 n);
  Poly_nth : forall (X0 : I.type) xi0 n k,
    contains (I.convert X0) xi0 -> k < n.+1 ->
    contains (I.convert (tnth (IP X0 n) k)) (PolX.tnth (XP xi0 n) k) }.

Context { validXPoly_ : validXPoly }.
Context { validPoly_ : validPoly }.

Lemma XPoly_nth0 x n : PolX.tnth (XP x n) 0 = XF0 x.

Lemma Poly_size' X0 n : tsize (IP X0 n) = n.+1.

Lemma not_and_impl : forall P Q : Prop, ~ (P /\ Q) -> P -> ~ Q.

Note that we define the following predicate using exists quantifiers, to be able to use not_ex_all_not later on, rather than not_all_ex_not.
Definition Xnan_ex (D : nat -> ExtendedR -> ExtendedR)
  (n : nat) (X : interval) : Prop :=
  exists k : 'I_n.+1,
  exists x : R, contains X (Xreal x) /\ D k (Xreal x) = Xnan.

Lemma Xnan_ex_or D n X : Xnan_ex D n X \/ ~ Xnan_ex D n X.

Lemma Xnan_ex_not D n X : ~ Xnan_ex D n X ->
  forall k : nat, forall x : R, k <= n ->
  contains X (Xreal x) ->
  D k (Xreal x) = Xreal (proj_val (D k (Xreal x))).

Theorem i_validTM_TLrem X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (IP X0 n) (TLrem prec IP X0 X n)) XF0.

Lemma teval_contains u fi fx (X : I.type) x :
  Link.contains_pointwise fi fx ->
  contains (I.convert X) x ->
  contains (I.convert (teval u fi X)) (PolX.teval tt fx x).

Lemma convert_teval pi px (n := tsize pi) :
  n = PolX.tsize px ->
  (forall k, k < n -> contains (I.convert (tnth pi k)) (PolX.tnth px k)) ->
  forall I x, contains (I.convert I) x ->
  contains (I.convert (teval prec pi I)) (PolX.teval tt px x).

Lemma in_poly_bound p X pr :
  PolX.tsize pr = tsize p ->
  (forall i, i < tsize p -> contains (I.convert (tnth p i)) (PolX.tnth pr i)) ->
  contains (I.convert X) (Xreal 0) ->
  I.subset_ (I.convert (tnth p 0)) (I.convert (teval prec p X)).

Lemma le_upper_or (x y : ExtendedR) : le_upper x y \/ le_upper y x.

Lemma le_lower_or (x y : ExtendedR) : le_lower x y \/ le_lower y x.

Lemma Ilower_bnd (l u : I.bound_type) :
  I.convert_bound (I.lower (I.bnd l u)) = I.convert_bound l.

Lemma Iupper_bnd (l u : I.bound_type) :
  I.convert_bound (I.upper (I.bnd l u)) = I.convert_bound u.

Lemma Iupper_Xreal (X : I.type) (r : R) :
  I.convert_bound (I.upper X) = Xreal r -> I.convert X <> IInan.

Lemma Ilower_Xreal (X : I.type) (r : R) :
  I.convert_bound (I.lower X) = Xreal r -> I.convert X <> IInan.

Lemma upper_le (X : I.type) (x : ExtendedR ) :
  contains (I.convert X) x -> le_upper x (I.convert_bound (I.upper X)).

Lemma lower_le (X : I.type) (x : ExtendedR ) :
  contains (I.convert X) x -> le_lower (I.convert_bound (I.lower X)) x.

Lemma contains_lower_le_upper (X : interval) x :
  contains X x ->
  match (Xlower X), (Xupper X) with
  | Xreal r, Xreal s => (r <= s)%Re
  | _, _ => True
  end.

Lemma contains_IIbnd_Xreal (a b x : ExtendedR) :
  contains (IIbnd a b) x ->
  x = Xreal (proj_val x).

Lemma contains_lower_or_upper_Xreal (X : I.type) (xi0 : ExtendedR) (r : R) :
  contains (I.convert X) (Xreal r) -> contains (I.convert X) xi0 ->
  contains (IIbnd (I.convert_bound (I.lower X)) xi0) (Xreal r) \/
  contains (IIbnd xi0 (I.convert_bound (I.upper X))) (Xreal r).

Lemma contains_lower_or_upper (X : I.type) (xi0 : ExtendedR) (x : ExtendedR) :
  I.convert X <> IInan ->
  contains (I.convert X) x -> contains (I.convert X) xi0 ->
  contains (IIbnd (I.convert_bound (I.lower X)) xi0) x \/
  contains (IIbnd xi0 (I.convert_bound (I.upper X))) x.

Lemma isNNegOrNPos_false (X : I.type) :
  I.convert X = IInan -> isNNegOrNPos X = false.

Lemma teval_in_nan p : PolX.teval tt p Xnan = Xnan.

Lemma bounded_singleton_contains_lower_upper (X : I.type) :
  I.bounded X = true ->
  contains (I.convert (Ibnd2 (I.lower X))) (I.convert_bound (I.lower X)) /\
  contains (I.convert (Ibnd2 (I.upper X))) (I.convert_bound (I.upper X)).

Lemma bounded_IInan (X : I.type) :
  I.bounded X = true -> I.convert X <> IInan.

Lemma bounded_IIbnd (X : I.type) :
  I.bounded X = true -> I.convert X =
  IIbnd (I.convert_bound (I.lower X)) (I.convert_bound (I.upper X)).

Lemma bounded_Ilower (X : I.type) :
  I.bounded X = true -> I.convert_bound (I.lower X) =
  Xreal (proj_val (I.convert_bound (I.lower X))).

Lemma bounded_Iupper (X : I.type) :
  I.bounded X = true -> I.convert_bound (I.upper X) =
  Xreal (proj_val (I.convert_bound (I.upper X))).

Lemma bounded_lower_Xnan (X : I.type) :
  I.bounded X = true -> I.convert_bound (I.lower X) <> Xnan.

Lemma bounded_upper_Xnan (X : I.type) :
  I.bounded X = true -> I.convert_bound (I.upper X) <> Xnan.

Lemma bounded_contains_lower (X : I.type) (x : ExtendedR) :
  I.bounded X = true -> contains (I.convert X) x ->
  contains (I.convert X) (I.convert_bound (I.lower X)).


Lemma bounded_contains_upper (X : I.type) (x : ExtendedR) :
  I.bounded X = true -> contains (I.convert X) x ->
  contains (I.convert X) (I.convert_bound (I.upper X)).


Lemma contains_Xgt (a b : ExtendedR) :
  (exists x, contains (IIbnd a b) x) ->
  Xcmp a b <> Xgt.

Lemma Xund_imp (a b : ExtendedR) :
  Xcmp a b = Xund -> a = Xnan \/ b = Xnan.

Lemma Xund_contra (a b : ExtendedR) :
  a <> Xnan -> b <> Xnan -> Xcmp a b <> Xund.

Lemma Xeq_imp (a b : ExtendedR) :
  Xcmp a b = Xeq -> a = b.

Definition Xincr (f : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  forall x y : R,
  contains (IIbnd a b) (Xreal x) -> contains (IIbnd a b) (Xreal y) ->
  (x <= y -> forall vx (vy := vx), proj_fun vx f x <= proj_fun vy f y)%Re.

Definition Xdecr (f : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  forall x y : R,
  contains (IIbnd a b) (Xreal x) -> contains (IIbnd a b) (Xreal y) ->
  (x <= y -> forall vx (vy := vx), proj_fun vy f y <= proj_fun vx f x)%Re.

Definition Xmonot (f : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  Xincr f a b \/ Xdecr f a b.

Definition Xpos_over (g : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  forall x : R, contains (IIbnd a b) (Xreal x) ->
  forall v : R, (0 <= proj_fun v g x)%Re.

Definition Xneg_over (g : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  forall x : R, contains (IIbnd a b) (Xreal x) ->
  forall v : R, (proj_fun v g x <= 0)%Re.

Definition Xcst_sign (g : ExtendedR -> ExtendedR) (a b : ExtendedR) :=
  Xpos_over g a b \/ Xneg_over g a b.

Lemma eq_Xcst_sign (f g : ExtendedR -> ExtendedR) (a b : ExtendedR) :
  f =1 g -> Xcst_sign g a b -> Xcst_sign f a b.

Lemma eq'_Xcst_sign (f g : ExtendedR -> ExtendedR) (a b : ExtendedR) :
  (forall x, contains (IIbnd a b) x -> f x = g x) ->
  Xcst_sign g a b -> Xcst_sign f a b.

Definition Xderive_over (X : interval) (f f' : ExtendedR -> ExtendedR) :=
  Xderive_pt f Xnan (f' Xnan) /\
  forall x : ExtendedR, contains X x -> Xderive_pt f x (f' x).

Lemma Xderive_over_propagate
  (f f' : ExtendedR -> ExtendedR) (X : interval) (x : ExtendedR) :
  contains X x ->
  Xderive_over X f f' -> f x = Xnan -> f' x = Xnan.

Lemma Xderive_over_propagate' (f f' : ExtendedR -> ExtendedR) (X : interval) :
  Xderive_over X f f' -> f' Xnan = Xnan.

Lemma Xcmp_same (x : ExtendedR) :
  match Xcmp x x with
  | Xeq | Xund => True
  | _ => False
  end.

Lemma Xcmp_gt_imp_fun (f g : ExtendedR -> ExtendedR) (a b v : R) :
  (Xcmp (f(Xreal a)) (g(Xreal b)) = Xgt -> proj_fun v f a > proj_fun v g b)%Re.

Lemma Xcmp_lt_imp_fun (f g : ExtendedR -> ExtendedR) (a b v : R) :
  (Xcmp (f(Xreal a)) (g(Xreal b)) = Xlt -> proj_fun v f a < proj_fun v g b)%Re.

Lemma proj_fun_id (v x : R) : proj_fun v (@id ExtendedR) x = x.

Lemma contains_IIbnd (a b x : ExtendedR) (X : interval) :
  contains X a -> contains X b -> contains (IIbnd a b) x ->
  contains X x.

Definition Xreal_over (f' : ExtendedR -> ExtendedR) (a b : ExtendedR) : Prop :=
  forall x : R, contains (IIbnd a b) (Xreal x) ->
  f' (Xreal x) = Xreal (proj_val (f' (Xreal x))).

Definition Xreal_over' (f' : ExtendedR -> ExtendedR) (a b : ExtendedR) : Prop :=
  forall x : R, contains (IIbnd a b) (Xreal x) ->
  forall v : R, f' (Xreal x) = Xreal (proj_fun v f' x).

Lemma Xreal_over_imp (f' : ExtendedR -> ExtendedR) (a b : ExtendedR) :
  Xreal_over f' a b -> Xreal_over' f' a b.

Lemma Xderive_pos_imp_incr
  (f f' : ExtendedR -> ExtendedR) (X : interval) (a b : ExtendedR) :
  Xreal_over f' a b ->
  contains X a -> contains X b ->
  Xderive_over X f f' -> Xpos_over f' a b -> Xincr f a b.

Lemma Xderive_neg_imp_decr
  (f f' : ExtendedR -> ExtendedR) (X : interval) (a b : ExtendedR) :
  Xreal_over f' a b ->
  contains X a -> contains X b ->
  Xderive_over X f f' -> Xneg_over f' a b -> Xdecr f a b.

Lemma Xderive_cst_sign
  (f f' : ExtendedR -> ExtendedR) (X : interval) (a b : ExtendedR) :
  Xreal_over f' a b ->
  contains X a -> contains X b ->
  Xderive_over X f f' -> Xcst_sign f' a b -> Xmonot f a b.

Definition Xdelta (n : nat) (xi0 x : ExtendedR) :=
  Xsub (XF0 x) (PolX.teval tt (XP xi0 n) (Xsub x xi0)).

Definition Xdelta_big (n : nat) (xi0 x : ExtendedR) :=
  Xsub (XF0 x) (\big[Xadd/FullXR.tzero]_(i < n.+1) (PolX.tnth (XP xi0 n) i * (x - xi0) ^ i)%XR).

Definition Xdelta_mask (n : nat) (xi0 x : ExtendedR) :=
  Xsub (XF0 x) (Xmask (\big[Xadd/FullXR.tzero]_(i < n.+1) (PolX.tnth (XP xi0 n) i * (x - xi0) ^ i)%XR
) (x - xi0)%XR).

Lemma Xdelta_idem (n : nat) (xi0 x : ExtendedR) :
  Xdelta n xi0 x = Xdelta_big n xi0 x.

Lemma Xdelta_idem' (n : nat) (xi0 x : ExtendedR) :
  Xdelta n xi0 x = Xdelta_mask n xi0 x.

Now let's define the derivative of (Xdelta n xi0)
Definition Xdelta'_big (n : nat) (xi0 x : ExtendedR) :=
  Xsub (XDn 1 x) (Xmask (\big[Xadd/Xreal 0]_(i < n) (PolX.tnth (XP xi0 n) i.+1 *
   (Xreal (INR i.+1) * (x - xi0) ^ i)))%XR (x - xi0)%XR).

Lemma Xderive_pt_sub' f g f' g' x :
  Xderive_pt f x f' ->
  (f x - g x <> Xnan \/ f' - g' <> Xnan -> Xderive_pt g x g')%XR ->
  Xderive_pt (fun x => Xsub (f x) (g x)) x (Xsub f' g').

Lemma Xderive_pt_sub_r f g f' g' x :
  Xderive_pt f x f' ->
  (g x <> Xnan \/ g' <> Xnan -> Xderive_pt g x g') ->
  Xderive_pt (fun x => Xsub (f x) (g x)) x (Xsub f' g').

Lemma Xnan_ex_S D n X : Xnan_ex D n X -> Xnan_ex D n.+1 X.

Lemma derivable_pt_lim_pow_sub (r s : R) (n : nat) :
  derivable_pt_lim (fun x : R => (x - s) ^ n)%Re r (INR n * (r - s) ^ n.-1)%Re.

Lemma bigXadd_nS (n : nat) (xi0 : ExtendedR) (r : R) :
  \big[Xadd/Xreal 0]_(i < n) (PolX.tnth (XP xi0 n.+1) i.+1 * (Xreal (INR i.+1) * Xreal r ^ i))%XR =
  \big[Xadd/Xreal 0]_(i < n) (PolX.tnth (XP xi0 n) i.+1 * (Xreal (INR i.+1) * Xreal r ^ i))%XR.

Lemma bigXadd_ni (n : nat) (xi0 : ExtendedR) (r : R) :
  \big[Xadd/Xreal 0]_(i < n) (PolX.tnth (XP xi0 n.+1) i * Xreal r ^ i)%XR =
  \big[Xadd/Xreal 0]_(i < n) (PolX.tnth (XP xi0 n) i * Xreal r ^ i)%XR.

Lemma bigXadd_Si (n : nat) (xi0 : ExtendedR) (r : R) :
  \big[Xadd/Xreal 0]_(i < n.+1) (PolX.tnth (XP xi0 n.+1) i * Xreal r ^ i)%XR =
  \big[Xadd/Xreal 0]_(i < n.+1) (PolX.tnth (XP xi0 n) i * Xreal r ^ i)%XR.


Lemma bigXadd_XP (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  \big[Xadd/Xreal 0]_(i < n.+1)
    (PolX.tnth (XP xi0 n) i * Xreal s ^ i)%XR =
  \big[Xadd/Xreal 0]_(i < n.+1)
    ((XDn i xi0) / Xreal (INR (fact i)) * Xreal s ^ i)%XR.

Lemma bigXadd_real (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  ~ Xnan_ex XDn n X ->
  \big[Xadd/Xreal 0]_(i < n.+1)
    ((XDn i xi0) / Xreal (INR (fact i)) * Xreal s ^ i)%XR =
  Xreal (\big[Rplus/R0]_(i < n.+1)
    (proj_val (XDn i xi0) / (INR (fact i)) * s ^ i)%Re).

Lemma bigXadd'_XP_S (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  \big[Xadd/Xreal 0]_(i < n)
    (PolX.tnth (XP xi0 n) i.+1 * (Xreal (INR i.+1) * Xreal s ^ i))%XR =
  \big[Xadd/Xreal 0]_(i < n)
    ((XDn i.+1 xi0) / Xreal (INR (fact i.+1)) * (Xreal (INR i.+1) * Xreal s ^ i))%XR.

Lemma bigXadd'_real_S (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  ~ Xnan_ex XDn n X ->
  \big[Xadd/Xreal 0]_(i < n)
    ((XDn i.+1 xi0) / Xreal (INR (fact i.+1)) * (Xreal (INR i.+1) *
     Xreal s ^ i))%XR =
  Xreal (\big[Rplus/R0]_(i < n)
    (proj_val (XDn i.+1 xi0) / (INR (fact i.+1)) * ((INR i.+1) * s ^ i)))%Re.

Lemma bigXadd'_XP (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  \big[Xadd/Xreal 0]_(i < n)
    (PolX.tnth (XP xi0 n) i.+1 * (Xreal (INR i.+1) * Xreal s ^ i))%XR =
  \big[Xadd/Xreal 0]_(i < n)
    ((XDn i.+1 xi0) / Xreal (INR (fact i)) * Xreal s ^ i)%XR.

Lemma bigXadd'_real (n : nat) (X : interval) (r s : R) (xi0 := Xreal r) :
  contains X xi0 ->
  ~ Xnan_ex XDn n X ->
  \big[Xadd/Xreal 0]_(i < n)
    ((XDn i.+1 xi0) / Xreal (INR (fact i)) * Xreal s ^ i)%XR =
  Xreal (\big[Rplus/R0]_(i < n)
    (proj_val (XDn i.+1 xi0) / INR (fact i) * s ^ i))%Re.

Lemma Xderive_delta (n : nat) (xi0 : ExtendedR) (X : interval) :
  ~ Xnan_ex XDn n X ->
  contains X xi0 ->
  Xderive_over X (Xdelta n xi0) (Xdelta'_big n xi0).

Lemma Xdelta_real (n : nat) (X : interval) (r0 : R) (xi0 := Xreal r0) :
  contains X xi0 ->
  ~ Xnan_ex XDn n X ->
  forall a b : ExtendedR,
  contains X a -> contains X b ->
  Xreal_over (Xdelta n xi0) a b.

Lemma Xdelta'_real (n : nat) (X : interval) (r0 : R) (xi0 := Xreal r0) :
  contains X xi0 ->
  ~ Xnan_ex XDn n.+1 X ->
  forall a b : ExtendedR,
  contains X a -> contains X b ->
  Xreal_over (Xdelta'_big n xi0) a b.

Lemma notIInan_IIbnd (X : interval) :
  X <> IInan -> exists a : ExtendedR, exists b : ExtendedR, X = IIbnd a b.

Lemma Xmonot_contains
  (g : ExtendedR -> ExtendedR) (ra rb : R) (a := Xreal ra) (b := Xreal rb) :
  Xreal_over g a b ->
  Xmonot g a b ->
  forall (r s t : R) (x := Xreal r) (y := Xreal s) (z := Xreal t),
  contains (IIbnd a b) x -> contains (IIbnd a b) y ->
  contains (IIbnd x y) z ->
  contains (IIbnd (g x) (g y)) (g z) \/ contains (IIbnd (g y) (g x)) (g z).

Corollary Xmonot_contains_weak (g : ExtendedR -> ExtendedR) (a b : ExtendedR) :
  Xreal_over g a b ->
  Xmonot g a b ->
  a <> Xnan -> b <> Xnan ->
  forall x : ExtendedR, contains (IIbnd a b) x ->
  contains (IIbnd (g a) (g b)) (g x) \/ contains (IIbnd (g b) (g a)) (g x).

Lemma Rdiv_pos_compat (x y : R) :
  (0 <= x -> 0 < y -> 0 <= x / y)%Re.

Lemma Rlt_neq (x y : R) :
  (x < y -> x <> y)%Re.

Lemma Rlt_neq_sym (x y : R) :
  (x < y -> y <> x)%Re.

Lemma Rdiv_pos_compat_rev (x y : R) :
  (0 <= x / y -> 0 < y -> 0 <= x)%Re.

Lemma Rdiv_neg_compat (x y : R) :
  (x <= 0 -> 0 < y -> x / y <= 0)%Re.

Lemma Rdiv_neg_compat_rev (x y : R) :
  (x / y <= 0 -> 0 < y -> x <= 0)%Re.

Lemma Rdiv_twice (x y z : R) :
  y <> R0 -> z <> R0 -> (x / (y / z) = (x / y) * z)%Re.

Lemma upper_Xpos_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xpos_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd xi0 u) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (0 <= proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n)%Re.

Lemma upper_Xneg_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xneg_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd xi0 u) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n <= 0)%Re.

Lemma powerRZ_pow (r : R) (n : nat) :
  (powerRZ r (Z_of_nat n) = r ^ n)%Re.

Lemma pow_Rabs_sign (r : R) (n : nat) :
  (r ^ n = powerRZ
    (if Rle_bool R0 r then 1 else -1) (Z_of_nat n) * (Rabs r) ^ n)%Re.

Lemma powerRZ_1_even (k : Z) :(0 <= powerRZ (-1) (2 * k))%Re.

Lemma ZEven_pow_le (r : R) (n : nat) :
  Z.Even (Z_of_nat n) ->
  (0 <= r ^ n)%Re.

Lemma Ropp_le_0 (x : R) :
  (0 <= x -> - x <= 0)%Re.

Lemma ZOdd_pow_le (r : R) (n : nat) :
  Z.Odd (Z_of_nat n) ->
  (r <= 0 -> r ^ n <= 0)%Re.

Lemma lower_even_Xpos_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  Z.Even (Z_of_nat n) ->
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xpos_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd l xi0) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (0 <= proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n)%Re.

Lemma lower_even_Xneg_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  Z.Even (Z_of_nat n) ->
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xneg_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd l xi0) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n <= 0)%Re.

Lemma lower_odd_Xpos_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  Z.Odd (Z_of_nat n) ->
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xpos_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd l xi0) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n <= 0)%Re.

Lemma lower_odd_Xneg_over
  (X : I.type) (l := Xlower (I.convert X)) (u := Xupper (I.convert X))
  (c r0 : R) (xi0 := Xreal r0) (nm1 : nat) (n := nm1.+1) :
  Z.Odd (Z_of_nat n) ->
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  Xneg_over (XDn n.+1) l u ->
  forall x : R, contains (IIbnd l xi0) (Xreal x) ->
  contains (IIbnd (Xreal x) xi0) (Xreal c) \/
    contains (IIbnd xi0 (Xreal x)) (Xreal c) ->
  (0 <= proj_val (XDn n.+1 (Xreal c)) / INR (fact n) * (x - r0) ^ n)%Re.

Lemma Xsub_0_r (x : ExtendedR) : Xsub x (Xreal 0) = x.

Lemma Rdiv_1_r (x : R) : (x / 1 = x)%Re.

Lemma Xdiv_1_r (x : ExtendedR) : (x / Xreal 1 = x)%XR.

Proposition 2.2.1 in Mioara Joldes' thesis, adapted from Lemma 5.12 in Roland Zumkeller's thesis
Theorem Zumkeller_monot_rem (X : I.type) (r0 : R) (xi0 := Xreal r0) (n : nat)
  (l := I.convert_bound (I.lower X)) (u := I.convert_bound (I.upper X)) :
  I.bounded X = true ->
  contains (I.convert X) xi0 ->
  ~ Xnan_ex XDn n.+1 (I.convert X) ->
  Xcst_sign (XDn n.+1) l u ->
  Xmonot (Xdelta n xi0) l xi0 /\
  Xmonot (Xdelta n xi0) xi0 u.

Lemma Ztech_derive_sign (X : I.type) (n : nat) :
  (exists t, contains (I.convert X) t) ->
  ~ Xnan_ex XDn n.+1 (I.convert X) ->
  I.bounded X = true ->
  isNNegOrNPos (tnth (IP X n.+1) n.+1) = true ->
  Xcst_sign (XDn n.+1) (I.convert_bound (I.lower X)) (I.convert_bound (I.upper X)).

Theorem i_validTM_Ztech X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (IP X0 n) (Ztech prec IP (IP X0 n) F0 X0 X n)) XF0.

End GenericProof.

Section rec1_correct.

Variable (XF_rec : FullXR.T -> FullXR.T -> nat -> FullXR.T).
Variable (F_rec : T -> T -> nat -> T).
Class extension2N := Extension2N :
  forall (ix iy : I.type) (x y : ExtendedR) (n : nat),
  contains (I.convert ix) x ->
  contains (I.convert iy) y ->
  contains (I.convert (F_rec ix iy n)) (XF_rec x y n).
Context { H_F_rec : extension2N }.

Class compat_rec1 := Compat_rec1 :
  forall (r : FullXR.T) (k : nat),
  XF_rec r (XDn k r / Xreal (INR (fact k)))%XR k.+1 =
  (XDn k.+1 r / Xreal (INR (fact k.+1)))%XR.
Context { H_XF_rec : compat_rec1 }.

Definition IP_rec1 X0 := trec1 (F_rec X0) (F0 X0).
Definition XP_rec1 xi0 := PolX.trec1 (XF_rec xi0) (XF0 xi0).

Instance validXPoly_rec1 :
  validXPoly XP_rec1.

Instance validPoly_rec1 :
  validPoly XP_rec1 IP_rec1.

Let's apply the generic proof (for testing purposes).

Corollary TM_rec1_correct X0 X n :
  let err := TLrem prec (fun t => trec1 (F_rec X) (F0 t)) X0 X n in
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (trec1 (F_rec X0) (F0 X0) n) err) XF0.

Corollary TM_rec1_correct' X0 X n :
  let err := Ztech prec (fun t => trec1 (F_rec X) (F0 t))
    (trec1 (F_rec X0) (F0 X0) n) F0 X0 X n in
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (trec1 (F_rec X0) (F0 X0) n) err) XF0.

End rec1_correct.

Section rec2_correct.

Variable (XF1 : FullXR.T -> FullXR.T). Variable (F1 : T -> T).
Hypothesis XDn_1 : forall x, XDn 1 x = XF1 x.
Hypothesis F1_contains : I.extension XF1 F1.

Variable (XF_rec : FullXR.T -> FullXR.T -> FullXR.T -> nat -> FullXR.T).
Variable (F_rec : T -> T -> T -> nat -> T).

Class extension3N := Extension3N :
  forall (ix iy iz : I.type) (x y z : ExtendedR) (n : nat),
  contains (I.convert ix) x ->
  contains (I.convert iy) y ->
  contains (I.convert iz) z ->
  contains (I.convert (F_rec ix iy iz n)) (XF_rec x y z n).
Context { H_F_rec : extension3N }.

Class compat_rec2 := Compat_rec2 :
  forall (r : FullXR.T) (k : nat),
  XF_rec r
  (XDn k r / Xreal (INR (fact k)))%XR
  (XDn k.+1 r / Xreal (INR (fact k.+1)))%XR
  k.+2 =
  (XDn k.+2 r / Xreal (INR (fact k.+2)))%XR.
Context { H_XF_rec : compat_rec2 }.

Definition IP_rec2 X0 := trec2 (F_rec X0) (F0 X0) (F1 X0).
Definition XP_rec2 xi0 := PolX.trec2 (XF_rec xi0) (XF0 xi0) (XF1 xi0).

Instance validXPoly_rec2 :
  validXPoly XP_rec2.

Instance validPoly_rec2 :
  validPoly XP_rec2 IP_rec2.

Let's apply the generic proof (for testing purposes).

Corollary TM_rec2_correct X0 X n :
  let err := TLrem prec
    (fun t => trec2 (F_rec X) (F0 t) (F1 t)) X0 X n in
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (trec2 (F_rec X0) (F0 X0) (F1 X0) n) err) XF0.

Corollary TM_rec2_correct' X0 X n :
  let err := Ztech prec (fun t i => trec2 (F_rec X) (F0 t) (F1 t) i)
    (trec2 (F_rec X0) (F0 X0) (F1 X0) n) F0 X0 X n in
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X)
  (RPA (trec2 (F_rec X0) (F0 X0) (F1 X0) n) err) XF0.

End rec2_correct.

End ProofOfRec.

Theorem TM_cst_correct n (icst X0 X : I.type) (cst : ExtendedR) :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  contains (I.convert icst) cst ->
  i_validTM (I.convert X0) (I.convert X) (TM_cst icst X0 X n) (Xmask cst).

Lemma TM_var_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_var X0 X n) (fun x => x).

Lemma TM_inv_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_inv X0 X n) (FullXR.tinv tt).

Lemma TM_exp_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_exp X0 X n) (FullXR.texp tt).

Lemma TM_sqrt_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_sqrt X0 X n) (FullXR.tsqrt tt).

Ltac Inc :=
  rewrite /tnat /FullXR.tnat INR_IZR_INZ -Z2R_IZR;
  apply: I.fromZ_correct.

Lemma TM_invsqrt_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_invsqrt X0 X n)(FullXR.tinvsqrt tt).

Lemma TM_sin_correct X X0 n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_sin X0 X n) (FullXR.tsin tt).

Lemma TM_cos_correct X0 X n :
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) (TM_cos X0 X n) (FullXR.tcos tt).

Local Notation "a + b" := (FullXR.tadd tt a b).
Local Notation "a - b" := (FullXR.tsub tt a b).

Lemma Xneg_Xadd (a b : ExtendedR) : Xneg (Xadd a b) = Xadd (Xneg a) (Xneg b).

Lemma TM_add_correct_gen (smallX0 : interval) (X : I.type) (TMf TMg : rpa) f g :
  I.subset_ smallX0 (I.convert X) ->
  Pol.tsize (approx TMf) = Pol.tsize (approx TMg) ->
  i_validTM smallX0 (I.convert X) TMf f ->
  i_validTM smallX0 (I.convert X) TMg g ->
  i_validTM smallX0 (I.convert X) (TM_add TMf TMg)
  (fun xr => FullXR.tadd tt (f xr) (g xr)).

Lemma TM_add_correct (X0 X : I.type) (TMf TMg : rpa) f g :
  Pol.tsize (approx TMf) = Pol.tsize (approx TMg) ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  i_validTM (I.convert X0) (I.convert X) TMg g ->
  i_validTM (I.convert X0) (I.convert X) (TM_add TMf TMg)
  (fun xr => FullXR.tadd tt (f xr) (g xr)).

Definition mul_error prec n (f g : rpa) X0 X :=
 let pf := approx f in
 let pg := approx g in
 let sx := (I.sub prec X X0) in
 let B := I.mul prec (Pol.teval prec (Pol.tmul_tail prec n pf pg) sx)
                (I.power_int prec sx (Z_of_nat n.+1)) in
 let Bf := Pol.teval prec pf sx in
 let Bg := Pol.teval prec pg sx in
   I.add prec B (I.add prec (I.mul prec (error f) Bg)
     (I.add prec (I.mul prec (error g) Bf)
       (I.mul prec (error f) (error g)))).

Definition TM_mul (Mf Mg : rpa) X0 X n : rpa :=
 RPA (Pol.tmul_trunc prec n (approx Mf) (approx Mg))
     (mul_error prec n Mf Mg X0 X).

Lemma teval_poly_nan k p x :
  k < PolX.tsize p -> PolX.tnth p k = Xnan ->
  PolX.teval tt p x = Xnan.

Lemma poly_mul n (a b : nat -> R) (x : R) :
  (forall i, n < i -> (a i = R0)/\ (b i = R0))->
  \big[Rplus/R0]_(k < (n + n).+1)
    \big[Rplus/R0]_(i < k.+1) (a i * b (k - i)%N * x ^ k)%Re =
  \big[Rplus/R0]_(i < n.+1)
    \big[Rplus/R0]_(j < n.+1) (a i * b j * x ^ (i + j))%Re.

Lemma TM_mul_correct_gen (smallX0: interval) (X0 X: I.type) (TMf TMg: rpa) f g :
  I.subset_ smallX0 (I.convert X0) ->
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains smallX0 t) ->
  let n := Pol.tsize (approx TMf) in
  n = Pol.tsize (approx TMg)-> 0 < n ->
  i_validTM smallX0 (I.convert X) TMf f ->
  i_validTM smallX0 (I.convert X) TMg g ->
  i_validTM smallX0 (I.convert X) (TM_mul TMf TMg X0 X n.-1)
  (fun xr => FullXR.tmul tt (f xr) (g xr)).

Lemma TM_mul_correct (X0 X : I.type) (TMf TMg : rpa) f g :
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  let n := Pol.tsize (approx TMf) in
  n = Pol.tsize (approx TMg) -> 0 < n ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  i_validTM (I.convert X0) (I.convert X) TMg g ->
  i_validTM (I.convert X0) (I.convert X) (TM_mul TMf TMg X0 X n.-1)
  (fun xr => FullXR.tmul tt (f xr) (g xr)).

Definition poly_eval_tm n p (Mf : rpa) (X0 X : I.type) : rpa :=
  @Pol.tfold rpa
  (fun a b => (TM_add (TM_cst a X0 X n) (TM_mul b Mf X0 X n)))
  (TM_cst (I.fromZ 0) X0 X n) p.

Lemma size_TM_add Mf Mg :
  tsize (approx (TM_add Mf Mg)) =
  maxn (tsize (approx Mf)) (tsize (approx Mg)).

Lemma size_TM_mul Mf Mg n X0 X :
  tsize (approx (TM_mul Mf Mg X0 X n)) = n.+1.

Lemma size_poly_eval_tm n p Mf X0 X :
  tsize (approx (poly_eval_tm n p Mf X0 X)) = n.+1.

Lemma TM_fun_eq X0 X TMf f g :
  (forall x, f x = g x) ->
  i_validTM X0 X TMf f -> i_validTM X0 X TMf g.


Lemma poly_eval_tm_correct_aux_gen smallX0 X0 X Mf f pi :
  I.subset_ smallX0 (I.convert X0) ->
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains smallX0 t) ->
  f Xnan = Xnan ->
  i_validTM smallX0 (I.convert X) Mf f ->
  let n := tsize pi in
  forall nf, tsize (approx Mf) = nf.+1 ->
  forall pr, n = PolX.tsize pr ->
  (forall k, k < n -> contains (I.convert (tnth pi k)) (PolX.tnth pr k)) ->
  i_validTM smallX0 (I.convert X) (poly_eval_tm nf pi Mf X0 X)
  (fun x => @PolX.tfold _
    (fun a b => FullXR.tadd tt a (FullXR.tmul tt b (f x)))
    (Xmask FullXR.tzero x) pr).

Lemma poly_eval_tm_correct_gen (smallX0 : interval) (X0 X : I.type) Mf f p :
  I.subset_ smallX0 (I.convert X0) ->
  I.subset_ (I.convert X0) (I.convert X) ->
  (exists t : ExtendedR, contains smallX0 t) ->
  f Xnan = Xnan ->
  i_validTM smallX0 (I.convert X) Mf f ->
  forall n, tsize p = n.+1 ->
  forall nf, tsize (approx Mf) = nf.+1 ->
  forall pr, tsize p = PolX.tsize pr ->
  (forall k, k < tsize p ->
  contains (I.convert (tnth p k)) (PolX.tnth pr k)) ->
  i_validTM smallX0 (I.convert X)
  (poly_eval_tm nf p Mf X0 X)
  (fun x => PolX.teval tt pr (f x)).

Lemma poly_eval_tm_correct_aux X0 X Mf f p :
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  f Xnan = Xnan ->
  i_validTM (I.convert X0) (I.convert X) Mf f ->
  let n := tsize p in
  forall nf, tsize (approx Mf) = nf.+1 ->
  forall pr, n = PolX.tsize pr ->
  (forall k, k < n -> contains (I.convert (tnth p k)) (PolX.tnth pr k)) ->
  i_validTM (I.convert X0) (I.convert X) (poly_eval_tm nf p Mf X0 X)
  (fun x => @PolX.tfold _
    (fun a b => FullXR.tadd tt a (FullXR.tmul tt b (f x)))
    (Xmask FullXR.tzero x) pr).

Lemma poly_eval_tm_correct X0 X Mf f p :
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  f Xnan = Xnan ->
  i_validTM (I.convert X0) (I.convert X) Mf f ->
  forall n, tsize p = n.+1 ->
  forall nf, tsize (approx Mf) = nf.+1 ->
  forall pr, PolX.tsize pr = tsize p ->
  (forall k, k < tsize p ->
  contains (I.convert (tnth p k)) (PolX.tnth pr k)) ->
  i_validTM (I.convert X0) (I.convert X) (poly_eval_tm nf p Mf X0 X)
  (fun x => PolX.teval tt pr (f x)).

Definition TM_type := I.type -> I.type -> nat -> rpa.

Definition TMset0 (Mf : rpa) t :=
  RPA (Pol.tset_nth (approx Mf) 0 t) (error Mf).

Definition TM_comp (TMg : TM_type) (Mf : rpa) X0 X n :=
  let Bf := Pol.teval prec (approx Mf) (I.sub prec X X0) in
  let Mg := TMg (Pol.tnth (approx Mf) 0) (I.add prec Bf (error Mf)) n in
  let M1 := TMset0 Mf (I.fromZ 0) in
  let M0 := poly_eval_tm n (approx Mg) M1 X0 X in
  RPA (approx M0) (I.add prec (error M0) (error Mg)).

REMARK: the TM is supposed not void

Lemma TMset0_correct X0 X TMf f :
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  forall nf, tsize (approx TMf) = nf.+1 ->
  forall fi0, contains (I.convert X0) fi0 ->
  exists a0, i_validTM (Interval_interval.Ibnd fi0 fi0) (I.convert X)
  (TMset0 TMf (I.fromZ 0)) (fun x => f x - a0).


Lemma TM_comp_correct (X0 X : I.type) (Tyg : TM_type) (TMf : rpa) g f :
  f Xnan = Xnan ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  (exists t : ExtendedR, contains (I.convert (tnth (approx TMf) 0)) t) ->
  forall n, Pol.tsize (approx TMf) = n.+1 ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  (forall Y0 Y k, I.subset_ (I.convert Y0) (I.convert Y) ->
    (exists t : ExtendedR, contains (I.convert Y0) t) ->
    i_validTM (I.convert Y0) (I.convert Y) (Tyg Y0 Y k) g
    /\ tsize (approx (Tyg Y0 Y k)) = k.+1) ->
  i_validTM (I.convert X0) (I.convert X)
  (TM_comp Tyg TMf X0 X n) (fun xr => (g (f xr))).

Definition TM_inv_comp Mf X0 X (n : nat) := TM_comp TM_inv Mf X0 X n.

Lemma TM_inv_comp_correct (X0 X : I.type) (TMf : rpa) f :
  f Xnan = Xnan ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  (exists t : ExtendedR, contains (I.convert (tnth (approx TMf) 0)) t) ->
  forall n, Pol.tsize (approx TMf) = n.+1 ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  i_validTM (I.convert X0) (I.convert X)
  (TM_inv_comp TMf X0 X n) (fun xr => FullXR.tinv tt (f xr)).

Definition TM_div Mf Mg X0 X n :=
   TM_mul Mf (TM_inv_comp Mg X0 X n) X0 X n.

Lemma TM_div_correct (X0 X : I.type) (TMf TMg : rpa) f g :
  g Xnan = Xnan->
  (exists t : ExtendedR, contains (I.convert (tnth (approx TMg) 0)) t) ->
  (exists t : ExtendedR, contains (I.convert X0) t) ->
  forall n, Pol.tsize (approx TMf) = n.+1 ->
  n.+1 = Pol.tsize (approx TMg) ->
  i_validTM (I.convert X0) (I.convert X) TMf f ->
  i_validTM (I.convert X0) (I.convert X) TMg g ->
  i_validTM (I.convert X0) (I.convert X)
  (TM_div TMf TMg X0 X n) (fun xr => FullXR.tdiv tt (f xr) (g xr)).

End PrecArgument.
End TaylorModel.