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 } 467