1
/* mpfr_const_pi -- compute Pi
3
Copyright 1999, 2000, 2001 Free Software Foundation.
5
This file is part of the MPFR Library.
7
The MPFR Library is free software; you can redistribute it and/or modify
8
it under the terms of the GNU Lesser General Public License as published by
9
the Free Software Foundation; either version 2.1 of the License, or (at your
10
option) any later version.
12
The MPFR Library is distributed in the hope that it will be useful, but
13
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15
License for more details.
17
You should have received a copy of the GNU Lesser General Public License
18
along with the MPFR Library; see the file COPYING.LIB. If not, write to
19
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20
MA 02111-1307, USA. */
28
#include "mpfr-impl.h"
30
static int mpfr_aux_pi _PROTO ((mpfr_ptr, mpz_srcptr, int, int));
31
static int mpfr_pi_machin3 _PROTO ((mpfr_ptr, mp_rnd_t));
40
#define GENERIC mpfr_aux_pi
47
mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode)
49
int prec, logn, prec_x;
50
int prec_i_want=MPFR_PREC(mylog);
52
mpfr_t tmp1, tmp2, result,tmp3,tmp4,tmp5,tmp6;
55
MPFR_CLEAR_FLAGS(mylog);
56
logn = _mpfr_ceil_log2 ((double) MPFR_PREC(mylog));
57
prec_x = prec_i_want + logn + 5;
60
prec = _mpfr_ceil_log2 ((double) prec_x);
62
mpfr_init2(tmp1, prec_x);
63
mpfr_init2(tmp2, prec_x);
64
mpfr_init2(tmp3, prec_x);
65
mpfr_init2(tmp4, prec_x);
66
mpfr_init2(tmp5, prec_x);
67
mpfr_init2(tmp6, prec_x);
68
mpfr_init2(result, prec_x);
71
mpfr_aux_pi(tmp1, cst, 268*268, prec - 4);
72
mpfr_div_ui(tmp1, tmp1, 268, GMP_RNDD);
73
mpfr_mul_ui(tmp1, tmp1, 61, GMP_RNDD);
75
mpfr_aux_pi(tmp2, cst, 343*343, prec - 4);
76
mpfr_div_ui(tmp2, tmp2, 343, GMP_RNDD);
77
mpfr_mul_ui(tmp2, tmp2, 122, GMP_RNDD);
79
mpfr_aux_pi(tmp3, cst, 557*557, prec - 4);
80
mpfr_div_ui(tmp3, tmp3, 557, GMP_RNDD);
81
mpfr_mul_ui(tmp3, tmp3, 115, GMP_RNDD);
83
mpfr_aux_pi(tmp4, cst, 1068*1068, prec - 4);
84
mpfr_div_ui(tmp4, tmp4, 1068, GMP_RNDD);
85
mpfr_mul_ui(tmp4, tmp4, 32, GMP_RNDD);
87
mpfr_aux_pi(tmp5, cst, 3458*3458, prec - 4);
88
mpfr_div_ui(tmp5, tmp5, 3458, GMP_RNDD);
89
mpfr_mul_ui(tmp5, tmp5, 83, GMP_RNDD);
91
mpfr_aux_pi(tmp6, cst, 27493*27493, prec - 4);
92
mpfr_div_ui(tmp6, tmp6, 27493, GMP_RNDD);
93
mpfr_mul_ui(tmp6, tmp6, 44, GMP_RNDD);
95
mpfr_add(result, tmp1, tmp2, GMP_RNDD);
96
mpfr_add(result, result, tmp3, GMP_RNDD);
97
mpfr_sub(result, result, tmp4, GMP_RNDD);
98
mpfr_add(result, result, tmp5, GMP_RNDD);
99
mpfr_add(result, result, tmp6, GMP_RNDD);
100
mpfr_mul_2ui(result, result, 2, GMP_RNDD);
107
if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)){
108
mpfr_set(mylog, result, rnd_mode);
122
Set x to the value of Pi to precision MPFR_PREC(x) rounded to direction rnd_mode.
123
Use the formula giving the binary representation of Pi found by Simon Plouffe
124
and the Borwein's brothers:
127
----- ------- - ------- - ------- - -------
128
\ 8 n + 1 8 n + 4 8 n + 5 8 n + 6
129
Pi = ) -------------------------------------
134
i.e. Pi*16^N = S(N) + R(N) where
135
S(N) = sum(16^(N-n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)), n=0..N-1)
136
R(N) = sum((4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))/16^(n-N), n=N..infinity)
138
Let f(n) = 4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6), we can show easily that
139
f(n) < 15/(64*n^2), so R(N) < sum(15/(64*n^2)/16^(n-N), n=N..infinity)
140
< 15/64/N^2*sum(1/16^(n-N), n=N..infinity)
143
Now let S'(N) = sum(floor(16^(N-n)*(120*n^2+151*n+47),
144
(512*n^4+1024*n^3+712*n^2+194*n+15)), n=0..N-1)
146
S(N)-S'(N) <= sum(1, n=0..N-1) = N
148
so Pi*16^N-S'(N) <= N+1 (as 1/4/N^2 < 1)
151
mpfr_t __mpfr_const_pi; /* stored value of Pi */
152
int __mpfr_const_pi_prec=0; /* precision of stored value */
153
mp_rnd_t __mpfr_const_pi_rnd; /* rounding mode of stored value */
156
mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode)
158
int N, oldN, n, prec; mpz_t pi, num, den, d3, d2, tmp; mpfr_t y;
162
/* has stored value enough precision ? */
163
if ((prec==__mpfr_const_pi_prec && rnd_mode==__mpfr_const_pi_rnd) ||
164
(prec<=__mpfr_const_pi_prec &&
165
mpfr_can_round(__mpfr_const_pi, __mpfr_const_pi_prec,
166
__mpfr_const_pi_rnd, rnd_mode, prec)))
168
mpfr_set(x, __mpfr_const_pi, rnd_mode); return;
172
/* need to recompute */
176
N = (prec+3)/4 + _mpfr_ceil_log2((double) N + 1.0);
178
mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2);
181
mpz_set_ui(num, 16); /* num(-1) */
182
mpz_set_ui(den, 21); /* den(-1) */
183
mpz_set_si(d3, -2454);
184
mpz_set_ui(d2, 14736);
185
/* invariants: num=120*n^2+151*n+47, den=512*n^4+1024*n^3+712*n^2+194*n+15
186
d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448
188
for (n=0; n<N; n++) {
189
/* num(n)-num(n-1) = 240*n+31 */
190
mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */
191
/* d2(n) - d2(n-1) = 12288*(n-1) */
192
if (n>0) mpz_add_ui(d2, d2, 12288*(n-1));
193
else mpz_sub_ui(d2, d2, 12288);
194
/* d3(n) - d3(n-1) = d2 */
196
/* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */
197
mpz_add(den, den, d3);
198
mpz_mul_2exp(tmp, num, 4*(N-n));
199
mpz_fdiv_q(tmp, tmp, den);
200
mpz_add(pi, pi, tmp);
202
mpfr_set_z(x, pi, rnd_mode);
203
mpfr_init2(y, mpfr_get_prec(x));
204
mpz_add_ui(pi, pi, N+1);
205
mpfr_set_z(y, pi, rnd_mode);
206
if (mpfr_cmp(x, y) != 0) {
207
fprintf(stderr, "does not converge\n"); exit(1);
210
mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2);
211
mpz_clear(tmp); mpfr_clear(y);
213
mpfr_pi_machin3(x, rnd_mode);
214
/* store computed value */
215
if (__mpfr_const_pi_prec==0) mpfr_init2(__mpfr_const_pi, prec);
216
else mpfr_set_prec(__mpfr_const_pi, prec);
217
mpfr_set(__mpfr_const_pi, x, rnd_mode);
218
__mpfr_const_pi_prec=prec;
219
__mpfr_const_pi_rnd=rnd_mode;
added
removed
src/gmp/mpfr/const_pi.c
