Annotation of OpenXM_contrib/gmp/mpfr/const_pi.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_const_pi -- compute Pi
! 2:
! 3: Copyright 1999, 2000, 2001 Free Software Foundation.
! 4:
! 5: This file is part of the MPFR Library.
! 6:
! 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.
! 11:
! 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.
! 16:
! 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. */
! 21:
! 22: #include <stdio.h>
! 23: #include <stdlib.h>
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "longlong.h"
! 27: #include "mpfr.h"
! 28: #include "mpfr-impl.h"
! 29:
! 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));
! 32:
! 33: #define A
! 34: #define A1 1
! 35: #define A2 2
! 36: #undef B
! 37: #define C
! 38: #define C1 3
! 39: #define C2 2
! 40: #define GENERIC mpfr_aux_pi
! 41: #define R_IS_RATIONAL
! 42: #define NO_FACTORIAL
! 43: #include "generic.c"
! 44:
! 45:
! 46: static int
! 47: mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode)
! 48: {
! 49: int prec, logn, prec_x;
! 50: int prec_i_want=MPFR_PREC(mylog);
! 51: int good = 0;
! 52: mpfr_t tmp1, tmp2, result,tmp3,tmp4,tmp5,tmp6;
! 53: mpz_t cst;
! 54:
! 55: MPFR_CLEAR_FLAGS(mylog);
! 56: logn = _mpfr_ceil_log2 ((double) MPFR_PREC(mylog));
! 57: prec_x = prec_i_want + logn + 5;
! 58: mpz_init(cst);
! 59: while (!good){
! 60: prec = _mpfr_ceil_log2 ((double) prec_x);
! 61:
! 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);
! 69: mpz_set_si(cst, -1);
! 70:
! 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);
! 74:
! 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);
! 78:
! 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);
! 82:
! 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);
! 86:
! 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);
! 90:
! 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);
! 94:
! 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);
! 101: mpfr_clear(tmp1);
! 102: mpfr_clear(tmp2);
! 103: mpfr_clear(tmp3);
! 104: mpfr_clear(tmp4);
! 105: mpfr_clear(tmp5);
! 106: mpfr_clear(tmp6);
! 107: if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)){
! 108: mpfr_set(mylog, result, rnd_mode);
! 109: mpfr_clear(result);
! 110: good = 1;
! 111: } else
! 112: {
! 113: mpfr_clear(result);
! 114: prec_x += logn;
! 115: }
! 116: }
! 117: mpz_clear(cst);
! 118: return 0;
! 119: }
! 120:
! 121: /*
! 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:
! 125:
! 126: infinity 4 2 1 1
! 127: ----- ------- - ------- - ------- - -------
! 128: \ 8 n + 1 8 n + 4 8 n + 5 8 n + 6
! 129: Pi = ) -------------------------------------
! 130: / n
! 131: ----- 16
! 132: n = 0
! 133:
! 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)
! 137:
! 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)
! 141: = 1/4/N^2
! 142:
! 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)
! 145:
! 146: S(N)-S'(N) <= sum(1, n=0..N-1) = N
! 147:
! 148: so Pi*16^N-S'(N) <= N+1 (as 1/4/N^2 < 1)
! 149: */
! 150:
! 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 */
! 154:
! 155: void
! 156: mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode)
! 157: {
! 158: int N, oldN, n, prec; mpz_t pi, num, den, d3, d2, tmp; mpfr_t y;
! 159:
! 160: prec=MPFR_PREC(x);
! 161:
! 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)))
! 167: {
! 168: mpfr_set(x, __mpfr_const_pi, rnd_mode); return;
! 169: }
! 170:
! 171: if (prec < 20000){
! 172: /* need to recompute */
! 173: N=1;
! 174: do {
! 175: oldN = N;
! 176: N = (prec+3)/4 + _mpfr_ceil_log2((double) N + 1.0);
! 177: } while (N != oldN);
! 178: mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2);
! 179: mpz_init(tmp);
! 180: mpz_set_ui(pi, 0);
! 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
! 187: */
! 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 */
! 195: mpz_add(d3, d3, 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);
! 201: }
! 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);
! 208: }
! 209: MPFR_EXP(x) -= 4*N;
! 210: mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2);
! 211: mpz_clear(tmp); mpfr_clear(y);
! 212: } else
! 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;
! 220: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to const_pi.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr