contrib/mpfr/src/zeta_ui.c

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