Require Import Reals Ranalysis5 Psatz ClassicalEpsilon.
Require Import Rcomplements RInt.
Require Import Coquelicot.

Lemma Rlt_1_2 : 1 < 2.  Proof. psatzl R. Qed.

Ltac lt0 :=
  match goal with
  | |- _ => assumption
  | |- 0 < exp _ => apply exp_pos
  | |- 0 < pos _ => apply cond_pos
  | |- _ > 0 => unfold Rgt; lt0
  | |- 0 < 1 =>  apply Rlt_0_1
  | |- 0 < 2 => apply Rlt_0_2
  | |- 0 < PI => apply PI_RGT_0
  | |- _ <> 0 => apply Rgt_not_eq; lt0
  | |- 0 < Rabs _ + _ => apply (Rplus_le_lt_0_compat _ _ (Rabs_pos _)); lt0
  | |- 0 < _ * _ => apply Rmult_lt_0_compat; lt0
  | |- 0 < _ + _ => apply Rplus_lt_0_compat; lt0
  | |- 0 < Rmin _ _ => apply Rmin_glb_lt; lt0
  | |- 0 < / _ => apply Rinv_0_lt_compat; lt0
  | |- 0 < sqrt _ => apply sqrt_lt_R0; lt0
  | |- 0 < _ / _ => unfold Rdiv; apply Rmult_lt_0_compat; lt0
  | |- 0 < _ ^ _ => apply pow_lt; lt0
  | |- 0 < _ ^ 2 + _ => apply Rplus_le_lt_0_compat;[apply pow2_ge_0 | lt0]
  | |- _ <= _ => apply Rlt_le; lt0
  | |- _ => psatzl R
  end.

(* Starting the work on algebraic geometric means. *)

Fixpoint ag (a b : R) (n : nat) :=
  match n with
    0%nat => (a, b) 
  | S p => let (a_p, b_p) := ag a b p in
             ((a_p + b_p)/2, sqrt(a_p * b_p))
  end.

Definition u_ (n : nat) (x : R) := fst (ag 1 x n).

Lemma Rpower_mult_distr :  (* standard *)
  forall x y z, 0 < x -> 0 < y ->
   Rpower x z * Rpower y z = Rpower (x * y) z.
intros x y z x0 y0; unfold Rpower.
rewrite <- exp_plus, ln_mult, Rmult_plus_distr_l; auto.
Qed.

Definition v_ (n : nat) (x : R) := snd (ag 1 x n).

Definition y_ n x := u_ n x / v_ n x.

Axiom y_0 : y_ 0 (/sqrt 2) = sqrt 2.

Definition z_ n x := Derive (v_ n) x / Derive (u_ n) x.

(* standard. *)
Lemma pow_lt_impl :
  forall n x y, 0 < x -> 0 < y -> x ^ n < y ^ n -> x < y.
Proof.
intros n x y x0 y0 a.
destruct (Rle_lt_dec y x) as [yx | xy];[ | auto].
assert (t := pow_incr y x n (conj (Rlt_le _ _ y0) yx)).
destruct (Rle_not_lt _ _ t); assumption.
Qed.

(* standard *)
Lemma pow_increasing n x y : (0 < n)%nat -> 0 < x -> x < y -> x ^ n < y ^ n.
induction 1 as [ | m mgt0 IH]; [simpl; rewrite !Rmult_1_r; auto | ].
intros x0 xy; simpl; apply Rmult_gt_0_lt_compat.
   apply pow_lt; assumption.
  apply (Rlt_trans _ _ _ x0 xy).
 assumption.
apply IH; assumption.
Qed.

(* standard *)
Lemma pow_eq_impl :
  forall n x y, (0 < n)%nat -> 0 < x -> 0 < y -> x ^ n = y ^ n -> x = y.
Proof.
intros n x y n0 x0 y0 a.
destruct (Rle_lt_dec y x) as [yx | xy].
 destruct yx as [yx | yx]; auto.
 assert (t := pow_increasing n y x n0 y0 yx); rewrite a in t.
 now case (Rlt_irrefl  _ t).
assert (t := pow_increasing n x y n0 x0 xy); rewrite a in t.
now case (Rlt_irrefl _ t).
Qed.

Lemma pow2_sqrt : forall x, 0 <= x -> sqrt x ^ 2 = x.
 intros x x0; simpl; rewrite Rmult_1_r, sqrt_sqrt; auto.
Qed.


(* standard.  *)
Lemma sqrt_Rinv x : 0 < x -> sqrt (/ x) = / sqrt x.
intros x0; apply (pow_eq_impl 2).
   apply le_S, le_n.
  apply sqrt_lt_R0, Rinv_0_lt_compat; assumption.
 apply Rinv_0_lt_compat, sqrt_lt_R0; assumption.
rewrite pow2_sqrt, <- Rinv_pow, pow2_sqrt.
   reflexivity.
  apply Rlt_le; assumption.
 apply Rgt_not_eq, sqrt_lt_R0; assumption.
apply Rlt_le, Rinv_0_lt_compat; assumption.
Qed.

(* to be added in coquelicot *)
Axiom Derive_sqrt : forall x, 0 < x -> Derive sqrt x = /(2 * sqrt x).

Axiom z_1 : z_ 1 (/sqrt 2) = sqrt (sqrt 2).

Axiom q4_1_1 : forall x n, 0 < x < 1 -> 1 < y_ n x.

Axiom q4_1_2y : forall x n, 0 < x < 1 -> 
   y_ (n + 1) x = (1 + y_ n x) /(2 * (sqrt (y_ n x))).

Axiom q4_1_2z : forall x n, (1 <= n)%nat -> 0 < x < 1 ->
   z_ (n + 1) x = (1 + y_ n x * z_ n x)/((1 + z_ n x) * sqrt (y_ n x)).

Axiom q4_1_3 : forall x n, 0 < x < 1 -> (1 <= n)%nat -> 1 <= z_ n x.

Axiom q4_1_4 : forall x n, 0 < x < 1 -> (1 <= n)%nat ->
  y_ (n + 1) x <= z_ (n + 1) x <= sqrt (y_ n x).

Fixpoint agmpi n :=
  match n with
    0%nat => (2 + sqrt 2)
  | S p => agmpi p * (1 + y_ n (/sqrt 2)) / (1 + z_ n (/sqrt 2))
  end.

Axiom q4_3_1 : forall n a, 0 < a < 1 ->
    0  <= y_ (n+1) a - 1 <= (y_ n a - 1)^2/8.

Axiom q4_3_2 : forall n, (1 <= n)%nat ->
  0 <= agmpi (n + 1) - PI <= 4 * agmpi 0 * Rpower 500 (- 2 ^ n).

Axiom iter_million : 0 <= agmpi 20 - PI <= Rpower 10 (-(10 ^ 6 + 10 ^ 4)).