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               if (n < 1UL << (BITS_PER_MP_LIMB / 2))
 181                 /* (n - k + 1) * (n + k - 1) < n^2 */
 182                 mpz_divexact_ui (t, t, (n - k + 1) * (n + k - 1));
 183               else
 184                 {
 185                   mpz_divexact_ui (t, t, n - k + 1);
 186                   mpz_divexact_ui (t, t, n + k - 1);
 187                 }
 188               mpz_add (d, d, t);
 189             }
 190
 191           /* multiply by 1/(1-2^(1-m)) = 1 + 2^(1-m) + 2^(2-m) + ... */
 192           mpz_div_2exp (t, s, m - 1);
 193           do
 194             {
 195               err ++;
 196               mpz_add (s, s, t);
 197               mpz_div_2exp (t, t, m - 1);
 198             }
 199           while (mpz_cmp_ui (t, 0) > 0);
 200
 201           /* divide by d[n] */
 202           mpz_mul_2exp (s, s, p);
 203           mpz_tdiv_q (s, s, d);
 204           mpfr_set_z (y, s, GMP_RNDN);
 205           mpfr_div_2ui (y, y, p, GMP_RNDN);
 206
 207           err = MPFR_INT_CEIL_LOG2 (err);
 208
 209           if (MPFR_LIKELY(MPFR_CAN_ROUND (y, p - err, MPFR_PREC(z), r)))
 210             break;
 211
 212           MPFR_ZIV_NEXT (loop, p);
 213         }
 214       MPFR_ZIV_FREE (loop);
 215
 216       mpz_clear (d);
 217       mpz_clear (t);
 218       mpz_clear (q);
 219       mpz_clear (s);
 220       inex = mpfr_set (z, y, r);
 221       mpfr_clear (y);
 222       return inex;
 223     }
 224 }