Library Interval.Integral.Bertrand

From Coq Require Import Reals ZArith Psatz Fourier_util.
From Flocq Require Import Raux.
From Coquelicot Require Import Coquelicot AutoDerive.
From mathcomp.ssreflect Require Import ssreflect ssrfun ssrbool ssrnat bigop.

Require Import Stdlib.
Require Import Coquelicot.
Require Import Xreal.
Require Import Interval.

Section powerRZMissing.

Lemma powerRZ_ind (P : Z (R R) Prop) :
  P 0%Z (fun x ⇒ 1)
  ( n, P (Z.of_nat n) (fun xx ^ n))
  ( n, P ((-(Z.of_nat n))%Z) (fun x/ x ^ n))
  ( f g, f =1 g n, P n f P n g)
   (m : Z), P m (fun xpowerRZ x m).

Lemma is_derive_powerRZ (n : Z) (x : R):
  ((0 n)%Z x 0)
  is_derive (fun x : RpowerRZ x n) x (IZR n × (powerRZ x (n - 1))).

Lemma is_derive_powerRZS (n : Z) (x : R):
  ((1 n)%Z x 0)
  is_derive (fun x : RpowerRZ x (n+1)) x (IZR (n+1) × (powerRZ x n)).

Lemma ex_derive_powerRZ (n : Z) (x : R):
  ((0 n)%Z x 0)
  ex_derive (fun x : RpowerRZ x n) x.

Lemma ex_derive_powerRZS (n : Z) (x : R):
  ((1 n)%Z x 0)
  ex_derive (fun x : RpowerRZ x (n+1)) x.

Lemma is_RInt_powerRZ (alpha : Z) (a b : R) (HneqN1 : alpha (-1)%Z) (H : 0 < a b) :
is_RInt (powerRZ^~ alpha) a b
    ((powerRZ b (alpha + 1) - powerRZ a (alpha + 1)) / IZR (alpha + 1)).

Lemma int_x_alpha alpha A B (H : 0 < A B) (Halpha: alpha (-1)%Z) :
is_RInt (powerRZ^~ alpha) A B
((powerRZ B (alpha + 1) - powerRZ A (alpha + 1)) / IZR (alpha + 1)).

End powerRZMissing.

Section CoquelicotMissing.

Lemma continuous_Rdiv_1_x x (H : x 0) : continuous (Rdiv 1) x.

End CoquelicotMissing.

Definition Bertrand alpha beta A B (I : R) :=
  is_RInt (fun xpowerRZ x alpha × (pow (ln x) beta)) A B I.


Fixpoint f (alpha : Z) (beta : nat) (A B : R) {struct beta} :=
  match beta with
    | 0 ⇒ (powerRZ B (alpha+1)- powerRZ A (alpha+1)) / (IZR (alpha + 1))
    | S m
      (powerRZ B (alpha+1) × (pow (ln B) (beta)) -
       powerRZ A (alpha+1) × (pow (ln A) beta)) / (IZR (alpha + 1)) -
      (INR beta) / (IZR (alpha+1)) × f alpha m A B end.

Definition Bertrand_lim alpha beta (A : R) (I : R) :=
  is_RInt_gen (fun xpowerRZ x alpha × (pow (ln x) beta)) (at_point A) (Rbar_locally p_infty) I.

Fixpoint f_lim (alpha : Z) (beta : nat) (A : R) {struct beta} :=
  match beta with
    | 0 ⇒ (- powerRZ A (alpha+1)) / (IZR (alpha + 1))
    | S m
       - ( powerRZ A (alpha+1) × (pow (ln A) beta)) / (IZR (alpha + 1)) -
      (INR beta) / (IZR (alpha+1)) × f_lim alpha m A end.



Lemma one_step_by_parts alpha beta (A : R) (B : R) (H : 0 < A B) (Halpha: alpha (-1)%Z) :
   I, Bertrand alpha beta A B I
  Bertrand alpha (S beta) A B
     ((powerRZ B (alpha+1) × (pow (ln B) (S beta)) -
       powerRZ A (alpha+1) × (pow (ln A) (S beta))) / (IZR (alpha + 1)) -
      (INR (S beta)) / (IZR (alpha+1)) × I).

Lemma f_correct alpha beta A B (H : 0 < A B) (Halpha: alpha (-1)%Z) :
 Bertrand alpha beta A B (f alpha beta A B).

Lemma f_correct_RInt alpha beta A B (H : 0 < A B) (Halpha: alpha (-1)%Z) : f alpha beta A B = RInt (fun xpowerRZ x alpha × (pow (ln x) beta)) A B.

Lemma is_lim_RInv_p_infty:
is_lim [eta Rinv] p_infty 0.

Lemma is_lim_RInv_m_infty:
is_lim [eta Rinv] m_infty 0.

Lemma is_lim_powerRZ_0 alpha (Halpha : (alpha < 0)%Z) :
  is_lim (powerRZ^~ (alpha)%Z) p_infty (0%R).

Lemma is_lim_pow_infty n : is_lim (fun xx^n.+1) p_infty p_infty.

Lemma is_lim_pow_0 (f : R R) n :
  is_lim f p_infty 0
  is_lim (fun x(f x)^n.+1) p_infty 0.

Lemma Rbar_mult_p_l y (Hy : 0 < y) :
  Rbar_mult y p_infty = p_infty.

Lemma Rbar_mult_p_r y (Hy : 0 < y) :
  Rbar_mult p_infty y = p_infty.

Lemma Rbar_mult_m_l y (Hy : 0 < y) :
  Rbar_mult y m_infty = m_infty.

Lemma Rbar_mult_m_r y (Hy : 0 < y) :
  Rbar_mult m_infty y = m_infty.

Lemma is_lim_Rpower y (Hy : 0 < y) :
  is_lim (fun xRpower x y) p_infty p_infty.

Lemma x_alpha_0 alpha (Halpha : (alpha < -1)%Z) :
  is_lim (powerRZ^~ (alpha + 1)%Z) p_infty (0%R).

Lemma Rpower_pos {x y} (Hx : 0 < x) (Hy : 0 y) : Rpower x y > 0.

Lemma is_derive_Rpower {x y} (Hx : 0 < x) (Hy : 0 y) :
  is_derive (fun tRpower t y) x (y × Rpower x (y - 1)).

Lemma ln_Rpower x y (Hx : 0 < x) (Hy : 0 y) : ln (Rpower x y) = y × ln x.

Lemma x_alpha_beta alpha beta (Halpha : (alpha < -1)%Z) :
  is_lim (fun xpowerRZ x (alpha + 1)%Z × (pow (ln x) beta.+1)) p_infty (0%R).

Lemma f_lim_is_lim alpha beta A (H : 0 < A) (Halpha : (alpha < -1)%Z):
   filterlim (f alpha beta A) (Rbar_locally p_infty)
   (locally (f_lim alpha beta A)).

Lemma f_lim_correct alpha beta A (H : 0 < A) (Halpha : (alpha < -1)%Z) :
 Bertrand_lim alpha beta A (f_lim alpha beta A).

Section BertrandLogNeg.


Definition f_neg (a : R) (beta : nat) :=
  - / ((INR beta) × (ln a)^beta).

Lemma continuous_f_neg x beta (H0x : 0 < x) (Hx1 : x 1):
  continuous (fun x0 : R/ (x0 × ln x0 ^ beta)) x.

Lemma is_derive_f_neg :
   beta x,
  beta 0%N
  0 < x x 1
  is_derive (fun x : Rf_neg x beta) x (/ (x × (ln x)^beta.+1)).

Lemma f_neg_correct_RInt_0_1 a b beta (Hab1 : 0 < a b) (Hb1 : b < 1) (Hbeta : beta 0%N) : is_RInt (fun x/ (x × (ln x)^beta.+1)) a b (f_neg b beta - f_neg a beta).

Lemma f_neg_correct_RInt_a_infty a b beta (Ha : 1 < a b) (Hbeta : beta 0%N) : is_RInt (fun x/ (x × (ln x)^beta.+1)) a b (f_neg b beta - f_neg a beta).

Lemma filterlim_sqr_m_infty:
  filterlim (pow^~ 2%N) (Rbar_locally m_infty) (Rbar_locally p_infty).

Lemma is_lim_sqr_infty:
  is_lim (pow^~ 2%N) m_infty p_infty.

Lemma filterlim_pow_infty n : filterlim (pow^~ n.+1) (Rbar_locally p_infty) (Rbar_locally p_infty).

Lemma filterlim_pow_m_even n : ~~ odd n.+1 filterlim (fun xpow x n.+1) (Rbar_locally m_infty) (Rbar_locally p_infty).

Lemma filterlim_pow_m_odd n : odd n.+1 filterlim (fun xpow x n.+1) (Rbar_locally m_infty) (Rbar_locally m_infty).


Lemma f_neg_correct_RInt_gen_0_a a beta (Ha : 0 < a < 1) (Hbeta : beta 0%N) : is_RInt_gen (fun x/ (x × (ln x)^beta.+1)) (at_right 0) (at_point a) (f_neg a beta).

Lemma f_neg_correct_RInt_gen_a_infty a beta (Ha : 1 < a) (Hbeta : beta 0%N) : is_RInt_gen (fun x/ (x × (ln x)^beta.+1)) (at_point a) (Rbar_locally p_infty) (- f_neg a beta).

End BertrandLogNeg.

Section ExponentInQ.

End ExponentInQ.

Section ZeroToEpsilon.


Definition f0eps (alpha : Z) (beta : nat) (epsilon : R) (B : R) :=
  (-1) ^ beta × f (- 2 - alpha) beta (/ epsilon) B.

Definition f0eps_lim (alpha : Z) (beta : nat) (epsilon : R) :=
  (-1) ^ beta × f_lim (- 2 - alpha) beta (/ epsilon).

Lemma pow_negx x n : pow (- x) n = (pow (-1) n) × pow x n.

Lemma subst_lemma alpha beta epsilon (eta : R) (Heps : 0 < epsilon) (Heta : 0 < eta epsilon) (Halpha : -1 < IZR alpha) :
  RInt_gen
    (fun xpowerRZ x alpha × pow (ln x) beta)
    (at_point eta)
    (at_point epsilon) =
  - RInt
      (fun x- (pow (-1) beta) × powerRZ x (- 2 - alpha) × (pow (ln x) beta))
      (1 / epsilon)
      (1 / eta).

Lemma f0eps_correct alpha beta epsilon (B : R) (Heps : 0 < / B epsilon) (HB : 0 < B) (Halpha : (-1 < alpha)%Z) :
  is_RInt_gen ((fun xpowerRZ x alpha × (pow (ln x) beta))) (at_point (/ B)) (at_point epsilon) (f0eps alpha beta epsilon B).

Lemma f0eps_correct_sing alpha beta epsilon sing (B : R) (Heps : 0 < / B epsilon) (HB : 0 < B) (Halpha : (-1 < alpha)%Z) :
  is_RInt_gen ((fun xpowerRZ (x - sing) alpha × (pow (ln (x - sing)) beta))) (at_point (sing + / B)) (at_point (sing + epsilon)) (f0eps alpha beta epsilon B).

Lemma f0eps_lim_is_lim alpha beta epsilon (Halpha : (-1 < alpha)%Z) (Heps : 0 < epsilon) :
  filterlim (fun x : Rf0eps alpha beta epsilon (/ x))
            (at_right 0) (locally (f0eps_lim alpha beta epsilon)).

Lemma f0eps_lim_is_lim_sing alpha beta epsilon sing (Halpha : (-1 < alpha)%Z) (Heps : 0 < epsilon) :
  filterlim (fun x : Rf0eps alpha beta epsilon (/ (x - sing)))
            (at_right sing) (locally (f0eps_lim alpha beta epsilon)).

Lemma f0eps_lim_correct alpha beta epsilon (Halpha : (-1 < alpha)%Z) (Heps : 0 < epsilon) :
  is_RInt_gen ((fun xpowerRZ x alpha × (pow (ln x) beta))) (at_right 0) (at_point epsilon) (f0eps_lim alpha beta epsilon).

Lemma f0eps_lim_correct_sing alpha beta epsilon sing (Halpha : (-1 < alpha)%Z) (Heps : 0 < epsilon) :
  is_RInt_gen ((fun xpowerRZ (x - sing) alpha × (pow (ln (x - sing)) beta))) (at_right sing) (at_point (sing + epsilon)) (f0eps_lim alpha beta epsilon).

End ZeroToEpsilon.

Module BertrandInterval (I : IntervalOps).

Module J := IntervalExt I.

Section EffectiveBertrand.

Variable prec : I.precision.

Section Infinity.

Variable a : R.
Variable A : I.type.
Let iA := I.convert A.

Hypothesis Hcontainsa : contains iA (Xreal a).

Section BertrandLogNegInt.

Definition f_neg_int beta :=
  I.inv prec (I.mul prec (I.fromZ prec (Z.of_nat beta)) (I.power_int prec (I.ln prec A) (Z.of_nat beta))).

Lemma f_neg_int_correct beta : contains (I.convert (f_neg_int beta)) (Xreal (- f_neg a beta)).

End BertrandLogNegInt.

Fixpoint f_int_aux (alpha : Z) (beta : nat) (A_pow_Salpha : I.type) (ln_A : I.type) {struct beta} : I.type :=
  let alphap1 := I.fromZ prec (alpha + 1) in
  match beta with
  | 0 ⇒ I.div prec (I.neg A_pow_Salpha) alphap1
  | S m
    let beta := Z.of_nat beta in
    I.sub prec (I.div prec (I.neg (I.mul prec A_pow_Salpha (I.power_int prec ln_A beta))) alphap1)
      (I.mul prec (I.div prec (I.fromZ prec beta) alphap1) (f_int_aux alpha m A_pow_Salpha ln_A))
  end.

Definition f_int_fast (alpha : Z) (beta : nat) :=
  let A_pow_Salpha := I.power_int prec A (alpha+1) in
  let ln_A := I.ln prec A in
  f_int_aux alpha beta A_pow_Salpha ln_A.

Fixpoint f_int (alpha : Z) (beta : nat) {struct beta} : I.type :=
  let alphap1' := (alpha + 1)%Z in
  let alphap1 := I.fromZ prec alphap1' in
  match beta with
  | 0 ⇒ I.div prec (I.neg (I.power_int prec A alphap1')) alphap1
  | S m
    let beta := Z.of_nat beta in
    I.sub prec (I.div prec (I.neg (I.mul prec (I.power_int prec A alphap1') (I.power_int prec (I.ln prec A) beta))) alphap1)
      (I.mul prec (I.div prec (I.fromZ prec beta) alphap1) (f_int alpha m))
  end.

Lemma f_int_correct alpha beta (H : 0 < a) (Halpha: alpha (-1)%Z) :
  contains (I.convert (f_int alpha beta)) (Xreal (f_lim alpha beta a)).

Lemma f_int_fast_f_int alpha beta : f_int_fast alpha beta = f_int alpha beta.

Lemma f_int_fast_correct alpha beta (H : 0 < a) (Halpha: alpha (-1)%Z) :
  contains (I.convert (f_int_fast alpha beta)) (Xreal (f_lim alpha beta a)).

End Infinity.

Section Sing.

Variable epsilon : R.
Variable Epsilon : I.type.
Let iEps := I.convert Epsilon.

Hypothesis HEps : contains iEps (Xreal epsilon).
Hypothesis eps_gt0 : 0 < epsilon.

Definition f0eps_int (alpha : Z) (beta : nat) :=
  let yi := f_int_fast (I.inv prec Epsilon) (- 2 - alpha) beta in
  if Nat.even beta then yi else I.neg yi.

Lemma f0eps_correct (alpha : Z) (beta : nat) (Halpha : (alpha -1)%Z) :
  contains (I.convert (f0eps_int alpha beta)) (Xreal (f0eps_lim alpha beta epsilon)).

End Sing.

End EffectiveBertrand.

End BertrandInterval.