Require Import Reals.
Require Import Fourier.

Require Import Machin.
Check PI_2_3_7_ineq.
Open Scope R_scope.

(* This tactic expands the various functions and does the computation with
   small natural numbers.  Then it replaces all R constants with Z constants
 *)

Ltac expand1 :=
  simpl; unfold tg_alt, PI_2_3_7_tg, Ratan_seq; simpl plus; simpl mult;
  rewrite <- !(tech_pow_Rmult (-1)); simpl ((-1)^0); rewrite !INR_IZR_INZ;
  repeat rewrite ?Ropp_mult_distr_l_reverse, ?Rmult_1_l,
    ?Ropp_involutive; simpl Z_of_nat; replace 3 with (IZR 3)
    by (simpl; ring); replace 2 with (IZR 2) by (simpl; ring);
    replace 1 with (IZR 1) by reflexivity;
    repeat rewrite <- ?mult_IZR, <- ?plus_IZR;
    repeat match goal with
      |- context[IZR (Zplus ?a ?b)] => let v := eval compute in (Zplus a b) in
         change (Zplus a b) with v
    | |- context[IZR (Zmult ?a ?b)] => let v := eval compute in (Zmult a b) in
         change (Zmult a b) with v
    end.

(* This tactic combines together all inverses of natural numbers, coming
  from the coefficients or from the arguments raised to a power. It only
  works for inverses of integers. *)
Ltac expand2 :=
assert (u : forall a, a / IZR 1 = a)
  by (intros; simpl (IZR 1); unfold Rdiv; rewrite Rinv_1, Rmult_1_r;
         reflexivity);
rewrite !u, !pow_1, <- !Rinv_pow, !pow_IZR;
    try (apply Rgt_not_eq; apply (IZR_lt 0); reflexivity);
   repeat match goal with |- context[Zpow_def.Zpower ?a ?b] =>
     let v := eval compute in (Zpow_def.Zpower a b) in
     change (Zpow_def.Zpower a b) with v
   end;
assert (u1 : forall a b, /IZR (Zpos a) /IZR (Zpos b) = /IZR (Zpos (a * b)));
    [intros a b; unfold Rdiv; rewrite <- Rinv_mult_distr, <- mult_IZR;
     simpl Zmult; [reflexivity | 
    apply Rgt_not_eq, (IZR_lt 0); reflexivity |
   apply Rgt_not_eq, (IZR_lt 0); reflexivity ] | ];
assert (u2 : forall a b, /IZR (Zpos a) * /IZR (Zpos b) = /IZR (Zpos (a * b)));
    [intros a b; unfold Rdiv; rewrite <- Rinv_mult_distr, <- mult_IZR;
     simpl Zmult; [reflexivity | 
    apply Rgt_not_eq, (IZR_lt 0); reflexivity |
   apply Rgt_not_eq, (IZR_lt 0); reflexivity ] | ];
  repeat rewrite ?u1, ?u2; simpl Zplus; simpl Zmult; clear u u1 u2.

(* This tactic processes the sums. *)
Ltac expand3 :=
assert (u : forall a b c, IZR a * / IZR (Zpos b) + / IZR (Zpos c) =
              IZR (a * Zpos c + Zpos b) / IZR (Zpos (b * c)));[
  intros a b c; change (Zpos (b * c)) with (Zpos b * Zpos c)%Z;
  rewrite plus_IZR, !mult_IZR; field; split;
    apply Rgt_not_eq, (IZR_lt 0); reflexivity|rewrite !u; simpl Zplus];
  assert (u1 : forall a b, -(IZR a / b) = IZR (-a) / b);[
   intros a b; unfold Rdiv; 
   rewrite <- Ropp_mult_distr_l_reverse, opp_IZR; reflexivity|
   rewrite ?u1]; clear u u1.

Ltac expand4 :=
assert (u : forall a b c d, a / IZR (Zpos b) + c / IZR (Zpos d) =
               (a * IZR (Zpos d) + c * IZR (Zpos b))/(IZR (Zpos (b * d))));
 [intros a b c d; change (Zpos (b * d)) with (Zpos b * Zpos d)%Z;
 rewrite mult_IZR; field; split; apply Rgt_not_eq, (IZR_lt 0); reflexivity |];
 rewrite !u; repeat rewrite <- ?plus_IZR, <- ?mult_IZR;
 apply Rminus_lt; unfold Rminus; rewrite <- Ropp_div, <- opp_IZR, u,
  <- !mult_IZR, <- !plus_IZR;
assert (u1 : forall a b, (a < 0)%Z -> IZR a / IZR (Zpos b) < 0);[
  intros a b aneg; rewrite <- (Rmult_0_l (/(IZR (Zpos b)))); unfold Rdiv;
   apply Rmult_lt_compat_r;
   [apply Rinv_0_lt_compat, (IZR_lt 0);reflexivity |
    apply (fun x => IZR_lt x 0); assumption] | apply u1; reflexivity].

Lemma PI_close : 314/100 < PI < 315/100 .
destruct (PI_2_3_7_ineq 2) as [_ t1].
destruct (PI_2_3_7_ineq 2) as [t2 _].
split.
 assert (314/400 < PI/4).
  apply Rlt_le_trans with (2 := t2).
  expand1; expand2; expand3; expand4.
fourier.
assert (PI/4 < 315/400).
 apply Rle_lt_trans with (1 := t1).
 expand1; expand2; expand3; expand4.
fourier.
Qed.

Lemma PI_close' : 3141/1000 < PI < 3142/1000 .
destruct (PI_2_3_7_ineq 2) as [_ t1].
destruct (PI_2_3_7_ineq 2) as [t2 _].
split.
 assert (3141/4000 < PI/4).
  apply Rlt_le_trans with (2 := t2).
  expand1; expand2; expand3; expand4.
fourier.
assert (PI/4 < 3142/4000).
 apply Rle_lt_trans with (1 := t1).
 expand1; expand2; expand3; expand4.
fourier.
Qed.