aboutsummaryrefslogtreecommitdiffhomepage
path: root/theories/Reals/Machin.v
blob: 8f7e07ac4d860a3d906dbd8a4980fb81bba1fb81 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
(************************************************************************)
(*         *   The Coq Proof Assistant / The Coq Development Team       *)
(*  v      *   INRIA, CNRS and contributors - Copyright 1999-2018       *)
(* <O___,, *       (see CREDITS file for the list of authors)           *)
(*   \VV/  **************************************************************)
(*    //   *    This file is distributed under the terms of the         *)
(*         *     GNU Lesser General Public License Version 2.1          *)
(*         *     (see LICENSE file for the text of the license)         *)
(************************************************************************)

Require Import Lra.
Require Import Rbase.
Require Import Rtrigo1.
Require Import Ranalysis_reg.
Require Import Rfunctions.
Require Import AltSeries.
Require Import Rseries.
Require Import SeqProp.
Require Import PartSum.
Require Import Ratan.
Require Import Omega.

Local Open Scope R_scope.

(* Proving a few formulas in the style of John Machin to compute Pi *)

Definition atan_sub u v := (u - v)/(1 + u * v).

Lemma atan_sub_correct :
 forall u v, 1 + u * v <> 0 -> -PI/2 < atan u - atan v < PI/2 ->
   -PI/2 < atan (atan_sub u v) < PI/2 ->
   atan u = atan v + atan (atan_sub u v).
Proof.
intros u v pn0 uvint aint.
assert (cos (atan u) <> 0).
 destruct (atan_bound u); apply Rgt_not_eq, cos_gt_0; auto.
 rewrite <- Ropp_div; assumption.
assert (cos (atan v) <> 0).
 destruct (atan_bound v); apply Rgt_not_eq, cos_gt_0; auto.
 rewrite <- Ropp_div; assumption.
assert (t : forall a b c, a - b = c -> a = b + c) by (intros; subst; field).
apply t, tan_is_inj; clear t; try assumption.
rewrite tan_minus; auto.
  rewrite !atan_right_inv; reflexivity.
apply Rgt_not_eq, cos_gt_0; rewrite <- ?Ropp_div; tauto.
rewrite !atan_right_inv; assumption.
Qed.

Lemma tech : forall x y , -1 <= x <= 1 -> -1 < y < 1 -> 
  -PI/2 < atan x - atan y < PI/2.
Proof.
assert (ut := PI_RGT_0).
intros x y [xm1 x1] [ym1 y1].
assert (-(PI/4) <= atan x).
 destruct xm1 as [xm1 | xm1].
  rewrite <- atan_1, <- atan_opp; apply Rlt_le, atan_increasing.
  assumption.
 solve[rewrite <- xm1; change (-1) with (-(1)); rewrite atan_opp, atan_1; apply Rle_refl].
assert (-(PI/4) < atan y).
 rewrite <- atan_1, <- atan_opp; apply atan_increasing.
 assumption.
assert (atan x <= PI/4).
 destruct x1 as [x1 | x1].
  rewrite <- atan_1; apply Rlt_le, atan_increasing.
  assumption.
 solve[rewrite x1, atan_1; apply Rle_refl].
assert (atan y < PI/4).
  rewrite <- atan_1; apply atan_increasing.
  assumption.
rewrite Ropp_div; split; lra.
Qed.

(* A simple formula, reasonably efficient. *)
Lemma Machin_2_3 : PI/4 = atan(/2) + atan(/3).
Proof.
assert (utility : 0 < PI/2) by (apply PI2_RGT_0).
rewrite <- atan_1.
rewrite (atan_sub_correct 1 (/2)).
   apply f_equal, f_equal; unfold atan_sub; field.
  apply Rgt_not_eq; lra.
 apply tech; try split; try lra.
apply atan_bound.
Qed.

Lemma Machin_4_5_239 : PI/4 = 4 * atan (/5) - atan(/239).
Proof.
rewrite <- atan_1.
rewrite (atan_sub_correct 1 (/5));
 [ | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
replace (4 * atan (/5) - atan (/239)) with
  (atan (/5) + (atan (/5) + (atan (/5) + (atan (/5) + -
     atan (/239))))) by ring.
apply f_equal.
replace (atan_sub 1 (/5)) with (2/3) by
  (unfold atan_sub; field).
rewrite (atan_sub_correct (2/3) (/5));
 [apply f_equal | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
replace (atan_sub (2/3) (/5)) with (7/17)  by
  (unfold atan_sub; field).
rewrite (atan_sub_correct (7/17) (/5));
 [apply f_equal | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
replace (atan_sub (7/17) (/5)) with (9/46) by
  (unfold atan_sub; field).
rewrite (atan_sub_correct (9/46) (/5));
 [apply f_equal | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
rewrite <- atan_opp; apply f_equal.
unfold atan_sub; field.
Qed.

Lemma Machin_2_3_7 : PI/4 = 2 * atan(/3) + (atan (/7)).
Proof.
rewrite <- atan_1.
rewrite (atan_sub_correct 1 (/3));
 [ | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
replace (2 * atan (/3) + atan (/7)) with
  (atan (/3) + (atan (/3) + atan (/7))) by ring.
apply f_equal.
replace (atan_sub 1 (/3)) with (/2) by
  (unfold atan_sub; field).
rewrite (atan_sub_correct (/2) (/3));
 [apply f_equal | apply Rgt_not_eq; lra | apply tech; try split; lra |
     apply atan_bound ].
apply f_equal; unfold atan_sub; field.
Qed.

(* More efficient way to compute approximations of PI. *)

Definition PI_2_3_7_tg n := 
  2 * Ratan_seq (/3) n + Ratan_seq (/7) n.

Lemma PI_2_3_7_ineq :
  forall N : nat,
    sum_f_R0 (tg_alt PI_2_3_7_tg) (S (2 * N)) <= PI / 4 <=
    sum_f_R0 (tg_alt PI_2_3_7_tg) (2 * N).
Proof.
assert (dec3 : 0 <= /3 <= 1) by (split; lra).
assert (dec7 : 0 <= /7 <= 1) by (split; lra).
assert (decr : Un_decreasing PI_2_3_7_tg).
 apply Ratan_seq_decreasing in dec3.
 apply Ratan_seq_decreasing in dec7.
 intros n; apply Rplus_le_compat.
  apply Rmult_le_compat_l; [ lra | exact (dec3 n)].
 exact (dec7 n).
assert (cv : Un_cv PI_2_3_7_tg 0).
 apply Ratan_seq_converging in dec3.
 apply Ratan_seq_converging in dec7.
 intros eps ep.
 assert (ep' : 0 < eps /3) by lra.
 destruct (dec3 _ ep') as [N1 Pn1]; destruct (dec7 _ ep') as [N2 Pn2].
 exists (N1 + N2)%nat; intros n Nn.
 unfold PI_2_3_7_tg.
 rewrite <- (Rplus_0_l 0).
 apply Rle_lt_trans with
    (1 := R_dist_plus (2 * Ratan_seq (/3) n) 0 (Ratan_seq (/7) n) 0).
 replace eps with (2 * eps/3 + eps/3) by field.
 apply Rplus_lt_compat.
  unfold R_dist, Rminus, Rdiv.
  rewrite <- (Rmult_0_r 2), <- Ropp_mult_distr_r_reverse.
  rewrite <- Rmult_plus_distr_l, Rabs_mult, (Rabs_pos_eq 2);[|lra].
  rewrite Rmult_assoc; apply Rmult_lt_compat_l;[lra | ].
   apply (Pn1 n); omega.
  apply (Pn2 n); omega.
rewrite Machin_2_3_7.
rewrite !atan_eq_ps_atan; try (split; lra).
unfold ps_atan; destruct (in_int (/3)); destruct (in_int (/7));
  try match goal with id : ~ _ |- _ => case id; split; lra end.
destruct (ps_atan_exists_1 (/3)) as [v3 Pv3].
destruct (ps_atan_exists_1 (/7)) as [v7 Pv7].
assert (main : Un_cv (sum_f_R0 (tg_alt PI_2_3_7_tg)) (2 * v3 + v7)).
 assert (main :Un_cv (fun n => 2 * sum_f_R0 (tg_alt (Ratan_seq (/3))) n +
                        sum_f_R0 (tg_alt (Ratan_seq (/7))) n) (2 * v3 + v7)).
  apply CV_plus;[ | assumption].
  apply CV_mult;[ | assumption].
  exists 0%nat; intros; rewrite R_dist_eq; assumption.
 apply Un_cv_ext with (2 := main).
 intros n; rewrite scal_sum, <- plus_sum; apply sum_eq; intros.
 rewrite Rmult_comm; unfold PI_2_3_7_tg, tg_alt; field.
intros N; apply (alternated_series_ineq _ _ _ decr cv main).
Qed.