contrib/mpfr/const_pi.c

   1 /* mpfr_const_pi -- compute Pi
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #include "mpfr-impl.h"
  24
  25 /* Declare the cache */
  26 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_pi, mpfr_const_pi_internal);
  27
  28 /* Set User Interface */
  29 #undef mpfr_const_pi
  30 int
  31 mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) {
  32   return mpfr_cache (x, __gmpfr_cache_const_pi, rnd_mode);
  33 }
  34
  35 /* Don't need to save/restore exponent range: the cache does it */
  36 int
  37 mpfr_const_pi_internal (mpfr_ptr x, mp_rnd_t rnd_mode)
  38 {
  39   mpfr_t a, A, B, D, S;
  40   mp_prec_t px, p, cancel, k, kmax;
  41   MPFR_ZIV_DECL (loop);
  42   int inex;
  43
  44   MPFR_LOG_FUNC (("rnd_mode=%d", rnd_mode), ("x[%#R]=%R inex=%d", x, x, inex));
  45
  46   px = MPFR_PREC (x);
  47
  48   /* we need 9*2^kmax - 4 >= px+2*kmax+8 */
  49   for (kmax = 2; ((px + 2 * kmax + 12) / 9) >> kmax; kmax ++);
  50
  51   p = px + 3 * kmax + 14; /* guarantees no recomputation for px <= 10000 */
  52
  53   mpfr_init2 (a, p);
  54   mpfr_init2 (A, p);
  55   mpfr_init2 (B, p);
  56   mpfr_init2 (D, p);
  57   mpfr_init2 (S, p);
  58
  59   MPFR_ZIV_INIT (loop, p);
  60   for (;;) {
  61     mpfr_set_ui (a, 1, GMP_RNDN);          /* a = 1 */
  62     mpfr_set_ui (A, 1, GMP_RNDN);          /* A = a^2 = 1 */
  63     mpfr_set_ui_2exp (B, 1, -1, GMP_RNDN); /* B = b^2 = 1/2 */
  64     mpfr_set_ui_2exp (D, 1, -2, GMP_RNDN); /* D = 1/4 */
  65
  66 #define b B
  67 #define ap a
  68 #define Ap A
  69 #define Bp B
  70     for (k = 0, cancel = 0; ; k++)
  71       {
  72         /* invariant: 1/2 <= B <= A <= a < 1 */
  73         mpfr_add (S, A, B, GMP_RNDN); /* 1 <= S <= 2 */
  74         mpfr_div_2ui (S, S, 2, GMP_RNDN); /* exact, 1/4 <= S <= 1/2 */
  75         mpfr_sqrt (b, B, GMP_RNDN); /* 1/2 <= b <= 1 */
  76         mpfr_add (ap, a, b, GMP_RNDN); /* 1 <= ap <= 2 */
  77         mpfr_div_2ui (ap, ap, 1, GMP_RNDN); /* exact, 1/2 <= ap <= 1 */
  78         mpfr_mul (Ap, ap, ap, GMP_RNDN); /* 1/4 <= Ap <= 1 */
  79         mpfr_sub (Bp, Ap, S, GMP_RNDN); /* -1/4 <= Bp <= 3/4 */
  80         mpfr_mul_2ui (Bp, Bp, 1, GMP_RNDN); /* -1/2 <= Bp <= 3/2 */
  81         mpfr_sub (S, Ap, Bp, GMP_RNDN);
  82         MPFR_ASSERTN (mpfr_cmp_ui (S, 1) < 0);
  83         cancel = mpfr_cmp_ui (S, 0) ? (mpfr_uexp_t) -mpfr_get_exp(S) : p;
  84         /* MPFR_ASSERTN (cancel >= px || cancel >= 9 * (1 << k) - 4); */
  85         mpfr_mul_2ui (S, S, k, GMP_RNDN);
  86         mpfr_sub (D, D, S, GMP_RNDN);
  87         /* stop when |A_k - B_k| <= 2^(k-p) i.e. cancel >= p-k */
  88         if (cancel + k >= p)
  89           break;
  90       }
  91 #undef b
  92 #undef ap
  93 #undef Ap
  94 #undef Bp
  95
  96       mpfr_div (A, B, D, GMP_RNDN);
  97
  98       /* MPFR_ASSERTN(p >= 2 * k + 8); */
  99       if (MPFR_LIKELY (MPFR_CAN_ROUND (A, p - 2 * k - 8, px, rnd_mode)))
 100         break;
 101
 102       p += kmax;
 103       MPFR_ZIV_NEXT (loop, p);
 104       mpfr_set_prec (a, p);
 105       mpfr_set_prec (A, p);
 106       mpfr_set_prec (B, p);
 107       mpfr_set_prec (D, p);
 108       mpfr_set_prec (S, p);
 109   }
 110   MPFR_ZIV_FREE (loop);
 111   inex = mpfr_set (x, A, rnd_mode);
 112
 113   mpfr_clear (a);
 114   mpfr_clear (A);
 115   mpfr_clear (B);
 116   mpfr_clear (D);
 117   mpfr_clear (S);
 118
 119   return inex;
 120 }