source/libs/mpfr/mpfr-3.1.3/src/zeta_ui.c

   1 /* mpfr_zeta_ui -- compute the Riemann Zeta function for integer argument.
   2
   3 Copyright 2005-2015 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 int
  27 mpfr_zeta_ui (mpfr_ptr z, unsigned long m, mpfr_rnd_t r)
  28 {
  29   MPFR_ZIV_DECL (loop);
  30
  31   MPFR_LOG_FUNC
  32     (("m=%lu rnd=%d prec=%Pu", m, r, mpfr_get_prec (z)),
  33      ("z[%Pu]=%.*Rg", mpfr_get_prec (z), mpfr_log_prec, z));
  34
  35   if (m == 0)
  36     {
  37       return mpfr_set_si_2exp (z, -1, -1, r);
  38     }
  39   else if (m == 1)
  40     {
  41       MPFR_SET_INF (z);
  42       MPFR_SET_POS (z);
  43       mpfr_set_divby0 ();
  44       return 0;
  45     }
  46   else /* m >= 2 */
  47     {
  48       mpfr_prec_t p = MPFR_PREC(z);
  49       unsigned long n, k, err, kbits;
  50       mpz_t d, t, s, q;
  51       mpfr_t y;
  52       int inex;
  53       MPFR_SAVE_EXPO_DECL (expo);
  54
  55       if (r == MPFR_RNDA)
  56         r = MPFR_RNDU; /* since the result is always positive */
  57
  58       MPFR_SAVE_EXPO_MARK (expo);
  59
  60       if (m >= p) /* 2^(-m) < ulp(1) = 2^(1-p). This means that
  61                      2^(-m) <= 1/2*ulp(1). We have 3^(-m)+4^(-m)+... < 2^(-m)
  62                      i.e. zeta(m) < 1+2*2^(-m) for m >= 3 */
  63         {
  64           if (m == 2) /* necessarily p=2 */
  65             inex = mpfr_set_ui_2exp (z, 13, -3, r);
  66           else if (r == MPFR_RNDZ || r == MPFR_RNDD ||
  67                    (r == MPFR_RNDN && m > p))
  68             {
  69               mpfr_set_ui (z, 1, r);
  70               inex = -1;
  71             }
  72           else
  73             {
  74               mpfr_set_ui (z, 1, r);
  75               mpfr_nextabove (z);
  76               inex = 1;
  77             }
  78           goto end;
  79         }
  80
  81       /* now treat also the case where zeta(m) - (1+1/2^m) < 1/2*ulp(1),
  82          and the result is either 1+2^(-m) or 1+2^(-m)+2^(1-p). */
  83       mpfr_init2 (y, 31);
  84
  85       if (m >= p / 2) /* otherwise 4^(-m) > 2^(-p) */
  86         {
  87           /* the following is a lower bound for log(3)/log(2) */
  88           mpfr_set_str_binary (y, "1.100101011100000000011010001110");
  89           mpfr_mul_ui (y, y, m, MPFR_RNDZ); /* lower bound for log2(3^m) */
  90           if (mpfr_cmp_ui (y, p + 2) >= 0)
  91             {
  92               mpfr_clear (y);
  93               mpfr_set_ui (z, 1, MPFR_RNDZ);
  94               mpfr_div_2ui (z, z, m, MPFR_RNDZ);
  95               mpfr_add_ui (z, z, 1, MPFR_RNDZ);
  96               if (r != MPFR_RNDU)
  97                 inex = -1;
  98               else
  99                 {
 100                   mpfr_nextabove (z);
 101                   inex = 1;
 102                 }
 103               goto end;
 104             }
 105         }
 106
 107       mpz_init (s);
 108       mpz_init (d);
 109       mpz_init (t);
 110       mpz_init (q);
 111
 112       p += MPFR_INT_CEIL_LOG2(p); /* account of the n term in the error */
 113
 114       p += MPFR_INT_CEIL_LOG2(p) + 15; /* initial value */
 115
 116       MPFR_ZIV_INIT (loop, p);
 117       for(;;)
 118         {
 119           /* 0.39321985067869744 = log(2)/log(3+sqrt(8)) */
 120           n = 1 + (unsigned long) (0.39321985067869744 * (double) p);
 121           err = n + 4;
 122
 123           mpfr_set_prec (y, p);
 124
 125           /* computation of the d[k] */
 126           mpz_set_ui (s, 0);
 127           mpz_set_ui (t, 1);
 128           mpz_mul_2exp (t, t, 2 * n - 1); /* t[n] */
 129           mpz_set (d, t);
 130           for (k = n; k > 0; k--)
 131             {
 132               count_leading_zeros (kbits, k);
 133               kbits = GMP_NUMB_BITS - kbits;
 134               /* if k^m is too large, use mpz_tdiv_q */
 135               if (m * kbits > 2 * GMP_NUMB_BITS)
 136                 {
 137                   /* if we know in advance that k^m > d, then floor(d/k^m) will
 138                      be zero below, so there is no need to compute k^m */
 139                   kbits = (kbits - 1) * m + 1;
 140                   /* k^m has at least kbits bits */
 141                   if (kbits > mpz_sizeinbase (d, 2))
 142                     mpz_set_ui (q, 0);
 143                   else
 144                     {
 145                       mpz_ui_pow_ui (q, k, m);
 146                       mpz_tdiv_q (q, d, q);
 147                     }
 148                 }
 149               else /* use several mpz_tdiv_q_ui calls */
 150                 {
 151                   unsigned long km = k, mm = m - 1;
 152                   while (mm > 0 && km < ULONG_MAX / k)
 153                     {
 154                       km *= k;
 155                       mm --;
 156                     }
 157                   mpz_tdiv_q_ui (q, d, km);
 158                   while (mm > 0)
 159                     {
 160                       km = k;
 161                       mm --;
 162                       while (mm > 0 && km < ULONG_MAX / k)
 163                         {
 164                           km *= k;
 165                           mm --;
 166                         }
 167                       mpz_tdiv_q_ui (q, q, km);
 168                     }
 169                 }
 170               if (k % 2)
 171                 mpz_add (s, s, q);
 172               else
 173                 mpz_sub (s, s, q);
 174
 175               /* we have d[k] = sum(t[i], i=k+1..n)
 176                  with t[i] = n*(n+i-1)!*4^i/(n-i)!/(2i)!
 177                  t[k-1]/t[k] = k*(2k-1)/(n-k+1)/(n+k-1)/2 */
 178 #if (GMP_NUMB_BITS == 32)
 179 #define KMAX 46341 /* max k such that k*(2k-1) < 2^32 */
 180 #elif (GMP_NUMB_BITS == 64)
 181 #define KMAX 3037000500
 182 #endif
 183 #ifdef KMAX
 184               if (k <= KMAX)
 185                 mpz_mul_ui (t, t, k * (2 * k - 1));
 186               else
 187 #endif
 188                 {
 189                   mpz_mul_ui (t, t, k);
 190                   mpz_mul_ui (t, t, 2 * k - 1);
 191                 }
 192               mpz_fdiv_q_2exp (t, t, 1);
 193               /* Warning: the test below assumes that an unsigned long
 194                  has no padding bits. */
 195               if (n < 1UL << ((sizeof(unsigned long) * CHAR_BIT) / 2))
 196                 /* (n - k + 1) * (n + k - 1) < n^2 */
 197                 mpz_divexact_ui (t, t, (n - k + 1) * (n + k - 1));
 198               else
 199                 {
 200                   mpz_divexact_ui (t, t, n - k + 1);
 201                   mpz_divexact_ui (t, t, n + k - 1);
 202                 }
 203               mpz_add (d, d, t);
 204             }
 205
 206           /* multiply by 1/(1-2^(1-m)) = 1 + 2^(1-m) + 2^(2-m) + ... */
 207           mpz_fdiv_q_2exp (t, s, m - 1);
 208           do
 209             {
 210               err ++;
 211               mpz_add (s, s, t);
 212               mpz_fdiv_q_2exp (t, t, m - 1);
 213             }
 214           while (mpz_cmp_ui (t, 0) > 0);
 215
 216           /* divide by d[n] */
 217           mpz_mul_2exp (s, s, p);
 218           mpz_tdiv_q (s, s, d);
 219           mpfr_set_z (y, s, MPFR_RNDN);
 220           mpfr_div_2ui (y, y, p, MPFR_RNDN);
 221
 222           err = MPFR_INT_CEIL_LOG2 (err);
 223
 224           if (MPFR_LIKELY(MPFR_CAN_ROUND (y, p - err, MPFR_PREC(z), r)))
 225             break;
 226
 227           MPFR_ZIV_NEXT (loop, p);
 228         }
 229       MPFR_ZIV_FREE (loop);
 230
 231       mpz_clear (d);
 232       mpz_clear (t);
 233       mpz_clear (q);
 234       mpz_clear (s);
 235       inex = mpfr_set (z, y, r);
 236       mpfr_clear (y);
 237
 238     end:
 239       MPFR_LOG_VAR (z);
 240       MPFR_LOG_MSG (("inex = %d before mpfr_check_range\n", inex));
 241       MPFR_SAVE_EXPO_FREE (expo);
 242       return mpfr_check_range (z, inex, r);
 243     }
 244 }