1/* mpfr_zeta -- compute the Riemann Zeta function 2 3Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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*/ 35static void 36mpfr_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*/ 84static void 85mpfr_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) */ 113static void 114mpfr_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*/ 141static int 142mpfr_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 */ 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) <= -(mpfr_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, MPFR_RNDN); 214 mpfr_const_euler (f, MPFR_RNDN); 215 mpfr_add (z_pre, z_pre, f, MPFR_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, MPFR_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, MPFR_RNDN); 270 mpfr_ui_pow (f, n, s1, MPFR_RNDN); 271 mpfr_div (c, c, f, MPFR_RNDN); 272 MPFR_TRACE (MPFR_DUMP (c)); 273 mpfr_add (z_pre, z_pre, c, MPFR_RNDN); 274 mpfr_add (z_pre, z_pre, b, MPFR_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 297int 298mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) 299{ 300 mpfr_t z_pre, s1, y, p; 301 double sd, eps, m1, c; 302 long add; 303 mpfr_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, MPFR_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 return mpfr_set_si_2exp (z, -1, -1, rnd_mode); 331 } 332 } 333 334 /* s is neither Nan, nor Inf, nor Zero */ 335 336 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0, 337 and for |s| <= 0.074, we have |zeta(s) + 1/2| <= |s|. 338 Thus if |s| <= 1/4*ulp(1/2), we can deduce the correct rounding 339 (the 1/4 covers the case where |zeta(s)| < 1/2 and rounding to nearest). 340 A sufficient condition is that EXP(s) + 1 < -PREC(z). */ 341 if (MPFR_EXP(s) + 1 < - (mpfr_exp_t) MPFR_PREC(z)) 342 { 343 int signs = MPFR_SIGN(s); 344 345 MPFR_SAVE_EXPO_MARK (expo); 346 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */ 347 if (rnd_mode == MPFR_RNDA) 348 rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */ 349 if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0) 350 { 351 mpfr_nextabove (z); /* z = -1/2 + epsilon */ 352 inex = 1; 353 } 354 else if (rnd_mode == MPFR_RNDD && signs > 0) 355 { 356 mpfr_nextbelow (z); /* z = -1/2 - epsilon */ 357 inex = -1; 358 } 359 else 360 { 361 if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */ 362 inex = 1; 363 else if (rnd_mode == MPFR_RNDD) 364 inex = -1; /* s < 0: z = -1/2 */ 365 else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */ 366 inex = (signs > 0) ? 1 : -1; 367 } 368 MPFR_SAVE_EXPO_FREE (expo); 369 return mpfr_check_range (z, inex, rnd_mode); 370 } 371 372 /* Check for case s= -2n */ 373 if (MPFR_IS_NEG (s)) 374 { 375 mpfr_t tmp; 376 tmp[0] = *s; 377 MPFR_EXP (tmp) = MPFR_EXP (s) - 1; 378 if (mpfr_integer_p (tmp)) 379 { 380 MPFR_SET_ZERO (z); 381 MPFR_SET_POS (z); 382 MPFR_RET (0); 383 } 384 } 385 386 MPFR_SAVE_EXPO_MARK (expo); 387 388 /* Compute Zeta */ 389 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ 390 inex = mpfr_zeta_pos (z, s, rnd_mode); 391 else /* use reflection formula 392 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */ 393 { 394 int overflow = 0; 395 396 precz = MPFR_PREC (z); 397 precs = MPFR_PREC (s); 398 399 /* Precision precs1 needed to represent 1 - s, and s + 2, 400 without any truncation */ 401 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); 402 sd = mpfr_get_d (s, MPFR_RNDN) - 1.0; 403 if (sd < 0.0) 404 sd = -sd; /* now sd = abs(s-1.0) */ 405 /* Precision prec1 is the precision on elementary computations; 406 it ensures a final precision prec1 - add for zeta(s) */ 407 /* eps = pow (2.0, - (double) precz - 14.0); */ 408 eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0); 409 m1 = 1.0 + MAX(1.0 / eps, 2.0 * sd) * (1.0 + eps); 410 c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1)); 411 /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */ 412 add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1)); 413 prec1 = precz + add; 414 prec1 = MAX (prec1, precs1) + 10; 415 416 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p); 417 MPFR_ZIV_INIT (loop, prec1); 418 for (;;) 419 { 420 mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN);/* s1 = 1-s */ 421 mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ 422 mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */ 423 if (MPFR_IS_INF (y)) /* Zeta(s) < 0 for -4k-2 < s < -4k, 424 Zeta(s) > 0 for -4k < s < -4k+2 */ 425 { 426 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ 427 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ 428 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1; 429 break; 430 } 431 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */ 432 mpfr_const_pi (p, MPFR_RNDD); 433 mpfr_mul (y, s, p, MPFR_RNDN); 434 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* s*Pi/2 */ 435 mpfr_sin (y, y, MPFR_RNDN); /* sin(Pi*s/2) */ 436 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 437 mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */ 438 mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */ 439 mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */ 440 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 441 mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN); 442 443 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz, 444 rnd_mode))) 445 break; 446 447 MPFR_ZIV_NEXT (loop, prec1); 448 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p); 449 } 450 MPFR_ZIV_FREE (loop); 451 if (overflow != 0) 452 { 453 inex = mpfr_overflow (z, rnd_mode, overflow); 454 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 455 } 456 else 457 inex = mpfr_set (z, z_pre, rnd_mode); 458 MPFR_GROUP_CLEAR (group); 459 } 460 461 MPFR_SAVE_EXPO_FREE (expo); 462 return mpfr_check_range (z, inex, rnd_mode); 463} 464