Library Interval.Poly.Taylor_model_sharp

From Coq Require Import ZArith Psatz Reals.
From Flocq Require Import Raux.
From Coquelicot Require Import Coquelicot.
From mathcomp.ssreflect Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq fintype bigop.

Require Import Stdlib.
Require Import MathComp.
Require Import Coquelicot.
Require Import Xreal.
Require Import Interval.
Require Import Basic.
Require Import Xreal_derive.
Require Import Interval_compl.
Require Import Datatypes.
Require Import Taylor_poly.
Require Import Basic_rec.
Require Import Bound.

Set Implicit Arguments.

Local Open Scope nat_scope.

Record rpa {pol int : Type} : Type := RPA { approx : pol ; error : int }.

Module TR := TaylorPoly FullR PolR.

Module TaylorModel (I : IntervalOps) (Pol : PolyIntOps I) (Bnd : PolyBound I Pol).
Import Pol.Notations.
Import PolR.Notations.
Local Open Scope ipoly_scope.
Module Export Aux := IntervalAux I.
Module Import TI := TaylorPoly Pol.Int Pol.
Module Import BndThm := PolyBoundThm I Pol Bnd.
Module J := IntervalExt I.


Ltac step_xr xr :=
  match goal with
  [ |- contains ?i ?g ] ⇒ rewrite -(_ : xr = g)
  end.
Ltac step_r r :=
  match goal with
  [ |- contains ?i (Xreal ?g) ] ⇒ rewrite -(_ : r = g)
  end.
Ltac step_xi xi :=
  match goal with
  [ |- contains ?g ?xr ] ⇒ rewrite -(_ : xi = g)
  end.
Ltac step_i i :=
  match goal with
  [ |- contains (I.convert ?g) ?xr ] ⇒ rewrite -(_ : i = g)
  end.


Notation rpa := (@rpa Pol.T I.type).

Section PrecArgument.

Variable prec : I.precision.

Section TaylorModel.

Variable M : rpa.

Variable Tcoeffs : I.type → nat → Pol.T.

Definition TLrem (X0 X : I.type) n :=
  let N := S n in
  let NthCoeff := Pol.nth (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.

Definition Ztech (P : Pol.T) F (X0 X : I.type) n :=
  let N := S n in
  let NthCoeff := Pol.nth (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) (Pol.horner prec P (I.sub prec A X0)) in
    let Db := I.sub prec (F B) (Pol.horner prec P (I.sub prec B X0)) in
    let Dx0 := I.sub prec (F X0) (Pol.nth P 0) 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.

End TaylorModel.

Definition TM_cst (X c : I.type) : rpa :=
  RPA (Pol.polyC c) (I.mask (I.mask I.zero X) c).

Definition TM_var (X X0 : I.type) :=
  RPA (Pol.set_nth Pol.polyX 0 X0) (I.mask I.zero X).

Definition TM_exp (X0 X : I.type) (n : nat) : rpa :=
  

  let P := (T_exp prec X0 n) in
  RPA P (Ztech (T_exp prec) P (I.exp prec) X0 X n).

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

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

Definition TM_atan (X0 X : I.type) (n : nat) : rpa :=
  let P := (T_atan prec X0 n) in
  RPA P (Ztech (T_atan prec) P (I.atan prec) X0 X n).

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

Definition TM_opp (M : rpa) : rpa :=
  RPA (Pol.opp (approx M)) (I.neg (error M)).

Definition TM_sub (Mf Mg : rpa) : rpa :=
  RPA (Pol.sub prec (approx Mf) (approx Mg))
      (I.sub prec (error Mf) (error Mg)).

Definition i_validTM (x0 : R) (X : interval)
  (M : rpa) (xf : R → ExtendedR) :=
  [/\
    ∀ x : R, contains X (Xreal x) → xf x = Xnan → I.convert (error M) = IInan,
    X = IInan → I.convert (error M) = IInan,
    contains (I.convert (error M)) (Xreal 0),
    contains X (Xreal x0) &
    exists2 Q, approx M >:: Q
    & ∀ x, contains X (Xreal x) →
      error M >: proj_val (xf x) - Q.[x - x0]].

Lemma TM_fun_eq f g (x0 : R) (X : interval) TMf :
  (∀ x, contains X (Xreal x) → f x = g x) →
  i_validTM x0 X TMf f → i_validTM x0 X TMf g.

Section TM_integral.

Local Open Scope R_scope.

Lemma Rpol_continuous p (x : R) : continuous (PolR.horner tt p) x.

Lemma Rpol_integrable p (x1 x2 : R) : ex_RInt (PolR.horner tt p) x1 x2.


Lemma ex_derive_big I (s : seq I) :
   ∀
   (f : I → R → R)
   (x : R),
   (∀ k : I, ex_derive (f k) x) →
 ex_derive (fun x0 : Rdefinitions.R ⇒ \big[Rplus/0]_(i <- s) f i x0) x.

Lemma derive_big I (s : seq I) :
   ∀
   (f : I → R → R)
   (x : R),
   (∀ k : I, ex_derive (f k) x) →
   Derive (fun y : Rdefinitions.R ⇒
             \big[Rplus/0]_(i <- s) (f i y)) x =
   \big[(fun _ : unit ⇒ Rplus) tt/0]_(i <- s) ((Derive (f i)) x).

Lemma horner_primitive (c : R) (p : PolR.T) (t : R):
  PolR.horner tt (PolR.primitive tt c p) t =
  c + \big[Rplus/0]_(0 ≤ i < (size p))
   (PolR.nth (PolR.primitive tt c p) i.+1 × powerRZ t (Z.of_nat i.+1)).

Lemma Rpol_derive p (c : R) (x : R) :
  Derive (PolR.horner tt (PolR.primitive tt c p)) x = PolR.horner tt p x.

Lemma Rpol_integral_0 (x1 x2 : R) (p : PolR.T) :
  RInt (PolR.horner tt p) x1 x2 =
    PolR.horner tt (PolR.primitive tt 0 p) x2 - PolR.horner tt (PolR.primitive tt 0 p) x1.

Local Notation Isub := (I.sub prec).
Local Notation Imul := (I.mul prec).
Local Notation Iadd := (I.add prec).

Variables (X0 X : I.type).
Variable xF : R → ExtendedR. Let f := fun x ⇒ proj_val (xF x).
Let iX0 := I.convert X0.
Let iX := I.convert X.

Hypothesis f_int : ∀ x x0 : R,
  contains iX (Xreal x) → contains iX0 (Xreal x0) →
  ex_RInt f x0 x.
Hypothesis x_Def : ∀ x : R, contains iX (Xreal x) → xF x ≠ Xnan.
Variable Mf : rpa.


Definition TM_integral_poly :=
  Pol.primitive prec (I.zero) (approx Mf).

Definition TM_integral_error :=
  Imul (Isub X X0) (error Mf).

Local Open Scope R_scope.

Lemma pol_int_sub pol x1 x2 x3 :
  ex_RInt (fun y : R ⇒ pol.[y - x3]) x1 x2.

Section NumericIntegration.

Local Open Scope R_scope.

Variables (x0 a b : R) (ia ib : I.type).

Hypothesis Hx0 : contains iX0 (Xreal x0).
Hypothesis Ha : contains iX (Xreal a).
Hypothesis Hb : contains iX (Xreal b).
Hypothesis Hia : contains (I.convert ia) (Xreal a).
Hypothesis Hib : contains (I.convert ib) (Xreal b).
Hypothesis f_int_numeric : ex_RInt f a b.

Definition polIntegral := Isub (Pol.horner prec TM_integral_poly (Isub ib X0)) (Pol.horner prec TM_integral_poly (Isub ia X0)).
Definition integralError := Imul (Isub ib ia) (error Mf).
Definition integralEnclosure := Iadd polIntegral integralError.

Lemma integralEnclosure_correct :
  i_validTM x0 iX Mf xF →
  contains (I.convert integralEnclosure) (Xreal (RInt f a b)).

End NumericIntegration.

Lemma contains_interval_float_integral (p : PolR.T) :
  approx Mf >:: p →
  TM_integral_poly >:: (PolR.primitive tt 0%R p).

Lemma TM_integral_error_0 (x0 : R) :
  contains iX0 (Xreal x0) →
  i_validTM x0 iX Mf xF →
  contains (I.convert TM_integral_error) (Xreal 0).

Definition TM_integral :=
  RPA TM_integral_poly TM_integral_error.

Lemma TM_integral_correct (x0 : Rdefinitions.R) :
  contains iX0 (Xreal x0) →
  i_validTM x0 iX Mf xF →
  i_validTM x0 iX TM_integral (fun x ⇒ Xreal (RInt f x0 x)).

End TM_integral.

Section Const_prelim.

Definition is_const (f : R → ExtendedR) (X c : interval) : Prop :=
  exists2 y : ExtendedR, contains c y
  & ∀ x : R, contains X (Xreal x) → f x = y.

Lemma is_const_ext (f g : R → ExtendedR) (X c : interval) :
  (∀ x : R, contains X (Xreal x) → f x = g x) →
  is_const f X c → is_const g X c.

Corollary is_const_ext_weak (f g : R → ExtendedR) (X c : interval) :
  (∀ x : R, f x = g x) →
  is_const f X c → is_const g X c.

End Const_prelim.

Section GenericProof.

Variable xf : R → ExtendedR.
Variable P : R → nat → PolR.T.

Let f0 := fun x ⇒ proj_val (xf x).
Let Dn n := Derive_n f0 n.

Let Rdelta (n : nat) (x0 x : R) :=
  (f0 x - (P x0 n).[x - x0])%R.

Let Rdelta' (n : nat) (x0 x : R) :=
  (Dn 1 x - (PolR.deriv tt (P x0 n)).[x - x0])%R.

Section aux.

Variable dom : R → Prop.
Hypothesis Hdom : connected dom.

Lemma Rmonot_contains (g : R → R) :
  Rmonot dom g →
  ∀ (x y z : R),
  dom x → dom y → intvl x y z →
  intvl (g x) (g y) (g z) ∨ intvl (g y) (g x) (g z).

Lemma upper_Rpos_over
  (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  dom x0 →
  Rpos_over dom (Dn n.+1) →
  ∀ x : R, (x0 ≤ x)%R → dom x → intvl x x0 c ∨ intvl x0 x c →
  (0 ≤ (Dn n.+1 c) / INR (fact n) × (x - x0) ^ n)%R.

Lemma upper_Rneg_over
  (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  dom x0 →
  Rneg_over dom (Dn n.+1) →
  ∀ x : R, (x0 ≤ x)%R → dom x → intvl x x0 c ∨ intvl x0 x c →
  (Dn n.+1 c / INR (fact n) × (x - x0) ^ n ≤ 0)%R.

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)%R.

Lemma powerRZ_1_even (k : Z) : (0 ≤ powerRZ (-1) (2 × k))%R.

Lemma ZEven_pow_le (r : R) (n : nat) :
  Z.Even (Z_of_nat n) →
  (0 ≤ r ^ n)%R.

Lemma Ropp_le_0 (x : R) :
  (0 ≤ x → - x ≤ 0)%R.

Lemma ZOdd_pow_le (r : R) (n : nat) :
  Z.Odd (Z_of_nat n) →
  (r ≤ 0 → r ^ n ≤ 0)%R.

Lemma lower_even_Rpos_over
  (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  Z.Even (Z_of_nat n) →
  dom x0 →
  Rpos_over dom (Dn n.+1) →
  ∀ x : R, (x ≤ x0)%R → dom x → intvl x x0 c ∨ intvl x0 x c →
  (0 ≤ Dn n.+1 c / INR (fact n) × (x - x0) ^ n)%R.

Lemma lower_even_Rneg_over
  (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  Z.Even (Z_of_nat n) →
  dom x0 →
  Rneg_over dom (Dn n.+1) →
  ∀ x : R, (x ≤ x0)%R → dom x → intvl x x0 c ∨ intvl x0 x c →
  (Dn n.+1 c / INR (fact n) × (x - x0) ^ n ≤ 0)%R.

Lemma lower_odd_Rpos_over
  (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  Z.Odd (Z_of_nat n) →
  dom x0 →
  Rpos_over dom (Dn n.+1) →
  ∀ x : R, (x ≤ x0)%R → dom x → intvl x x0 c ∨ intvl x0 x c →
  (Dn n.+1 c / INR (fact n) × (x - x0) ^ n ≤ 0)%R.

Lemma lower_odd_Rneg_over (c x0 : R) (nm1 : nat) (n := nm1.+1) :
  Z.Odd (Z_of_nat n) →
  dom x0 →
  Rneg_over dom (Dn n.+1) →
  ∀ x : R, dom x → (x ≤ x0)%R → intvl x x0 c ∨ intvl x0 x c →
  (0 ≤ (Dn n.+1 c) / INR (fact n) × (x - x0) ^ n)%R.

Hypothesis Hder_n : ∀ n x, dom x → ex_derive_n f0 n x.

Lemma Rderive_delta (Pr : R → Prop) (n : nat) (x0 : R) :
  dom x0 →
  Pr x0 →
  Rderive_over (fun t ⇒ dom t ∧ Pr t) (Rdelta n x0) (Rdelta' n x0).

Hypothesis Poly_size : ∀ (x0 : R) n, PolR.size (P x0 n) = n.+1.
Hypothesis Poly_nth : ∀ (x : R) n k, dom x → k ≤ n →
  PolR.nth (P x n) k = Rdiv (Dn k x) (INR (fact k)).

Lemma bigXadd'_P (m n : nat) (x0 s : R) :
  dom x0 →
  ex_derive_n f0 n x0 →
  m ≤ n →
  \big[Rplus/R0]_(0 ≤ i < m) (PolR.nth (P x0 n) i.+1 × INR i.+1 × s ^ i)%R =
  \big[Rplus/R0]_(0 ≤ i < m) ((Dn i.+1 x0) / INR (fact i) × s ^ i)%R.

Theorem Zumkeller_monot_rem (x0 : R) (n : nat) :
  dom x0 →
  Rcst_sign dom (Dn n.+1) →
  Rmonot (fun t ⇒ dom t ∧ (t ≤ x0)%R) (Rdelta n x0) ∧
  Rmonot (fun t ⇒ dom t ∧ (x0 ≤ t)%R) (Rdelta n x0).

End aux.

Variable F : I.type → I.type.
Variable IP : I.type → nat → Pol.T.

Hypothesis F_contains : I.extension (Xbind xf) F.

Variables X : I.type.

Class validPoly : Prop := ValidPoly {
  Poly_size : ∀ (x0 : R) n, PolR.size (P x0 n) = n.+1;
  Poly_nth :
    ∀ (x : R) n k,
    X >: x →
    k ≤ n →
    PolR.nth (P x n) k = Rdiv (Dn k x) (INR (fact k)) }.

Class validIPoly : Prop := ValidIPoly {
  IPoly_size :
    ∀ (X0 : I.type) x0 n, eq_size (IP X0 n) (P x0 n);
  IPoly_nth : ∀ (X0 : I.type) x0 n, X0 >: x0 → IP X0 n >:: P x0 n;
  IPoly_nai :
    ∀ X, ∀ r : R, contains (I.convert X) (Xreal r) → xf r = Xnan →
    ∀ n k, k ≤ n → I.convert (Pol.nth (IP X n) k) = IInan
}.

Context { validPoly_ : validPoly }.
Context { validIPoly_ : validIPoly }.

Hypothesis Hder_n : ∀ n x, X >: x → ex_derive_n f0 n x.

Lemma Poly_nth0 x n : X >: x → PolR.nth (P x n) 0 = f0 x.

Theorem i_validTM_TLrem (x0 : R) (X0 : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (RPA (IP X0 n) (TLrem IP X0 X n)) xf.

Lemma Ztech_derive_sign (n : nat) :
  isNNegOrNPos (Pol.nth (IP X n.+1) n.+1) = true →
  Rcst_sign (fun t ⇒ contains (I.convert X) (Xreal t)) (Dn n.+1).

Lemma F_Rcontains : ∀ X x, X >: x → F X >: f0 x.

Theorem i_validTM_Ztech (x0 : R) (X0 : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X)
  (RPA (IP X0 n) (Ztech IP (IP X0 n) F X0 X n)) xf.

End GenericProof.

Lemma size_TM_cst (X c : I.type) : Pol.size (approx (TM_cst X c)) = 1.

Theorem TM_cst_correct (x0 : R) (ci X : I.type) (c : ExtendedR) :
  contains (I.convert X) (Xreal x0) →
  contains (I.convert ci) c →
  i_validTM x0 (I.convert X) (TM_cst X ci) (fun _ ⇒ c).

Theorem TM_cst_correct_strong (x0 : R) (ci X : I.type) (f : R → ExtendedR) :
  contains (I.convert X) (Xreal x0) →
  is_const f (I.convert X) (I.convert ci) →
  i_validTM x0 (I.convert X) (TM_cst X ci) f.

Definition TM_any (Y : I.type) (X : I.type) (n : nat) :=
  let mid := J.midpoint Y in
  let pol := Pol.polyC mid in
  {| approx := if n == 0 then pol
               else Pol.set_nth pol n Pol.Int.zero;
     error := I.mask (I.sub prec Y mid) X
  |}.

Definition sizes := (Pol.size_polyNil, Pol.size_polyCons,
                     PolR.size_polyNil, PolR.size_polyCons,
                     Pol.size_set_nth, PolR.size_set_nth,
                     Pol.polyCE).

Lemma size_TM_any (c X : I.type) n : Pol.size (approx (TM_any c X n)) = n.+1.

Theorem TM_any_correct
  (x0 : R) (Y X : I.type) (n : nat) (f : R → ExtendedR) :
  contains (I.convert X) (Xreal x0) →
  (∀ x : R, contains (I.convert X) (Xreal x) →
    contains (I.convert Y) (f x)) →
  i_validTM x0 (I.convert X) (TM_any Y X n) f.

Lemma size_TM_var (X X0 : I.type) : Pol.size (approx (TM_var X X0)) = 2.

Lemma TM_var_correct (x0 : R) (X0 X : I.type) :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_var X X0) Xreal.

Theorem TM_var_correct_strong (x0 : R) (X0 X : I.type) (f : R → ExtendedR) :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  (∀ x : R, contains (I.convert X) (Xreal x) → f x = Xreal x) →
  i_validTM x0 (I.convert X) (TM_var X X0) f.

Lemma size_TM_exp (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_exp X0 X n)) = n.+1.

Lemma TM_exp_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_exp X0 X n) (fun x ⇒ Xreal (exp x)).

Lemma nat_ind2 (P : nat → Prop) :
  P 0 → P 1 →
  (∀ k, P k → P k.+1 → P k.+2) →
  ∀ k, P k.

Lemma size_TM_sin (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_sin X0 X n)) = n.+1.

Lemma TM_sin_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_sin X0 X n) (fun x ⇒ Xreal (sin x)).

Lemma size_TM_cos (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_cos X0 X n)) = n.+1.

Lemma TM_cos_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_cos X0 X n) (fun x ⇒ Xreal (cos x)).

Lemma is_derive_n_atan n (x : R) :
  let q n := iteri n
    (fun i c ⇒ PolR.div_mixed_r tt (PolR.sub tt (PolR.add tt c^`() (PolR.lift 2 c^`()))
      (PolR.mul_mixed tt (INR (i.+1).*2) (PolR.lift 1 c))) (INR i.+2))
    PolR.one in
  is_derive_n atan n x
  (if n is n.+1 then PolR.horner tt (q n) x / (1 + x²) ^ (n.+1) × INR (fact n.+1)
   else atan x)%R.

Lemma size_TM_atan (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_atan X0 X n)) = n.+1.

Lemma TM_atan_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_atan X0 X n) (fun x ⇒ Xreal (atan x)).

Definition TM_tan (X0 X : I.type) (n : nat) : rpa :=
  let P := (T_tan prec X0 n) in
  let ic := I.cos prec X in
  if apart0 ic
  then RPA P (Ztech (T_tan prec) P (I.tan prec) X0 X n)
  else RPA P I.nai.

Lemma is_derive_n_tan n x :
  cos x ≠ 0%R →
  let q n := iteri n
    (fun i c ⇒ PolR.div_mixed_r tt (PolR.add tt c^`() (PolR.lift 2 c^`())) (INR i.+1))
    (PolR.lift 1 PolR.one) in
  is_derive_n tan n x ((q n).[tan x] × INR (fact n)).

Lemma size_TM_tan (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_tan X0 X n)) = n.+1.

Lemma TM_tan_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_tan X0 X n) Xtan'.


Definition TM_sqrt (X0 X : I.type) (n : nat) : rpa :=
  
  let P := (T_sqrt prec X0 n) in
  if gt0 X
  then RPA P (Ztech (T_sqrt prec) P (I.sqrt prec) X0 X n)
  else RPA P I.nai.

Lemma size_TM_sqrt (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_sqrt X0 X n)) = n.+1.

Lemma TM_sqrt_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_sqrt X0 X n) Xsqrt'.

Definition I_invsqrt prec x := I.inv prec (I.sqrt prec x).

Definition TM_invsqrt (X0 X : I.type) (n : nat) : rpa :=
  
  let P := (T_invsqrt prec X0 n) in
  if gt0 X
  then RPA P (Ztech (T_invsqrt prec) P (I_invsqrt prec) X0 X n)
  else RPA P I.nai.

Lemma size_TM_invsqrt (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_invsqrt X0 X n)) = n.+1.

Ltac Inc :=
  rewrite INR_IZR_INZ;
  apply: I.fromZ_correct.

Lemma TM_invsqrt_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_invsqrt X0 X n) (fun x ⇒ Xinv (Xsqrt (Xreal x))).

Definition TM_power_int (p : Z) (X0 X : I.type) (n : nat) :=
  let P := (T_power_int prec p X0 n) in
  if p is Z.neg _ then
    if apart0 X then
      RPA P (Ztech (T_power_int prec p) P
                   (fun x ⇒ I.power_int prec x p) X0 X n)
    else RPA P I.nai
  else RPA P (Ztech (T_power_int prec p) P
                    (fun x ⇒ I.power_int prec x p) X0 X n).

Lemma size_TM_power_int (p : Z) (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_power_int p X0 X n)) = n.+1.

Lemma toR_power_int p x : (0 ≤ p)%Z ∨ x ≠ R0 →
  powerRZ x p = proj_val (Xpower_int' x p).

Lemma toR_power_int_loc p x : (0 ≤ p)%Z ∨ x ≠ R0 →
  locally x (fun t ⇒ powerRZ t p = proj_val (Xpower_int' t p)).

Lemma TM_power_int_correct_aux (p : Z) (x0 : R) (X0 X : I.type) n :
  (0 ≤ p)%Z ∨ apart0 X →
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X)
    (let P := (T_power_int prec p X0 n) in
      RPA P (Ztech (T_power_int prec p) P (fun x ⇒ I.power_int prec x p) X0 X n))
    (fun x ⇒ Xpower_int' x p).

Lemma TM_power_int_correct (p : Z) (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_power_int p X0 X n) (fun x ⇒ Xpower_int' x p).

Definition TM_inv (X0 X : I.type) (n : nat) :=
  let P := (T_inv prec X0 n) in
  if apart0 X then
    RPA P (Ztech (T_inv prec) P (I.inv prec) X0 X n)
  else RPA P I.nai.

Lemma size_TM_inv (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_inv X0 X n)) = n.+1.

Lemma TM_inv_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_inv X0 X n) Xinv'.

Definition TM_ln (X0 X : I.type) (n : nat) : rpa :=
  let P := (T_ln prec X0 n) in
  if gt0 X then
    RPA P (Ztech (T_ln prec) P (I.ln prec) X0 X n)
  else RPA P I.nai.

Lemma size_TM_ln (X0 X : I.type) (n : nat) :
  Pol.size (approx (TM_ln X0 X n)) = n.+1.

Lemma powerRZ_opp x n :
  x ≠ 0%R → powerRZ x (- n) = / (powerRZ x n).

Lemma TM_ln_correct (x0 : R) (X0 X : I.type) n :
  subset' (I.convert X0) (I.convert X) →
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) (TM_ln X0 X n) Xln'.

Local Notation "a + b" := (Xadd a b).
Local Notation "a - b" := (Xsub a b).

Lemma TM_add_correct (x0 : R) (X : interval) (TMf TMg : rpa) f g :
  i_validTM x0 X TMf f →
  i_validTM x0 X TMg g →
  i_validTM x0 X (TM_add TMf TMg) (fun xr ⇒ Xadd (f xr) (g xr)).

Lemma TM_opp_correct (x0 : R) (X : interval) (TMf : rpa) f :
  i_validTM x0 X TMf f →
  i_validTM x0 X (TM_opp TMf) (fun xr ⇒ Xneg (f xr)).

Lemma TM_sub_correct (x0 : R) (X : interval) (TMf TMg : rpa) f g :
  i_validTM x0 X TMf f →
  i_validTM x0 X TMg g →
  i_validTM x0 X (TM_sub TMf TMg) (fun xr ⇒ Xsub (f xr) (g xr)).

Definition TM_mul_mixed (a : I.type) (M : rpa) : rpa :=
  RPA (Pol.map (I.mul prec a) (approx M))
      (I.mul prec a (error M)).

Definition TM_div_mixed_r (M : rpa) (b : I.type) : rpa :=
  RPA (Pol.map (I.div prec ^~ b) (approx M))
      (I.div prec (error M) b).

Lemma size_TM_mul_mixed (a : I.type) M :
  Pol.size (approx (TM_mul_mixed a M)) = Pol.size (approx M).

Lemma size_TM_div_mixed_r M (b : I.type) :
  Pol.size (approx (TM_div_mixed_r M b)) = Pol.size (approx M).

Lemma TM_mul_mixed_correct (a : I.type) M (x0 : R) (X : interval) f (y : R) :
  a >: y →
  i_validTM x0 X M f →
  i_validTM x0 X (TM_mul_mixed a M) (fun x ⇒ Xmul (Xreal y) (f x)).

Lemma TM_mul_mixed_nai M (x0 : R) (X : interval) a f g :
  I.convert a = Inan →
  i_validTM x0 X M f →
  i_validTM x0 X (TM_mul_mixed a M) g.

Lemma TM_mul_mixed_correct_strong (a : I.type) M (x0 : R) (X : interval) f g :
  is_const f X (I.convert a) →
  i_validTM x0 X M g →
  i_validTM x0 X (TM_mul_mixed a M) (fun x ⇒ Xmul (f x) (g x)).

Lemma TM_div_mixed_r_aux0 M (b : I.type) (x0 : R) (X : interval) f :
  b >: 0%R →
  i_validTM x0 X M f →
  i_validTM x0 X (TM_div_mixed_r M b) (fun x ⇒ Xdiv (f x) (Xreal 0)).

Lemma TM_div_mixed_r_correct M (b : I.type) (x0 : R) (X : interval) f (y : R) :
  b >: y →
  i_validTM x0 X M f →
  i_validTM x0 X (TM_div_mixed_r M b) (fun x ⇒ Xdiv (f x) (Xreal y)).

Lemma TM_div_mixed_r_nai M (x0 : R) (X : interval) a f g :
  I.convert a = Inan →
  i_validTM x0 X M f →
  i_validTM x0 X (TM_div_mixed_r M a) g.

Lemma TM_div_mixed_r_correct_strong M (b : I.type) (x0 : R) (X : interval) f g :
  i_validTM x0 X M f →
  is_const g X (I.convert b) →
  i_validTM x0 X (TM_div_mixed_r M b) (fun x ⇒ Xdiv (f x) (g x)).

Definition mul_error prec n (f g : rpa) (X0 X : I.type) :=
 let pf := approx f in
 let pg := approx g in
 let sx := (I.sub prec X X0) in
 let B := I.mul prec (Bnd.ComputeBound prec (Pol.mul_tail prec n pf pg) sx)
                (I.power_int prec sx (Z_of_nat n.+1)) in
 let Bf := Bnd.ComputeBound prec pf sx in
 let Bg := Bnd.ComputeBound 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 : I.type) n : rpa :=
 RPA (Pol.mul_trunc prec n (approx Mf) (approx Mg))
     (mul_error prec n Mf Mg X0 X).

Lemma TM_mul_correct (x0 : R) (X0 X : I.type) (TMf TMg : rpa) f g n :
  contains (I.convert X0) (Xreal x0) →
  i_validTM x0 (I.convert X) TMf f →
  i_validTM x0 (I.convert X) TMg g →
  i_validTM x0 (I.convert X) (TM_mul TMf TMg X0 X n) (fun xr ⇒ Xmul (f xr) (g xr)).

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

Lemma size_TM_mul Mf Mg n (X0 X : I.type) :
  Pol.size (approx (TM_mul Mf Mg X0 X n)) = n.+1.

Lemma size_TM_sub Mf Mg :
  Pol.size (approx (TM_sub Mf Mg)) =
  maxn (Pol.size (approx Mf)) (Pol.size (approx Mg)).

Lemma size_TM_opp Mf :
  Pol.size (approx (TM_opp Mf)) = Pol.size (approx Mf).

Definition TM_horner n p (Mf : rpa) (X0 X : I.type) : rpa :=
  @Pol.fold rpa
  (fun a b ⇒ (TM_add (TM_cst X a) (TM_mul b Mf X0 X n)))
  (TM_cst X I.zero) p.

Lemma size_TM_horner n p Mf (X0 X : I.type) :
  Pol.size (approx (TM_horner n p Mf X0 X)) = (if 0 < Pol.size p then n else 0).+1.