Return to pi.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/pi.c, Revision 1.1

1.1     ! maekawa     1: /* mpfr_pi -- compute Pi
        !             2:
        !             3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
        !             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 Library General Public License as published by
        !             9: the Free Software Foundation; either version 2 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 Library General Public
        !            15: License for more details.
        !            16:
        !            17: You should have received a copy of the GNU Library 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 <math.h>
        !            24: #include "gmp.h"
        !            25: #include "gmp-impl.h"
        !            26: #include "longlong.h"
        !            27: #include "mpfr.h"
        !            28:
        !            29: /*
        !            30: Set x to the value of Pi to precision PREC(x) rounded to direction rnd_mode.
        !            31: Use the formula giving the binary representation of Pi found by Simon Plouffe
        !            32: and the Borwein's brothers:
        !            33:
        !            34:                    infinity    4         2         1         1
        !            35:                     -----   ------- - ------- - ------- - -------
        !            36:                      \      8 n + 1   8 n + 4   8 n + 5   8 n + 6
        !            37:               Pi =    )     -------------------------------------
        !            38:                      /                         n
        !            39:                     -----                    16
        !            40:                     n = 0
        !            41:
        !            42: i.e. Pi*16^N = S(N) + R(N) where
        !            43: 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)
        !            44: 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)
        !            45:
        !            46: Let f(n) = 4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6), we can show easily that
        !            47: f(n) < 15/(64*n^2), so R(N) < sum(15/(64*n^2)/16^(n-N), n=N..infinity)
        !            48:                             < 15/64/N^2*sum(1/16^(n-N), n=N..infinity)
        !            49:                            = 1/4/N^2
        !            50:
        !            51: Now let S'(N) = sum(floor(16^(N-n)*(120*n^2+151*n+47),
        !            52:   (512*n^4+1024*n^3+712*n^2+194*n+15)), n=0..N-1)
        !            53:
        !            54: S(N)-S'(N) <= sum(1, n=0..N-1) = N
        !            55:
        !            56: so Pi*16^N-S'(N) <= N+1 (as 1/4/N^2 < 1)
        !            57: */
        !            58:
        !            59: mpfr_t __mpfr_pi; /* stored value of Pi */
        !            60: int __mpfr_pi_prec=0; /* precision of stored value */
        !            61: char __mpfr_pi_rnd; /* rounding mode of stored value */
        !            62:
        !            63: void
        !            64: #if __STDC__
        !            65: mpfr_pi(mpfr_ptr x, unsigned char rnd_mode)
        !            66: #else
        !            67: mpfr_pi(x, rnd_mode)
        !            68:      mpfr_ptr x;
        !            69:      unsigned char rnd_mode;
        !            70: #endif
        !            71: {
        !            72:   int N, oldN, n, prec; mpz_t pi, num, den, d3, d2, tmp; mpfr_t y;
        !            73:
        !            74:   prec=PREC(x);
        !            75:
        !            76:   /* has stored value enough precision ? */
        !            77:   if ((prec==__mpfr_pi_prec && rnd_mode==__mpfr_pi_rnd) ||
        !            78:       (prec<=__mpfr_pi_prec &&
        !            79:       mpfr_can_round(__mpfr_pi, __mpfr_pi_prec, __mpfr_pi_rnd, rnd_mode, prec)))
        !            80:     {
        !            81:       mpfr_set(x, __mpfr_pi, rnd_mode); return;
        !            82:     }
        !            83:
        !            84:   /* need to recompute */
        !            85:   N=1;
        !            86:   do {
        !            87:     oldN = N;
        !            88:     N = (prec+3)/4 + (int)ceil(log((double)N+1.0)/log(2.0));
        !            89:   } while (N != oldN);
        !            90:   mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2);
        !            91:   mpz_init(tmp);
        !            92:   mpz_set_ui(pi, 0);
        !            93:   mpz_set_ui(num, 16); /* num(-1) */
        !            94:   mpz_set_ui(den, 21); /* den(-1) */
        !            95:   mpz_set_si(d3, -2454);
        !            96:   mpz_set_ui(d2, 14736);
        !            97:   /* invariants: num=120*n^2+151*n+47, den=512*n^4+1024*n^3+712*n^2+194*n+15
        !            98:                  d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448
        !            99:    */
        !           100:   for (n=0; n<N; n++) {
        !           101:     /* num(n)-num(n-1) = 240*n+31 */
        !           102:     mpz_add_ui(num, num, 240*n+31); /* no overflow up to PREC=71M */
        !           103:     /* d2(n) - d2(n-1) = 12288*(n-1) */
        !           104:     if (n>0) mpz_add_ui(d2, d2, 12288*(n-1));
        !           105:     else mpz_sub_ui(d2, d2, 12288);
        !           106:     /* d3(n) - d3(n-1) = d2 */
        !           107:     mpz_add(d3, d3, d2);
        !           108:     /* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */
        !           109:     mpz_add(den, den, d3);
        !           110:     mpz_mul_2exp(tmp, num, 4*(N-n));
        !           111:     mpz_fdiv_q(tmp, tmp, den);
        !           112:     mpz_add(pi, pi, tmp);
        !           113:   }
        !           114:   mpfr_set_z(x, pi, rnd_mode);
        !           115:   mpfr_init2(y, mpfr_get_prec(x));
        !           116:   mpz_add_ui(pi, pi, N+1);
        !           117:   mpfr_set_z(y, pi, rnd_mode);
        !           118:   if (mpfr_cmp(x, y) != 0) {
        !           119:     fprintf(stderr, "does not converge\n"); exit(1);
        !           120:   }
        !           121:   EXP(x) -= 4*N;
        !           122:   mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2);
        !           123:   mpz_clear(tmp); mpfr_clear(y);
        !           124:
        !           125:   /* store computed value */
        !           126:   if (__mpfr_pi_prec==0) mpfr_init2(__mpfr_pi, prec);
        !           127:   else mpfr_set_prec(__mpfr_pi, prec);
        !           128:   mpfr_set(__mpfr_pi, x, rnd_mode);
        !           129:   __mpfr_pi_prec=prec;
        !           130:   __mpfr_pi_rnd=rnd_mode;
        !           131: }