contrib/mpfr/zeta_ui.c

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