source/libs/mpfr/mpfr-src/src/const_catalan.c

   1 /* mpfr_const_catalan -- compute Catalan's constant.
   2
   3 Copyright 2005-2016 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramba 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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* Declare the cache */
  27 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_catalan, mpfr_const_catalan_internal);
  28
  29 #ifdef MPFR_WIN_THREAD_SAFE_DLL
  30 mpfr_cache_t *
  31 __gmpfr_cache_const_catalan_f()
  32 {
  33   return &__gmpfr_cache_const_catalan;
  34 }
  35 #endif
  36
  37 /* Set User Interface */
  38 #undef mpfr_const_catalan
  39 int
  40 mpfr_const_catalan (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
  41   return mpfr_cache (x, __gmpfr_cache_const_catalan, rnd_mode);
  42 }
  43
  44 /* return T, Q such that T/Q = sum(k!^2/(2k)!/(2k+1)^2, k=n1..n2-1) */
  45 static void
  46 S (mpz_t T, mpz_t P, mpz_t Q, unsigned long n1, unsigned long n2)
  47 {
  48   if (n2 == n1 + 1)
  49     {
  50       if (n1 == 0)
  51         {
  52           mpz_set_ui (P, 1);
  53           mpz_set_ui (Q, 1);
  54         }
  55       else
  56         {
  57           mpz_set_ui (P, 2 * n1 - 1);
  58           mpz_mul_ui (P, P, n1);
  59           mpz_ui_pow_ui (Q, 2 * n1 + 1, 2);
  60           mpz_mul_2exp (Q, Q, 1);
  61         }
  62       mpz_set (T, P);
  63     }
  64   else
  65     {
  66       unsigned long m = (n1 + n2) / 2;
  67       mpz_t T2, P2, Q2;
  68       S (T, P, Q, n1, m);
  69       mpz_init (T2);
  70       mpz_init (P2);
  71       mpz_init (Q2);
  72       S (T2, P2, Q2, m, n2);
  73       mpz_mul (T, T, Q2);
  74       mpz_mul (T2, T2, P);
  75       mpz_add (T, T, T2);
  76       mpz_mul (P, P, P2);
  77       mpz_mul (Q, Q, Q2);
  78       mpz_clear (T2);
  79       mpz_clear (P2);
  80       mpz_clear (Q2);
  81     }
  82 }
  83
  84 /* Don't need to save/restore exponent range: the cache does it.
  85    Catalan's constant is G = sum((-1)^k/(2*k+1)^2, k=0..infinity).
  86    We compute it using formula (31) of Victor Adamchik's page
  87    "33 representations for Catalan's constant"
  88    http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm
  89
  90    G = Pi/8*log(2+sqrt(3)) + 3/8*sum(k!^2/(2k)!/(2k+1)^2,k=0..infinity)
  91 */
  92 int
  93 mpfr_const_catalan_internal (mpfr_ptr g, mpfr_rnd_t rnd_mode)
  94 {
  95   mpfr_t x, y, z;
  96   mpz_t T, P, Q;
  97   mpfr_prec_t pg, p;
  98   int inex;
  99   MPFR_ZIV_DECL (loop);
 100   MPFR_GROUP_DECL (group);
 101
 102   MPFR_LOG_FUNC (("rnd_mode=%d", rnd_mode),
 103     ("g[%Pu]=%.*Rg inex=%d", mpfr_get_prec (g), mpfr_log_prec, g, inex));
 104
 105   /* Here are the WC (max prec = 100.000.000)
 106      Once we have found a chain of 11, we only look for bigger chain.
 107      Found 3 '1' at 0
 108      Found 5 '1' at 9
 109      Found 6 '0' at 34
 110      Found 9 '1' at 176
 111      Found 11 '1' at 705
 112      Found 12 '0' at 913
 113      Found 14 '1' at 12762
 114      Found 15 '1' at 152561
 115      Found 16 '0' at 171725
 116      Found 18 '0' at 525355
 117      Found 20 '0' at 529245
 118      Found 21 '1' at 6390133
 119      Found 22 '0' at 7806417
 120      Found 25 '1' at 11936239
 121      Found 27 '1' at 51752950
 122   */
 123   pg = MPFR_PREC (g);
 124   p = pg + MPFR_INT_CEIL_LOG2 (pg) + 7;
 125
 126   MPFR_GROUP_INIT_3 (group, p, x, y, z);
 127   mpz_init (T);
 128   mpz_init (P);
 129   mpz_init (Q);
 130
 131   MPFR_ZIV_INIT (loop, p);
 132   for (;;) {
 133     mpfr_sqrt_ui (x, 3, MPFR_RNDU);
 134     mpfr_add_ui (x, x, 2, MPFR_RNDU);
 135     mpfr_log (x, x, MPFR_RNDU);
 136     mpfr_const_pi (y, MPFR_RNDU);
 137     mpfr_mul (x, x, y, MPFR_RNDN);
 138     S (T, P, Q, 0, (p - 1) / 2);
 139     mpz_mul_ui (T, T, 3);
 140     mpfr_set_z (y, T, MPFR_RNDU);
 141     mpfr_set_z (z, Q, MPFR_RNDD);
 142     mpfr_div (y, y, z, MPFR_RNDN);
 143     mpfr_add (x, x, y, MPFR_RNDN);
 144     mpfr_div_2ui (x, x, 3, MPFR_RNDN);
 145
 146     if (MPFR_LIKELY (MPFR_CAN_ROUND (x, p - 5, pg, rnd_mode)))
 147       break;
 148
 149     MPFR_ZIV_NEXT (loop, p);
 150     MPFR_GROUP_REPREC_3 (group, p, x, y, z);
 151   }
 152   MPFR_ZIV_FREE (loop);
 153   inex = mpfr_set (g, x, rnd_mode);
 154
 155   MPFR_GROUP_CLEAR (group);
 156   mpz_clear (T);
 157   mpz_clear (P);
 158   mpz_clear (Q);
 159
 160   return inex;
 161 }