X-Git-Url: https://oss.titaniummirror.com/gitweb/?a=blobdiff_plain;f=mpfr%2Fconst_pi.c;fp=mpfr%2Fconst_pi.c;h=19599c009ef0d48b410ef170225402385ab9bc4e;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git

diff --git a/mpfr/const_pi.c b/mpfr/const_pi.c
new file mode 100644
index 00000000..19599c00
--- /dev/null
+++ b/mpfr/const_pi.c
@@ -0,0 +1,120 @@
+/* mpfr_const_pi -- compute Pi
+
+Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+/* Declare the cache */
+MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_pi, mpfr_const_pi_internal);
+
+/* Set User Interface */
+#undef mpfr_const_pi
+int
+mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode) {
+  return mpfr_cache (x, __gmpfr_cache_const_pi, rnd_mode);
+}
+
+/* Don't need to save/restore exponent range: the cache does it */
+int
+mpfr_const_pi_internal (mpfr_ptr x, mp_rnd_t rnd_mode)
+{
+  mpfr_t a, A, B, D, S;
+  mp_prec_t px, p, cancel, k, kmax;
+  MPFR_ZIV_DECL (loop);
+  int inex;
+
+  MPFR_LOG_FUNC (("rnd_mode=%d", rnd_mode), ("x[%#R]=%R inex=%d", x, x, inex));
+
+  px = MPFR_PREC (x);
+
+  /* we need 9*2^kmax - 4 >= px+2*kmax+8 */
+  for (kmax = 2; ((px + 2 * kmax + 12) / 9) >> kmax; kmax ++);
+
+  p = px + 3 * kmax + 14; /* guarantees no recomputation for px <= 10000 */
+
+  mpfr_init2 (a, p);
+  mpfr_init2 (A, p);
+  mpfr_init2 (B, p);
+  mpfr_init2 (D, p);
+  mpfr_init2 (S, p);
+
+  MPFR_ZIV_INIT (loop, p);
+  for (;;) {
+    mpfr_set_ui (a, 1, GMP_RNDN);          /* a = 1 */
+    mpfr_set_ui (A, 1, GMP_RNDN);          /* A = a^2 = 1 */
+    mpfr_set_ui_2exp (B, 1, -1, GMP_RNDN); /* B = b^2 = 1/2 */
+    mpfr_set_ui_2exp (D, 1, -2, GMP_RNDN); /* D = 1/4 */
+
+#define b B
+#define ap a
+#define Ap A
+#define Bp B
+    for (k = 0, cancel = 0; ; k++)
+      {
+        /* invariant: 1/2 <= B <= A <= a < 1 */
+        mpfr_add (S, A, B, GMP_RNDN); /* 1 <= S <= 2 */
+        mpfr_div_2ui (S, S, 2, GMP_RNDN); /* exact, 1/4 <= S <= 1/2 */
+        mpfr_sqrt (b, B, GMP_RNDN); /* 1/2 <= b <= 1 */
+        mpfr_add (ap, a, b, GMP_RNDN); /* 1 <= ap <= 2 */
+        mpfr_div_2ui (ap, ap, 1, GMP_RNDN); /* exact, 1/2 <= ap <= 1 */
+        mpfr_mul (Ap, ap, ap, GMP_RNDN); /* 1/4 <= Ap <= 1 */
+        mpfr_sub (Bp, Ap, S, GMP_RNDN); /* -1/4 <= Bp <= 3/4 */
+        mpfr_mul_2ui (Bp, Bp, 1, GMP_RNDN); /* -1/2 <= Bp <= 3/2 */
+        mpfr_sub (S, Ap, Bp, GMP_RNDN);
+        MPFR_ASSERTN (mpfr_cmp_ui (S, 1) < 0);
+        cancel = mpfr_cmp_ui (S, 0) ? (mpfr_uexp_t) -mpfr_get_exp(S) : p;
+        /* MPFR_ASSERTN (cancel >= px || cancel >= 9 * (1 << k) - 4); */
+        mpfr_mul_2ui (S, S, k, GMP_RNDN);
+        mpfr_sub (D, D, S, GMP_RNDN);
+        /* stop when |A_k - B_k| <= 2^(k-p) i.e. cancel >= p-k */
+        if (cancel + k >= p)
+          break;
+      }
+#undef b
+#undef ap
+#undef Ap
+#undef Bp
+
+      mpfr_div (A, B, D, GMP_RNDN);
+
+      /* MPFR_ASSERTN(p >= 2 * k + 8); */
+      if (MPFR_LIKELY (MPFR_CAN_ROUND (A, p - 2 * k - 8, px, rnd_mode)))
+        break;
+
+      p += kmax;
+      MPFR_ZIV_NEXT (loop, p);
+      mpfr_set_prec (a, p);
+      mpfr_set_prec (A, p);
+      mpfr_set_prec (B, p);
+      mpfr_set_prec (D, p);
+      mpfr_set_prec (S, p);
+  }
+  MPFR_ZIV_FREE (loop);
+  inex = mpfr_set (x, A, rnd_mode);
+
+  mpfr_clear (a);
+  mpfr_clear (A);
+  mpfr_clear (B);
+  mpfr_clear (D);
+  mpfr_clear (S);
+
+  return inex;
+}