contrib/mpfr/zeta.c

   1 /* mpfr_zeta -- compute the Riemann Zeta function
   2
   3 Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
   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 /*
  27    Parameters:
  28    s - the input floating-point number
  29    n, p - parameters from the algorithm
  30    tc - an array of p floating-point numbers tc[1]..tc[p]
  31    Output:
  32    b is the result, i.e.
  33    sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
  34 */
  35 static void
  36 mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
  37 {
  38   mpfr_t s1, d, u;
  39   unsigned long n2;
  40   int l, t;
  41   MPFR_GROUP_DECL (group);
  42
  43   if (p == 0)
  44     {
  45       MPFR_SET_ZERO (b);
  46       MPFR_SET_POS (b);
  47       return;
  48     }
  49
  50   n2 = n * n;
  51   MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u);
  52
  53   /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
  54   t = 2 * p - 2;
  55   mpfr_set (d, tc[p], GMP_RNDN);
  56   for (l = 1; l < p; l++)
  57     {
  58       mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l) */
  59       mpfr_mul (d, d, s1, GMP_RNDN);
  60       t = t - 1;
  61       mpfr_add_ui (s1, s, t, GMP_RNDN); /* s + (2p-2l-1) */
  62       mpfr_mul (d, d, s1, GMP_RNDN);
  63       t = t - 1;
  64       mpfr_div_ui (d, d, n2, GMP_RNDN);
  65       mpfr_add (d, d, tc[p-l], GMP_RNDN);
  66       /* since s is positive and the tc[i] have alternate signs,
  67          the following is unlikely */
  68       if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
  69         mpfr_set (d, tc[p-l], GMP_RNDN);
  70     }
  71   mpfr_mul (d, d, s, GMP_RNDN);
  72   mpfr_add (s1, s, __gmpfr_one, GMP_RNDN);
  73   mpfr_neg (s1, s1, GMP_RNDN);
  74   mpfr_ui_pow (u, n, s1, GMP_RNDN);
  75   mpfr_mul (b, d, u, GMP_RNDN);
  76
  77   MPFR_GROUP_CLEAR (group);
  78 }
  79
  80 /* Input: p - an integer
  81    Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
  82    tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
  83 */
  84 static void
  85 mpfr_zeta_c (int p, mpfr_t *tc)
  86 {
  87   mpfr_t d;
  88   int k, l;
  89
  90   if (p > 0)
  91     {
  92       mpfr_init2 (d, MPFR_PREC (tc[1]));
  93       mpfr_div_ui (tc[1], __gmpfr_one, 12, GMP_RNDN);
  94       for (k = 2; k <= p; k++)
  95         {
  96           mpfr_set_ui (d, k-1, GMP_RNDN);
  97           mpfr_div_ui (d, d, 12*k+6, GMP_RNDN);
  98           for (l=2; l < k; l++)
  99             {
 100               mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), GMP_RNDN);
 101               mpfr_add (d, d, tc[l], GMP_RNDN);
 102             }
 103           mpfr_div_ui (tc[k], d, 24, GMP_RNDN);
 104           MPFR_CHANGE_SIGN (tc[k]);
 105         }
 106       mpfr_clear (d);
 107     }
 108 }
 109
 110 /* Input: s - a floating-point number
 111           n - an integer
 112    Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
 113 static void
 114 mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
 115 {
 116   mpfr_t u, s1;
 117   int i;
 118   MPFR_GROUP_DECL (group);
 119
 120   MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);
 121
 122   mpfr_neg (s1, s, GMP_RNDN);
 123   mpfr_ui_pow (u, n, s1, GMP_RNDN);
 124   mpfr_div_2ui (u, u, 1, GMP_RNDN);
 125   mpfr_set (sum, u, GMP_RNDN);
 126   for (i=n-1; i>1; i--)
 127     {
 128       mpfr_ui_pow (u, i, s1, GMP_RNDN);
 129       mpfr_add (sum, sum, u, GMP_RNDN);
 130     }
 131   mpfr_add (sum, sum, __gmpfr_one, GMP_RNDN);
 132
 133   MPFR_GROUP_CLEAR (group);
 134 }
 135
 136 /* Input: s - a floating-point number >= 1/2.
 137           rnd_mode - a rounding mode.
 138           Assumes s is neither NaN nor Infinite.
 139    Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
 140 */
 141 static int
 142 mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
 143 {
 144   mpfr_t b, c, z_pre, f, s1;
 145   double beta, sd, dnep;
 146   mpfr_t *tc1;
 147   mp_prec_t precz, precs, d, dint;
 148   int p, n, l, add;
 149   int inex;
 150   MPFR_GROUP_DECL (group);
 151   MPFR_ZIV_DECL (loop);
 152
 153   MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);
 154
 155   precz = MPFR_PREC (z);
 156   precs = MPFR_PREC (s);
 157
 158   /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
 159      so with 2^(EXP(x)-1) <= x < 2^EXP(x)
 160      So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
 161      Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
 162              = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
 163             <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
 164      And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
 165      So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
 166      The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
 167   if (MPFR_GET_EXP (s) > 3)
 168     {
 169       mp_exp_t err;
 170       err = MPFR_GET_EXP (s) - 1;
 171       if (err > (mp_exp_t) (sizeof (mp_exp_t)*CHAR_BIT-2))
 172         err = MPFR_EMAX_MAX;
 173       else
 174         err = ((mp_exp_t)1) << err;
 175       err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
 176       MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1,
 177                                         rnd_mode, {});
 178     }
 179
 180   d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;
 181
 182   /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
 183   dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
 184   mpfr_init2 (s1, MAX (precs, dint));
 185   inex = mpfr_sub (s1, s, __gmpfr_one, GMP_RNDN);
 186   MPFR_ASSERTD (inex == 0);
 187
 188   /* case s=1 */
 189   if (MPFR_IS_ZERO (s1))
 190     {
 191       MPFR_SET_INF (z);
 192       MPFR_SET_POS (z);
 193       MPFR_ASSERTD (inex == 0);
 194       goto clear_and_return;
 195     }
 196
 197   MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
 198
 199   MPFR_ZIV_INIT (loop, d);
 200   for (;;)
 201     {
 202       /* Principal loop: we compute, in z_pre,
 203          an approximation of Zeta(s), that we send to can_round */
 204       if (MPFR_GET_EXP (s1) <= -(mp_exp_t) ((mpfr_prec_t) (d-3)/2))
 205         /* Branch 1: when s-1 is very small, one
 206            uses the approximation Zeta(s)=1/(s-1)+gamma,
 207            where gamma is Euler's constant */
 208         {
 209           dint = MAX (d + 3, precs);
 210           MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n",
 211                               (unsigned long) dint));
 212           MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
 213           mpfr_div (z_pre, __gmpfr_one, s1, GMP_RNDN);
 214           mpfr_const_euler (f, GMP_RNDN);
 215           mpfr_add (z_pre, z_pre, f, GMP_RNDN);
 216         }
 217       else /* Branch 2 */
 218         {
 219           size_t size;
 220
 221           MPFR_TRACE (printf ("branch 2\n"));
 222           /* Computation of parameters n, p and working precision */
 223           dnep = (double) d * LOG2;
 224           sd = mpfr_get_d (s, GMP_RNDN);
 225           /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
 226              but a larger value is ok */
 227 #define LOG6dot2832 1.83787940484160805532
 228           beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
 229                                      __gmpfr_floor_log2 (sd));
 230           if (beta <= 0.0)
 231             {
 232               p = 0;
 233               /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
 234               n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
 235             }
 236           else
 237             {
 238               p = 1 + (int) beta / 2;
 239               n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
 240             }
 241           MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
 242           /* add = 4 + floor(1.5 * log(d) / log (2)).
 243              We should have add >= 10, which is always fulfilled since
 244              d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
 245           add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
 246           MPFR_ASSERTD(add >= 10);
 247           dint = d + add;
 248           if (dint < precs)
 249             dint = precs;
 250
 251           MPFR_TRACE (printf ("internal precision=%lu\n",
 252                               (unsigned long) dint));
 253
 254           size = (p + 1) * sizeof(mpfr_t);
 255           tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
 256           for (l=1; l<=p; l++)
 257             mpfr_init2 (tc1[l], dint);
 258           MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
 259
 260           MPFR_TRACE (printf ("precision of z = %lu\n",
 261                               (unsigned long) precz));
 262
 263           /* Computation of the coefficients c_k */
 264           mpfr_zeta_c (p, tc1);
 265           /* Computation of the 3 parts of the fonction Zeta. */
 266           mpfr_zeta_part_a (z_pre, s, n);
 267           mpfr_zeta_part_b (b, s, n, p, tc1);
 268           /* s1 = s-1 is already computed above */
 269           mpfr_div (c, __gmpfr_one, s1, GMP_RNDN);
 270           mpfr_ui_pow (f, n, s1, GMP_RNDN);
 271           mpfr_div (c, c, f, GMP_RNDN);
 272           MPFR_TRACE (MPFR_DUMP (c));
 273           mpfr_add (z_pre, z_pre, c, GMP_RNDN);
 274           mpfr_add (z_pre, z_pre, b, GMP_RNDN);
 275           for (l=1; l<=p; l++)
 276             mpfr_clear (tc1[l]);
 277           (*__gmp_free_func) (tc1, size);
 278           /* End branch 2 */
 279         }
 280
 281       MPFR_TRACE (MPFR_DUMP (z_pre));
 282       if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
 283         break;
 284       MPFR_ZIV_NEXT (loop, d);
 285     }
 286   MPFR_ZIV_FREE (loop);
 287
 288   inex = mpfr_set (z, z_pre, rnd_mode);
 289
 290   MPFR_GROUP_CLEAR (group);
 291  clear_and_return:
 292   mpfr_clear (s1);
 293
 294   return inex;
 295 }
 296
 297 int
 298 mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
 299 {
 300   mpfr_t z_pre, s1, y, p;
 301   double sd, eps, m1, c;
 302   long add;
 303   mp_prec_t precz, prec1, precs, precs1;
 304   int inex;
 305   MPFR_GROUP_DECL (group);
 306   MPFR_ZIV_DECL (loop);
 307   MPFR_SAVE_EXPO_DECL (expo);
 308
 309   MPFR_LOG_FUNC (("s[%#R]=%R rnd=%d", s, s, rnd_mode),
 310                  ("z[%#R]=%R inexact=%d", z, z, inex));
 311
 312   /* Zero, Nan or Inf ? */
 313   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
 314     {
 315       if (MPFR_IS_NAN (s))
 316         {
 317           MPFR_SET_NAN (z);
 318           MPFR_RET_NAN;
 319         }
 320       else if (MPFR_IS_INF (s))
 321         {
 322           if (MPFR_IS_POS (s))
 323             return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */
 324           MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
 325           MPFR_RET_NAN;
 326         }
 327       else /* s iz zero */
 328         {
 329           MPFR_ASSERTD (MPFR_IS_ZERO (s));
 330           mpfr_set_ui (z, 1, rnd_mode);
 331           mpfr_div_2ui (z, z, 1, rnd_mode);
 332           MPFR_CHANGE_SIGN (z);
 333           MPFR_RET (0);
 334         }
 335     }
 336
 337   /* s is neither Nan, nor Inf, nor Zero */
 338
 339   /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
 340      and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|.
 341      Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding
 342      (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest).
 343      A sufficient condition is that EXP(s) + 1 < -PREC(z). */
 344   if (MPFR_EXP(s) + 1 < - (mp_exp_t) MPFR_PREC(z))
 345     {
 346       int signs = MPFR_SIGN(s);
 347       mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
 348       if ((rnd_mode == GMP_RNDU || rnd_mode == GMP_RNDZ) && signs < 0)
 349         {
 350           mpfr_nextabove (z); /* z = -1/2 + epsilon */
 351           inex = 1;
 352         }
 353       else if (rnd_mode == GMP_RNDD && signs > 0)
 354         {
 355           mpfr_nextbelow (z); /* z = -1/2 - epsilon */
 356           inex = -1;
 357         }
 358       else
 359         {
 360           if (rnd_mode == GMP_RNDU) /* s > 0: z = -1/2 */
 361             inex = 1;
 362           else if (rnd_mode == GMP_RNDD)
 363             inex = -1;              /* s < 0: z = -1/2 */
 364           else /* (GMP_RNDZ and s > 0) or GMP_RNDN: z = -1/2 */
 365             inex = (signs > 0) ? 1 : -1;
 366         }
 367       return mpfr_check_range (z, inex, rnd_mode);
 368     }
 369
 370   /* Check for case s= -2n */
 371   if (MPFR_IS_NEG (s))
 372     {
 373       mpfr_t tmp;
 374       tmp[0] = *s;
 375       MPFR_EXP (tmp) = MPFR_EXP (s) - 1;
 376       if (mpfr_integer_p (tmp))
 377         {
 378           MPFR_SET_ZERO (z);
 379           MPFR_SET_POS (z);
 380           MPFR_RET (0);
 381         }
 382     }
 383
 384   MPFR_SAVE_EXPO_MARK (expo);
 385
 386   /* Compute Zeta */
 387   if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
 388     inex = mpfr_zeta_pos (z, s, rnd_mode);
 389   else /* use reflection formula
 390           zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
 391     {
 392       int overflow = 0;
 393
 394       precz = MPFR_PREC (z);
 395       precs = MPFR_PREC (s);
 396
 397       /* Precision precs1 needed to represent 1 - s, and s + 2,
 398          without any truncation */
 399       precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
 400       sd = mpfr_get_d (s, GMP_RNDN) - 1.0;
 401       if (sd < 0.0)
 402         sd = -sd; /* now sd = abs(s-1.0) */
 403       /* Precision prec1 is the precision on elementary computations;
 404          it ensures a final precision prec1 - add for zeta(s) */
 405       /* eps = pow (2.0, - (double) precz - 14.0); */
 406       eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
 407       m1 = 1.0 + MAX(1.0 / eps,  2.0 * sd) * (1.0 + eps);
 408       c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
 409       /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
 410       add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
 411       prec1 = precz + add;
 412       prec1 = MAX (prec1, precs1) + 10;
 413
 414       MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
 415       MPFR_ZIV_INIT (loop, prec1);
 416       for (;;)
 417         {
 418           mpfr_sub (s1, __gmpfr_one, s, GMP_RNDN);/* s1 = 1-s */
 419           mpfr_zeta_pos (z_pre, s1, GMP_RNDN);   /* zeta(1-s)  */
 420           mpfr_gamma (y, s1, GMP_RNDN);          /* gamma(1-s) */
 421           if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
 422                                   Zeta(s) > 0 for -4k < s < -4k+2 */
 423             {
 424               mpfr_div_2ui (s1, s, 2, GMP_RNDN); /* s/4, exact */
 425               mpfr_frac (s1, s1, GMP_RNDN); /* exact, -1 < s1 < 0 */
 426               overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
 427               break;
 428             }
 429           mpfr_mul (z_pre, z_pre, y, GMP_RNDN);  /* gamma(1-s)*zeta(1-s) */
 430           mpfr_const_pi (p, GMP_RNDD);
 431           mpfr_mul (y, s, p, GMP_RNDN);
 432           mpfr_div_2ui (y, y, 1, GMP_RNDN);      /* s*Pi/2 */
 433           mpfr_sin (y, y, GMP_RNDN);             /* sin(Pi*s/2) */
 434           mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
 435           mpfr_mul_2ui (y, p, 1, GMP_RNDN);      /* 2*Pi */
 436           mpfr_neg (s1, s1, GMP_RNDN);           /* s-1 */
 437           mpfr_pow (y, y, s1, GMP_RNDN);         /* (2*Pi)^(s-1) */
 438           mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
 439           mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN);
 440
 441           if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
 442                                            rnd_mode)))
 443             break;
 444
 445           MPFR_ZIV_NEXT (loop, prec1);
 446           MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
 447         }
 448       MPFR_ZIV_FREE (loop);
 449       if (overflow != 0)
 450         {
 451           inex = mpfr_overflow (z, rnd_mode, overflow);
 452           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 453         }
 454       else
 455         inex = mpfr_set (z, z_pre, rnd_mode);
 456       MPFR_GROUP_CLEAR (group);
 457     }
 458
 459   MPFR_SAVE_EXPO_FREE (expo);
 460   return mpfr_check_range (z, inex, rnd_mode);
 461 }