Library Interval.Interval.Interval

From Coq Require Import Bool Reals Psatz.
From Coquelicot Require Import Coquelicot.

Require Import Stdlib.
Require Import Xreal.
Require Import Basic.

Inductive interval : Set :=
  | Inan : interval
  | Ibnd (l u : ExtendedR) : interval.

Definition Xlower (xi : interval) : ExtendedR :=
  match xi with
  | Ibnd xl _ ⇒ xl
  | _ ⇒ Xnan
  end.

Definition Xupper (xi : interval) : ExtendedR :=
  match xi with
  | Ibnd _ xu ⇒ xu
  | _ ⇒ Xnan
  end.

Definition contains i v :=
  match i, v with
  | Ibnd l u, Xreal x ⇒
    match l with
    | Xreal r ⇒ Rle r x
    | Xnan ⇒ True
    end ∧
    match u with
    | Xreal r ⇒ Rle x r
    | Xnan ⇒ True
    end
  | Inan, _ ⇒ True
  | _, Xnan ⇒ False
  end.

Inductive output_bound : Set :=
  | BInteger : Z → output_bound
  | BDecimal : QArith_base.Q → output_bound
  | BFraction : Z → Z → output_bound.

Import Stdlib.Compatibility Rdefinitions.

Definition convert_bound b :=
  match b with
  | BInteger n ⇒ IZR n
  | BDecimal q ⇒ Q2R q
  | BFraction n d ⇒ (IZR n / IZR d)%R
  end.

Definition contains_output xi x :=
  match xi with
  | (None, None) ⇒ True
  | (None, Some xu) ⇒ (x ≤ convert_bound xu)%R
  | (Some xl, None) ⇒ (convert_bound xl ≤ x)%R
  | (Some xl, Some xu) ⇒ (convert_bound xl ≤ x ≤ convert_bound xu)%R
  end.

Lemma contains_connected :
  ∀ xi, connected (fun x ⇒ contains xi (Xreal x)).

Lemma contains_Xnan :
  ∀ xi, contains xi Xnan ↔ xi = Inan.

Lemma contains_Inan :
  ∀ xi x, xi = Inan → contains xi x.

Lemma contains_le :
  ∀ xl xu v,
  contains (Ibnd xl xu) v →
  le_lower xl v ∧ le_upper v xu.

Lemma le_contains :
  ∀ xl xu v,
  le_lower xl (Xreal v) → le_upper (Xreal v) xu →
  contains (Ibnd xl xu) (Xreal v).

Definition subset xi yi :=
  match xi, yi with
  | Ibnd xl xu, Ibnd yl yu ⇒
    match xl, xu with
    | Xreal xl, Xreal xu ⇒ (xu < xl)%R
    | _, _ ⇒ False
    end
    ∨ (le_lower yl xl ∧ le_upper xu yu)
  | _, Inan ⇒ True
  | _, _ ⇒ False
  end.

Definition subset' xi yi :=
  ∀ x, contains xi x → contains yi x.

Theorem subset_contains :
  ∀ xi yi,
  subset xi yi →
  subset' xi yi.

Definition domain' P b :=
  ∀ x, contains b (Xreal x) → P x.

Theorem bisect' :
  ∀ P xl xm xu,
  domain' P (Ibnd xl xm) →
  domain' P (Ibnd xm xu) →
  domain' P (Ibnd xl xu).

Definition not_empty xi :=
  ∃ v, contains xi (Xreal v).

Lemma contains_Xreal :
  ∀ xi x,
  contains xi x →
  contains xi (Xreal (proj_val x)).

Lemma not_empty_contains :
  ∀ xi x,
  contains xi x →
  not_empty xi.

Lemma contains_lower :
  ∀ l u x,
  contains (Ibnd (Xreal l) u) x →
  contains (Ibnd (Xreal l) u) (Xreal l).

Lemma contains_upper :
  ∀ l u x,
  contains (Ibnd l (Xreal u)) x →
  contains (Ibnd l (Xreal u)) (Xreal u).

Module Type IntervalBounds.

Parameter type : Type.
Parameter nan : type.
Parameter convert : type → ExtendedR.
Parameter precision : Type.
Parameter PtoP : positive → precision.

End IntervalBounds.

Module Type IntervalBasicOps.

Declare Module F : IntervalBounds.

Parameter type : Type.
Parameter valid_lb : F.type → Prop.
Parameter valid_ub : F.type → Prop.
Parameter convert : type → interval.
Parameter zero : type.
Parameter nai : type.
Parameter empty : type.
Parameter bnd : F.type → F.type → type.
Parameter singleton : F.type → type.
Parameter real : type → bool.
Parameter is_empty : type → bool.

Parameter valid_lb_real :
  ∀ b, F.convert b = Xreal (proj_val (F.convert b)) → valid_lb b.

Parameter valid_ub_real :
  ∀ b, F.convert b = Xreal (proj_val (F.convert b)) → valid_ub b.

Parameter valid_lb_nan : valid_lb F.nan.

Parameter valid_ub_nan : valid_ub F.nan.

Parameter bnd_correct :
  ∀ l u, valid_lb l → valid_ub u →
  convert (bnd l u) = Ibnd (F.convert l) (F.convert u).

Parameter singleton_correct :
  ∀ b,
  contains (convert (singleton b)) (Xreal (proj_val (F.convert b))).

Parameter zero_correct :
  convert zero = Ibnd (Xreal 0) (Xreal 0).

Parameter nai_correct :
  convert nai = Inan.

Parameter empty_correct :
  ∀ x, contains (convert empty) x → False.

Parameter real_correct :
  ∀ xi, real xi = match convert xi with Inan ⇒ false | _ ⇒ true end.

Parameter is_empty_correct :
  ∀ xi x, contains (convert xi) x →
  is_empty xi = true → False.

Parameter output : bool → type → option output_bound × option output_bound.

Parameter output_correct :
  ∀ fmt xi x, contains (convert xi) (Xreal x) → contains_output (output fmt xi) x.

Parameter subset : type → type → bool.

Parameter subset_correct :
  ∀ xi yi v,
  contains (convert xi) v →
  subset xi yi = true →
  contains (convert yi) v.

Parameter join : type → type → type.
Parameter meet : type → type → type.
Parameter sign_large : type → Xcomparison.
Parameter sign_strict : type → Xcomparison.

Parameter sign_large_correct :
  ∀ xi,
  match sign_large xi with
  | Xeq ⇒ ∀ x, contains (convert xi) x → x = Xreal 0
  | Xlt ⇒ ∀ x, contains (convert xi) x → x = Xreal (proj_val x) ∧ Rle (proj_val x) 0
  | Xgt ⇒ ∀ x, contains (convert xi) x → x = Xreal (proj_val x) ∧ Rle 0 (proj_val x)
  | Xund ⇒ True
  end.

Parameter sign_strict_correct :
  ∀ xi,
  match sign_strict xi with
  | Xeq ⇒ ∀ x, contains (convert xi) x → x = Xreal 0
  | Xlt ⇒ ∀ x, contains (convert xi) x → x = Xreal (proj_val x) ∧ Rlt (proj_val x) 0
  | Xgt ⇒ ∀ x, contains (convert xi) x → x = Xreal (proj_val x) ∧ Rlt 0 (proj_val x)
  | Xund ⇒ True
  end.

Parameter join_correct :
  ∀ xi yi v,
  contains (convert xi) v ∨ contains (convert yi) v →
  contains (convert (join xi yi)) v.

Parameter meet_correct :
  ∀ xi yi v,
  contains (convert xi) v → contains (convert yi) v →
  contains (convert (meet xi yi)) v.

Parameter meet_correct' :
  ∀ xi yi v,
  contains (convert (meet xi yi)) v →
  contains (convert xi) v ∧ contains (convert yi) v.

Parameter midpoint : type → F.type.

Parameter midpoint_correct :
  ∀ xi,
  not_empty (convert xi) →
  F.convert (midpoint xi) = Xreal (proj_val (F.convert (midpoint xi))) ∧
  contains (convert xi) (F.convert (midpoint xi)).

Parameter bisect : type → type × type.

Parameter bisect_correct :
  ∀ xi x,
  contains (convert xi) x →
  contains (convert (fst (bisect xi))) x ∨ contains (convert (snd (bisect xi))) x.

Definition extension f fi := ∀ b x,
  contains (convert b) x → contains (convert (fi b)) (f x).

Definition extension_2 f fi := ∀ ix iy x y,
  contains (convert ix) x →
  contains (convert iy) y →
  contains (convert (fi ix iy)) (f x y).

Parameter mask : type → type → type.

Parameter mask_correct : extension_2 Xmask mask.
Parameter mask_correct' :
  ∀ xi yi x,
  contains (convert xi) x →
  contains (convert (mask xi yi)) x.

Definition precision := F.precision.

Parameter wider : precision → type → type → bool.

Parameter neg : type → type.
Parameter abs : type → type.
Parameter inv : precision → type → type.
Parameter invnz : precision → type → type.
Parameter sqr : precision → type → type.
Parameter sqrt : precision → type → type.
Parameter add : precision → type → type → type.
Parameter sub : precision → type → type → type.
Parameter mul : precision → type → type → type.
Parameter div : precision → type → type → type.
Parameter power_int : precision → type → Z → type.
Parameter nearbyint : rounding_mode → type → type.
Parameter error_fix : precision → rounding_mode → Z → type → type.
Parameter error_flt : precision → rounding_mode → Z → positive → type → type.

Parameter neg_correct : extension Xneg neg.
Parameter abs_correct : extension Xabs abs.
Parameter inv_correct : ∀ prec, extension Xinv (inv prec).
Parameter sqr_correct : ∀ prec, extension Xsqr (sqr prec).
Parameter sqrt_correct : ∀ prec, extension Xsqrt (sqrt prec).
Parameter add_correct : ∀ prec, extension_2 Xadd (add prec).
Parameter sub_correct : ∀ prec, extension_2 Xsub (sub prec).
Parameter mul_correct : ∀ prec, extension_2 Xmul (mul prec).
Parameter div_correct : ∀ prec, extension_2 Xdiv (div prec).
Parameter power_int_correct : ∀ prec n, extension (fun x ⇒ Xpower_int x n) (fun x ⇒ power_int prec x n).
Parameter nearbyint_correct : ∀ mode, extension (Xnearbyint mode) (nearbyint mode).
Parameter error_fix_correct : ∀ prec mode emin, extension (Xerror_fix mode emin) (error_fix prec mode emin).
Parameter error_flt_correct : ∀ prec mode emin p, extension (Xerror_flt mode emin p) (error_flt prec mode emin p).

Parameter neg_correct' :
  ∀ xi x,
  contains (convert (neg xi)) (Xneg x) →
  contains (convert xi) x.

Parameter invnz_correct :
  ∀ prec xi x,
  x ≠ Xreal 0 → contains (convert xi) x → contains (convert (invnz prec xi)) (Xinv x).

Parameter cancel_add : precision → type → type → type.
Parameter cancel_sub : precision → type → type → type.

Parameter bounded : type → bool.
Parameter lower_bounded : type → bool.
Parameter upper_bounded : type → bool.

Parameter lower_extent : type → type.
Parameter upper_extent : type → type.
Parameter whole : type.

Parameter lower_extent_correct :
  ∀ xi x y,
  contains (convert xi) (Xreal y) →
  (x ≤ y)%R →
  contains (convert (lower_extent xi)) (Xreal x).

Parameter upper_extent_correct :
  ∀ xi x y,
  contains (convert xi) (Xreal y) →
  (y ≤ x)%R →
  contains (convert (upper_extent xi)) (Xreal x).

Parameter whole_correct :
  ∀ x,
  contains (convert whole) (Xreal x).

Parameter lower_complement : type → type.
Parameter upper_complement : type → type.

Parameter lower_complement_correct :
  ∀ xi x y,
  contains (convert xi) (Xreal x) →
  contains (convert (lower_complement xi)) (Xreal y) →
  (y ≤ x)%R.

Parameter upper_complement_correct :
  ∀ xi x y,
  contains (convert xi) (Xreal x) →
  contains (convert (upper_complement xi)) (Xreal y) →
  (x ≤ y)%R.

Parameter lower : type → F.type.
Parameter upper : type → F.type.

Parameter lower_correct :
  ∀ xi : type,
  not_empty (convert xi) →
  F.convert (lower xi) = Xlower (convert xi).

Parameter valid_lb_lower :
  ∀ xi : type,
  not_empty (convert xi) →
  valid_lb (lower xi).

Parameter upper_correct :
  ∀ xi : type,
  not_empty (convert xi) →
  F.convert (upper xi) = Xupper (convert xi).

Parameter valid_ub_upper :
  ∀ xi : type,
  not_empty (convert xi) →
  valid_ub (upper xi).

Definition bounded_prop xi :=
  not_empty (convert xi) →
  convert xi = Ibnd (F.convert (lower xi)) (F.convert (upper xi)).

Parameter lower_bounded_correct :
  ∀ xi,
  lower_bounded xi = true →
  F.convert (lower xi) = Xreal (proj_val (F.convert (lower xi))) ∧
  bounded_prop xi.

Parameter upper_bounded_correct :
  ∀ xi,
  upper_bounded xi = true →
  F.convert (upper xi) = Xreal (proj_val (F.convert (upper xi))) ∧
  bounded_prop xi.

Parameter bounded_correct :
  ∀ xi,
  bounded xi = true →
  lower_bounded xi = true ∧ upper_bounded xi = true.

Parameter fromZ_small : Z → type.

Parameter fromZ_small_correct :
  ∀ v,
  (Z.abs v ≤ 256)%Z →
  contains (convert (fromZ_small v)) (Xreal (IZR v)).

Parameter fromZ : precision → Z → type.

Parameter fromZ_correct :
  ∀ prec v, contains (convert (fromZ prec v)) (Xreal (IZR v)).

Definition propagate_l fi :=
  ∀ xi yi : type, convert xi = Inan → convert (fi xi yi) = Inan.
Definition propagate_r fi :=
  ∀ xi yi : type, convert yi = Inan → convert (fi xi yi) = Inan.

Parameter mask_propagate_l : propagate_l mask.
Parameter mask_propagate_r : propagate_r mask.
Parameter add_propagate_l : ∀ prec, propagate_l (add prec).
Parameter sub_propagate_l : ∀ prec, propagate_l (sub prec).
Parameter mul_propagate_l : ∀ prec, propagate_l (mul prec).
Parameter div_propagate_l : ∀ prec, propagate_l (div prec).
Parameter add_propagate_r : ∀ prec, propagate_r (add prec).
Parameter sub_propagate_r : ∀ prec, propagate_r (sub prec).
Parameter mul_propagate_r : ∀ prec, propagate_r (mul prec).
Parameter div_propagate_r : ∀ prec, propagate_r (div prec).

End IntervalBasicOps.

Module Type IntervalOps.

Include IntervalBasicOps.

Parameter pi : precision → type.
Parameter cos : precision → type → type.
Parameter sin : precision → type → type.
Parameter tan : precision → type → type.
Parameter atan : precision → type → type.
Parameter exp : precision → type → type.
Parameter ln : precision → type → type.

Parameter pi_correct : ∀ prec, contains (convert (pi prec)) (Xreal PI).
Parameter cos_correct : ∀ prec, extension Xcos (cos prec).
Parameter sin_correct : ∀ prec, extension Xsin (sin prec).
Parameter tan_correct : ∀ prec, extension Xtan (tan prec).
Parameter atan_correct : ∀ prec, extension Xatan (atan prec).
Parameter exp_correct : ∀ prec, extension Xexp (exp prec).
Parameter ln_correct : ∀ prec, extension Xln (ln prec).

End IntervalOps.

Module IntervalBasicExt (I : IntervalBasicOps).

Lemma nai_correct :
  ∀ x, contains (I.convert I.nai) x.

Lemma contains_le :
  ∀ xi x,
  contains (I.convert xi) x →
  le_lower (I.F.convert (I.lower xi)) x ∧ le_upper x (I.F.convert (I.upper xi)).

Definition propagate fi :=
  ∀ xi, I.convert xi = Inan → I.convert (fi xi) = Inan.

Lemma propagate_extension :
  ∀ fi f, I.extension (Xbind f) fi → propagate fi.

Lemma neg_propagate : propagate I.neg.

Lemma inv_propagate : ∀ prec, propagate (I.inv prec).

Lemma sqrt_propagate : ∀ prec, propagate (I.sqrt prec).

Lemma power_int_propagate : ∀ prec n, propagate (fun x ⇒ I.power_int prec x n).

Lemma error_fix_propagate : ∀ prec mode emin, propagate (I.error_fix prec mode emin).

Lemma error_flt_propagate : ∀ prec mode emin p, propagate (I.error_flt prec mode emin p).

Definition extension f fi :=
  ∀ (xi : I.type) (x : R),
  contains (I.convert xi) (Xreal x) →
  contains (I.convert (fi xi)) (Xreal (f x)).

Definition extension_2 f fi :=
  ∀ (xi yi : I.type) (x y : R),
  contains (I.convert xi) (Xreal x) →
  contains (I.convert yi) (Xreal y) →
  contains (I.convert (fi xi yi)) (Xreal (f x y)).

Lemma neg_correct : extension Ropp I.neg.

Lemma abs_correct : extension Rabs I.abs.

Lemma inv_correct : ∀ prec, extension Rinv (I.inv prec).

Lemma invnz_correct :
  ∀ prec xi x,
  x ≠ 0%R →
  contains (I.convert xi) (Xreal x) →
  contains (I.convert (I.invnz prec xi)) (Xreal (/ x)).

Lemma sqr_correct : ∀ prec, extension Rsqr (I.sqr prec).

Lemma sqrt_correct : ∀ prec, extension sqrt (I.sqrt prec).

Lemma add_correct : ∀ prec, extension_2 Rplus (I.add prec).

Lemma sub_correct : ∀ prec, extension_2 Rminus (I.sub prec).

Lemma mul_correct : ∀ prec, extension_2 Rmult (I.mul prec).

Lemma div_correct : ∀ prec, extension_2 Rdiv (I.div prec).

Lemma nearbyint_correct :
  ∀ mode, extension (Rnearbyint mode) (I.nearbyint mode).

Definition round_fix prec mode emin xi :=
  I.add prec xi (I.error_fix prec mode emin xi).

Lemma round_fix_correct :
  ∀ mode prec emin, extension (Basic.round_fix mode emin) (round_fix prec mode emin).

Lemma round_fix_correct' :
  ∀ mode prec emin, I.extension (Xround_fix mode emin) (round_fix prec mode emin).

Definition round_flt prec mode emin p xi :=
  I.add prec xi (I.error_flt prec mode emin p xi).

Lemma round_flt_correct :
  ∀ mode prec emin p, extension (Basic.round_flt mode emin p) (round_flt prec mode emin p).

Lemma round_flt_correct' :
  ∀ mode prec emin p, I.extension (Xround_flt mode emin p) (round_flt prec mode emin p).

Lemma power_int_correct :
  ∀ prec n, extension (fun x ⇒ powerRZ x n) (fun xi ⇒ I.power_int prec xi n).

Lemma zero_correct : contains (I.convert I.zero) (Xreal 0).

Lemma contains_only_0 r :
  contains (I.convert I.zero) (Xreal r) → r = 0%R.

Lemma join_correct :
  ∀ ui vi u v x,
  contains (I.convert ui) (Xreal u) →
  contains (I.convert vi) (Xreal v) →
  (u ≤ x ≤ v)%R →
  contains (I.convert (I.join ui vi)) (Xreal x).

Lemma contains_RInt prec (f3 : R → R) x1 x2 Y X1 X2 :
  ex_RInt f3 x1 x2→
  contains (I.convert X1) (Xreal x1) →
  contains (I.convert X2) (Xreal x2) →
  (∀ x, (Rmin x1 x2 ≤ x ≤ Rmax x1 x2)%R → contains (I.convert Y) (Xreal (f3 x))) →
  contains (I.convert (I.mul prec (I.sub prec X2 X1) Y)) (Xreal (RInt f3 x1 x2)).

Definition midpoint xi :=
  let m := I.midpoint xi in I.bnd m m.

Lemma midpoint_correct :
  ∀ xi, not_empty (I.convert xi) →
  I.convert (midpoint xi) = let m := Xreal (proj_val (I.F.convert (I.midpoint xi))) in Ibnd m m.

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

Lemma subset_midpoint :
  ∀ xi, not_empty (I.convert xi) →
  subset' (I.convert (midpoint xi)) (I.convert xi).

End IntervalBasicExt.

Module IntervalExt (I : IntervalOps).

Include (IntervalBasicExt I).

Lemma cos_correct : ∀ prec, extension cos (I.cos prec).

Lemma sin_correct : ∀ prec, extension sin (I.sin prec).

Lemma tan_correct : ∀ prec, extension tan (I.tan prec).

Lemma atan_correct : ∀ prec, extension atan (I.atan prec).

Lemma exp_correct : ∀ prec, extension exp (I.exp prec).

Lemma ln_correct : ∀ prec, extension ln (I.ln prec).

End IntervalExt.