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