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