contrib/mpfr/src/zeta.c

   1 /* mpfr_zeta -- compute the Riemann Zeta function
   2
   3 Copyright 2003, 2004, 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 /*
  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], MPFR_RNDN);
  56   for (l = 1; l < p; l++)
  57     {
  58       mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */
  59       mpfr_mul (d, d, s1, MPFR_RNDN);
  60       t = t - 1;
  61       mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */
  62       mpfr_mul (d, d, s1, MPFR_RNDN);
  63       t = t - 1;
  64       mpfr_div_ui (d, d, n2, MPFR_RNDN);
  65       mpfr_add (d, d, tc[p-l], MPFR_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], MPFR_RNDN);
  70     }
  71   mpfr_mul (d, d, s, MPFR_RNDN);
  72   mpfr_add (s1, s, __gmpfr_one, MPFR_RNDN);
  73   mpfr_neg (s1, s1, MPFR_RNDN);
  74   mpfr_ui_pow (u, n, s1, MPFR_RNDN);
  75   mpfr_mul (b, d, u, MPFR_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, MPFR_RNDN);
  94       for (k = 2; k <= p; k++)
  95         {
  96           mpfr_set_ui (d, k-1, MPFR_RNDN);
  97           mpfr_div_ui (d, d, 12*k+6, MPFR_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), MPFR_RNDN);
 101               mpfr_add (d, d, tc[l], MPFR_RNDN);
 102             }
 103           mpfr_div_ui (tc[k], d, 24, MPFR_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, MPFR_RNDN);
 123   mpfr_ui_pow (u, n, s1, MPFR_RNDN);
 124   mpfr_div_2ui (u, u, 1, MPFR_RNDN);
 125   mpfr_set (sum, u, MPFR_RNDN);
 126   for (i=n-1; i>1; i--)
 127     {
 128       mpfr_ui_pow (u, i, s1, MPFR_RNDN);
 129       mpfr_add (sum, sum, u, MPFR_RNDN);
 130     }
 131   mpfr_add (sum, sum, __gmpfr_one, MPFR_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, mpfr_rnd_t rnd_mode)
 143 {
 144   mpfr_t b, c, z_pre, f, s1;
 145   double beta, sd, dnep;
 146   mpfr_t *tc1;
 147   mpfr_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       mpfr_exp_t err;
 170       err = MPFR_GET_EXP (s) - 1;
 171       if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2))
 172         err = MPFR_EMAX_MAX;
 173       else
 174         err = ((mpfr_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, MPFR_RNDN);
 186   MPFR_ASSERTD (inex == 0);
 187
 188   /* case s=1 should have already been handled */
 189   MPFR_ASSERTD (!MPFR_IS_ZERO (s1));
 190
 191   MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
 192
 193   MPFR_ZIV_INIT (loop, d);
 194   for (;;)
 195     {
 196       /* Principal loop: we compute, in z_pre,
 197          an approximation of Zeta(s), that we send to can_round */
 198       if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2))
 199         /* Branch 1: when s-1 is very small, one
 200            uses the approximation Zeta(s)=1/(s-1)+gamma,
 201            where gamma is Euler's constant */
 202         {
 203           dint = MAX (d + 3, precs);
 204           MPFR_TRACE (printf ("branch 1\ninternal precision=%lu\n",
 205                               (unsigned long) dint));
 206           MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
 207           mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN);
 208           mpfr_const_euler (f, MPFR_RNDN);
 209           mpfr_add (z_pre, z_pre, f, MPFR_RNDN);
 210         }
 211       else /* Branch 2 */
 212         {
 213           size_t size;
 214
 215           MPFR_TRACE (printf ("branch 2\n"));
 216           /* Computation of parameters n, p and working precision */
 217           dnep = (double) d * LOG2;
 218           sd = mpfr_get_d (s, MPFR_RNDN);
 219           /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
 220              but a larger value is ok */
 221 #define LOG6dot2832 1.83787940484160805532
 222           beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
 223                                      __gmpfr_floor_log2 (sd));
 224           if (beta <= 0.0)
 225             {
 226               p = 0;
 227               /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
 228               n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
 229             }
 230           else
 231             {
 232               p = 1 + (int) beta / 2;
 233               n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
 234             }
 235           MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
 236           /* add = 4 + floor(1.5 * log(d) / log (2)).
 237              We should have add >= 10, which is always fulfilled since
 238              d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
 239           add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
 240           MPFR_ASSERTD(add >= 10);
 241           dint = d + add;
 242           if (dint < precs)
 243             dint = precs;
 244
 245           MPFR_TRACE (printf ("internal precision=%lu\n",
 246                               (unsigned long) dint));
 247
 248           size = (p + 1) * sizeof(mpfr_t);
 249           tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
 250           for (l=1; l<=p; l++)
 251             mpfr_init2 (tc1[l], dint);
 252           MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
 253
 254           MPFR_TRACE (printf ("precision of z = %lu\n",
 255                               (unsigned long) precz));
 256
 257           /* Computation of the coefficients c_k */
 258           mpfr_zeta_c (p, tc1);
 259           /* Computation of the 3 parts of the fonction Zeta. */
 260           mpfr_zeta_part_a (z_pre, s, n);
 261           mpfr_zeta_part_b (b, s, n, p, tc1);
 262           /* s1 = s-1 is already computed above */
 263           mpfr_div (c, __gmpfr_one, s1, MPFR_RNDN);
 264           mpfr_ui_pow (f, n, s1, MPFR_RNDN);
 265           mpfr_div (c, c, f, MPFR_RNDN);
 266           MPFR_TRACE (MPFR_DUMP (c));
 267           mpfr_add (z_pre, z_pre, c, MPFR_RNDN);
 268           mpfr_add (z_pre, z_pre, b, MPFR_RNDN);
 269           for (l=1; l<=p; l++)
 270             mpfr_clear (tc1[l]);
 271           (*__gmp_free_func) (tc1, size);
 272           /* End branch 2 */
 273         }
 274
 275       MPFR_TRACE (MPFR_DUMP (z_pre));
 276       if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
 277         break;
 278       MPFR_ZIV_NEXT (loop, d);
 279     }
 280   MPFR_ZIV_FREE (loop);
 281
 282   inex = mpfr_set (z, z_pre, rnd_mode);
 283
 284   MPFR_GROUP_CLEAR (group);
 285   mpfr_clear (s1);
 286
 287   return inex;
 288 }
 289
 290 int
 291 mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
 292 {
 293   mpfr_t z_pre, s1, y, p;
 294   double sd, eps, m1, c;
 295   long add;
 296   mpfr_prec_t precz, prec1, precs, precs1;
 297   int inex;
 298   MPFR_GROUP_DECL (group);
 299   MPFR_ZIV_DECL (loop);
 300   MPFR_SAVE_EXPO_DECL (expo);
 301
 302   MPFR_LOG_FUNC (
 303     ("s[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode),
 304     ("z[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex));
 305
 306   /* Zero, Nan or Inf ? */
 307   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
 308     {
 309       if (MPFR_IS_NAN (s))
 310         {
 311           MPFR_SET_NAN (z);
 312           MPFR_RET_NAN;
 313         }
 314       else if (MPFR_IS_INF (s))
 315         {
 316           if (MPFR_IS_POS (s))
 317             return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */
 318           MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
 319           MPFR_RET_NAN;
 320         }
 321       else /* s iz zero */
 322         {
 323           MPFR_ASSERTD (MPFR_IS_ZERO (s));
 324           return mpfr_set_si_2exp (z, -1, -1, rnd_mode);
 325         }
 326     }
 327
 328   /* s is neither Nan, nor Inf, nor Zero */
 329
 330   /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
 331      and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|.
 332      Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding
 333      (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest).
 334      A sufficient condition is that EXP(s) + 1 < -PREC(z). */
 335   if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z))
 336     {
 337       int signs = MPFR_SIGN(s);
 338
 339       MPFR_SAVE_EXPO_MARK (expo);
 340       mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
 341       if (rnd_mode == MPFR_RNDA)
 342         rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */
 343       if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0)
 344         {
 345           mpfr_nextabove (z); /* z = -1/2 + epsilon */
 346           inex = 1;
 347         }
 348       else if (rnd_mode == MPFR_RNDD && signs > 0)
 349         {
 350           mpfr_nextbelow (z); /* z = -1/2 - epsilon */
 351           inex = -1;
 352         }
 353       else
 354         {
 355           if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */
 356             inex = 1;
 357           else if (rnd_mode == MPFR_RNDD)
 358             inex = -1;              /* s < 0: z = -1/2 */
 359           else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */
 360             inex = (signs > 0) ? 1 : -1;
 361         }
 362       MPFR_SAVE_EXPO_FREE (expo);
 363       return mpfr_check_range (z, inex, rnd_mode);
 364     }
 365
 366   /* Check for case s= -2n */
 367   if (MPFR_IS_NEG (s))
 368     {
 369       mpfr_t tmp;
 370       tmp[0] = *s;
 371       MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1;
 372       if (mpfr_integer_p (tmp))
 373         {
 374           MPFR_SET_ZERO (z);
 375           MPFR_SET_POS (z);
 376           MPFR_RET (0);
 377         }
 378     }
 379
 380   /* Check for case s= 1 before changing the exponent range */
 381   if (mpfr_cmp (s, __gmpfr_one) ==0)
 382     {
 383       MPFR_SET_INF (z);
 384       MPFR_SET_POS (z);
 385       mpfr_set_divby0 ();
 386       MPFR_RET (0);
 387     }
 388
 389   MPFR_SAVE_EXPO_MARK (expo);
 390
 391   /* Compute Zeta */
 392   if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
 393     inex = mpfr_zeta_pos (z, s, rnd_mode);
 394   else /* use reflection formula
 395           zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
 396     {
 397       int overflow = 0;
 398
 399       precz = MPFR_PREC (z);
 400       precs = MPFR_PREC (s);
 401
 402       /* Precision precs1 needed to represent 1 - s, and s + 2,
 403          without any truncation */
 404       precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
 405       sd = mpfr_get_d (s, MPFR_RNDN) - 1.0;
 406       if (sd < 0.0)
 407         sd = -sd; /* now sd = abs(s-1.0) */
 408       /* Precision prec1 is the precision on elementary computations;
 409          it ensures a final precision prec1 - add for zeta(s) */
 410       /* eps = pow (2.0, - (double) precz - 14.0); */
 411       eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
 412       m1 = 1.0 + MAX(1.0 / eps,  2.0 * sd) * (1.0 + eps);
 413       c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
 414       /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
 415       add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
 416       prec1 = precz + add;
 417       prec1 = MAX (prec1, precs1) + 10;
 418
 419       MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
 420       MPFR_ZIV_INIT (loop, prec1);
 421       for (;;)
 422         {
 423           mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN);/* s1 = 1-s */
 424           mpfr_zeta_pos (z_pre, s1, MPFR_RNDN);   /* zeta(1-s)  */
 425           mpfr_gamma (y, s1, MPFR_RNDN);          /* gamma(1-s) */
 426           if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k,
 427                                   Zeta(s) > 0 for -4k < s < -4k+2 */
 428             {
 429               mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
 430               mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
 431               overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
 432               break;
 433             }
 434           mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);  /* gamma(1-s)*zeta(1-s) */
 435           mpfr_const_pi (p, MPFR_RNDD);
 436           mpfr_mul (y, s, p, MPFR_RNDN);
 437           mpfr_div_2ui (y, y, 1, MPFR_RNDN);      /* s*Pi/2 */
 438           mpfr_sin (y, y, MPFR_RNDN);             /* sin(Pi*s/2) */
 439           mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
 440           mpfr_mul_2ui (y, p, 1, MPFR_RNDN);      /* 2*Pi */
 441           mpfr_neg (s1, s1, MPFR_RNDN);           /* s-1 */
 442           mpfr_pow (y, y, s1, MPFR_RNDN);         /* (2*Pi)^(s-1) */
 443           mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
 444           mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);
 445
 446           if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
 447                                            rnd_mode)))
 448             break;
 449
 450           MPFR_ZIV_NEXT (loop, prec1);
 451           MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
 452         }
 453       MPFR_ZIV_FREE (loop);
 454       if (overflow != 0)
 455         {
 456           inex = mpfr_overflow (z, rnd_mode, overflow);
 457           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 458         }
 459       else
 460         inex = mpfr_set (z, z_pre, rnd_mode);
 461       MPFR_GROUP_CLEAR (group);
 462     }
 463
 464   MPFR_SAVE_EXPO_FREE (expo);
 465   return mpfr_check_range (z, inex, rnd_mode);
 466 }