zeta.c revision 1.1.1.4
1/* mpfr_zeta -- compute the Riemann Zeta function 2 3Copyright 2003-2020 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 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 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include <float.h> /* for DBL_MAX */ 24 25#define MPFR_NEED_LONGLONG_H 26#include "mpfr-impl.h" 27 28/* 29 Parameters: 30 s - the input floating-point number 31 n, p - parameters from the algorithm 32 tc - an array of p floating-point numbers tc[1]..tc[p] 33 Output: 34 b is the result, i.e. 35 sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1) 36*/ 37static void 38mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc) 39{ 40 mpfr_t s1, d, u; 41 unsigned long n2; 42 int l, t; 43 MPFR_GROUP_DECL (group); 44 45 if (p == 0) 46 { 47 MPFR_SET_ZERO (b); 48 MPFR_SET_POS (b); 49 return; 50 } 51 52 n2 = n * n; 53 MPFR_GROUP_INIT_3 (group, MPFR_PREC (b), s1, d, u); 54 55 /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */ 56 t = 2 * p - 2; 57 mpfr_set (d, tc[p], MPFR_RNDN); 58 for (l = 1; l < p; l++) 59 { 60 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */ 61 mpfr_mul (d, d, s1, MPFR_RNDN); 62 t = t - 1; 63 mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */ 64 mpfr_mul (d, d, s1, MPFR_RNDN); 65 t = t - 1; 66 mpfr_div_ui (d, d, n2, MPFR_RNDN); 67 mpfr_add (d, d, tc[p-l], MPFR_RNDN); 68 /* since s is positive and the tc[i] have alternate signs, 69 the following is unlikely */ 70 if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0)) 71 mpfr_set (d, tc[p-l], MPFR_RNDN); 72 } 73 mpfr_mul (d, d, s, MPFR_RNDN); 74 mpfr_add (s1, s, __gmpfr_one, MPFR_RNDN); 75 mpfr_neg (s1, s1, MPFR_RNDN); 76 mpfr_ui_pow (u, n, s1, MPFR_RNDN); 77 mpfr_mul (b, d, u, MPFR_RNDN); 78 79 MPFR_GROUP_CLEAR (group); 80} 81 82/* Input: p - an integer 83 Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)! 84 tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ... 85 Assumes all the tc[i] have the same precision. 86 87 Uses the recurrence (4.60) from the book "Modern Computer Arithmetic" 88 by Brent and Zimmermann for C_k = bernoulli(2k)/(2k)!: 89 sum(C_k/(2k+1-2j)!/4^(k-j), j=0..k) = 1/(2k)!/4^k 90 If we put together the terms involving C_0 and C_1 we get: 91 sum(D_k/(2k+1-2j)!/4^(k-j), j=1..k) = 0 92 with D_1 = C_0/4/(2k+1)/(2k)+C_1-1/(2k)/4=(k-1)/(12k+6), 93 and D_k = C_k for k >= 2. 94 95 FIXME: we have C_k = (-1)^(k-1) 2/(2pi)^(2k) * zeta(2k), 96 see for example formula (4.65) from the above book, 97 thus since |zeta(2k)-1| < 2^(1-2k) for k >= 2, we have: 98 |C_k - E_k| < E_k * 2^(1-2k) for k >= 2 and E_k := (-1)^(k-1) 2/(2pi)^(2k). 99 Then if 2k-1 >= prec we can evaluate E_k instead, which only requires one 100 multiplication per term, instead of O(k) small divisions. 101*/ 102static void 103mpfr_zeta_c (int p, mpfr_t *tc) 104{ 105 if (p > 0) 106 { 107 mpfr_t d; 108 int k, l; 109 mpfr_prec_t prec = MPFR_PREC (tc[1]); 110 111 mpfr_init2 (d, prec); 112 mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN); 113 for (k = 2; k <= p; k++) 114 { 115 mpfr_set_ui (d, k-1, MPFR_RNDN); 116 mpfr_div_ui (d, d, 12*k+6, MPFR_RNDN); 117 for (l=2; l < k; l++) 118 { 119 mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), MPFR_RNDN); 120 mpfr_add (d, d, tc[l], MPFR_RNDN); 121 } 122 mpfr_div_ui (tc[k], d, 24, MPFR_RNDN); 123 MPFR_CHANGE_SIGN (tc[k]); 124 } 125 mpfr_clear (d); 126 } 127} 128 129/* Input: s - a floating-point number 130 n - an integer 131 Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */ 132static void 133mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n) 134{ 135 mpfr_t u, s1; 136 int i; 137 MPFR_GROUP_DECL (group); 138 139 MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1); 140 141 mpfr_neg (s1, s, MPFR_RNDN); 142 mpfr_ui_pow (u, n, s1, MPFR_RNDN); 143 mpfr_div_2ui (u, u, 1, MPFR_RNDN); 144 mpfr_set (sum, u, MPFR_RNDN); 145 for (i=n-1; i>1; i--) 146 { 147 mpfr_ui_pow (u, i, s1, MPFR_RNDN); 148 mpfr_add (sum, sum, u, MPFR_RNDN); 149 } 150 mpfr_add (sum, sum, __gmpfr_one, MPFR_RNDN); 151 152 MPFR_GROUP_CLEAR (group); 153} 154 155/* Input: s - a floating-point number >= 1/2. 156 rnd_mode - a rounding mode. 157 Assumes s is neither NaN nor Infinite. 158 Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode 159*/ 160static int 161mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) 162{ 163 mpfr_t b, c, z_pre, f, s1; 164 double beta, sd, dnep; 165 mpfr_t *tc1; 166 mpfr_prec_t precz, precs, d, dint; 167 int p, n, l, add; 168 int inex; 169 MPFR_GROUP_DECL (group); 170 MPFR_ZIV_DECL (loop); 171 172 MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0); 173 174 precz = MPFR_PREC (z); 175 precs = MPFR_PREC (s); 176 177 /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x) 178 so with 2^(EXP(x)-1) <= x < 2^EXP(x) 179 So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8 180 Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...) 181 = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity)) 182 <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity)) 183 And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035 184 So Zeta(x) <= 1 + 1/2^x*2 for x >= 8 185 The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */ 186 if (MPFR_GET_EXP (s) > 3) 187 { 188 mpfr_exp_t err; 189 err = MPFR_GET_EXP (s) - 1; 190 if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2)) 191 err = MPFR_EMAX_MAX; 192 else 193 err = ((mpfr_exp_t)1) << err; 194 err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */ 195 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (z, __gmpfr_one, err, 0, 1, 196 rnd_mode, {}); 197 } 198 199 d = precz + MPFR_INT_CEIL_LOG2(precz) + 10; 200 201 /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */ 202 dint = (mpfr_uexp_t) MPFR_GET_EXP (s); 203 mpfr_init2 (s1, MAX (precs, dint)); 204 inex = mpfr_sub (s1, s, __gmpfr_one, MPFR_RNDN); 205 MPFR_ASSERTD (inex == 0); 206 207 /* case s=1 should have already been handled */ 208 MPFR_ASSERTD (!MPFR_IS_ZERO (s1)); 209 210 MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f); 211 212 MPFR_ZIV_INIT (loop, d); 213 for (;;) 214 { 215 /* Principal loop: we compute, in z_pre, 216 an approximation of Zeta(s), that we send to can_round */ 217 if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2)) 218 /* Branch 1: when s-1 is very small, one 219 uses the approximation Zeta(s)=1/(s-1)+gamma, 220 where gamma is Euler's constant */ 221 { 222 dint = MAX (d + 3, precs); 223 /* branch 1, with internal precision dint */ 224 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f); 225 mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN); 226 mpfr_const_euler (f, MPFR_RNDN); 227 mpfr_add (z_pre, z_pre, f, MPFR_RNDN); 228 } 229 else /* Branch 2 */ 230 { 231 size_t size; 232 233 /* branch 2 */ 234 /* Computation of parameters n, p and working precision */ 235 dnep = (double) d * LOG2; 236 sd = mpfr_get_d (s, MPFR_RNDN); 237 /* beta = dnep + 0.61 + sd * log (6.2832 / sd); 238 but a larger value is OK */ 239#define LOG6dot2832 1.83787940484160805532 240 beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 * 241 __gmpfr_floor_log2 (sd)); 242 if (beta <= 0.0) 243 { 244 p = 0; 245 /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */ 246 n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd); 247 } 248 else 249 { 250 p = 1 + (int) beta / 2; 251 n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832); 252 } 253 /* add = 4 + floor(1.5 * log(d) / log (2)). 254 We should have add >= 10, which is always fulfilled since 255 d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */ 256 add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2; 257 MPFR_ASSERTD(add >= 10); 258 dint = d + add; 259 if (dint < precs) 260 dint = precs; 261 262 /* internal precision is dint */ 263 264 size = (p + 1) * sizeof(mpfr_t); 265 tc1 = (mpfr_t*) mpfr_allocate_func (size); 266 for (l=1; l<=p; l++) 267 mpfr_init2 (tc1[l], dint); 268 MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f); 269 270 /* precision of z is precz */ 271 272 /* Computation of the coefficients c_k */ 273 mpfr_zeta_c (p, tc1); 274 /* Computation of the 3 parts of the function Zeta. */ 275 mpfr_zeta_part_a (z_pre, s, n); 276 mpfr_zeta_part_b (b, s, n, p, tc1); 277 /* s1 = s-1 is already computed above */ 278 mpfr_div (c, __gmpfr_one, s1, MPFR_RNDN); 279 mpfr_ui_pow (f, n, s1, MPFR_RNDN); 280 mpfr_div (c, c, f, MPFR_RNDN); 281 mpfr_add (z_pre, z_pre, c, MPFR_RNDN); 282 mpfr_add (z_pre, z_pre, b, MPFR_RNDN); 283 for (l=1; l<=p; l++) 284 mpfr_clear (tc1[l]); 285 mpfr_free_func (tc1, size); 286 /* End branch 2 */ 287 } 288 289 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode))) 290 break; 291 MPFR_ZIV_NEXT (loop, d); 292 } 293 MPFR_ZIV_FREE (loop); 294 295 inex = mpfr_set (z, z_pre, rnd_mode); 296 297 MPFR_GROUP_CLEAR (group); 298 mpfr_clear (s1); 299 300 return inex; 301} 302 303/* return add = 1 + floor(log(c^3*(13+m1))/log(2)) 304 where c = (1+eps)*(1+eps*max(8,m1)), 305 m1 = 1 + max(1/eps,2*sd)*(1+eps), 306 eps = 2^(-precz-14) 307 sd = abs(s-1) 308 */ 309static long 310compute_add (mpfr_srcptr s, mpfr_prec_t precz) 311{ 312 mpfr_t t, u, m1; 313 long add; 314 315 mpfr_inits2 (64, t, u, m1, (mpfr_ptr) 0); 316 if (mpfr_cmp_ui (s, 1) >= 0) 317 mpfr_sub_ui (t, s, 1, MPFR_RNDU); 318 else 319 mpfr_ui_sub (t, 1, s, MPFR_RNDU); 320 /* now t = sd = abs(s-1), rounded up */ 321 mpfr_set_ui_2exp (u, 1, - precz - 14, MPFR_RNDU); 322 /* u = eps */ 323 /* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then 324 sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */ 325 if (mpfr_get_exp (t) >= precz + 14) 326 mpfr_mul_2ui (t, t, 1, MPFR_RNDU); 327 else 328 mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU); 329 /* now t = max(1/eps,2*sd) */ 330 mpfr_add_ui (u, u, 1, MPFR_RNDU); /* u = 1+eps, rounded up */ 331 mpfr_mul (t, t, u, MPFR_RNDU); /* t = max(1/eps,2*sd)*(1+eps) */ 332 mpfr_add_ui (m1, t, 1, MPFR_RNDU); 333 if (mpfr_get_exp (m1) <= 3) 334 mpfr_set_ui (t, 8, MPFR_RNDU); 335 else 336 mpfr_set (t, m1, MPFR_RNDU); 337 /* now t = max(8,m1) */ 338 mpfr_div_2ui (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */ 339 mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */ 340 mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */ 341 mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */ 342 mpfr_mul (u, u, t, MPFR_RNDU); /* c*(13+m1) */ 343 mpfr_sqr (t, t, MPFR_RNDU); /* c^2 */ 344 mpfr_mul (u, u, t, MPFR_RNDU); /* c^3*(13+m1) */ 345 add = mpfr_get_exp (u); 346 mpfr_clears (t, u, m1, (mpfr_ptr) 0); 347 return add; 348} 349 350/* return in z a lower bound (for rnd = RNDD) or upper bound (for rnd = RNDU) 351 of |zeta(s)|/2, using: 352 log(|zeta(s)|/2) = (s-1)*log(2*Pi) + lngamma(1-s) 353 + log(|sin(Pi*s/2)| * zeta(1-s)). 354 Assumes s < 1/2 and s1 = 1-s exactly, thus s1 > 1/2. 355 y and p are temporary variables. 356 At input, p is Pi rounded down. 357 The comments in the code are for rnd = RNDD. */ 358static void 359mpfr_reflection_overflow (mpfr_t z, mpfr_t s1, const mpfr_t s, mpfr_t y, 360 mpfr_t p, mpfr_rnd_t rnd) 361{ 362 mpz_t sint; 363 364 MPFR_ASSERTD (rnd == MPFR_RNDD || rnd == MPFR_RNDU); 365 366 /* Since log is increasing, we want lower bounds on |sin(Pi*s/2)| and 367 zeta(1-s). */ 368 mpz_init (sint); 369 mpfr_get_z (sint, s, MPFR_RNDD); /* sint = floor(s) */ 370 /* We first compute a lower bound of |sin(Pi*s/2)|, which is a periodic 371 function of period 2. Thus: 372 if 2k < s < 2k+1, then |sin(Pi*s/2)| is increasing; 373 if 2k-1 < s < 2k, then |sin(Pi*s/2)| is decreasing. 374 These cases are distinguished by testing bit 0 of floor(s) as if 375 represented in two's complement (or equivalently, as an unsigned 376 integer mod 2): 377 0: sint = 0 mod 2, thus 2k < s < 2k+1 and |sin(Pi*s/2)| is increasing; 378 1: sint = 1 mod 2, thus 2k-1 < s < 2k and |sin(Pi*s/2)| is decreasing. 379 Let's recall that the comments are for rnd = RNDD. */ 380 if (mpz_tstbit (sint, 0) == 0) /* |sin(Pi*s/2)| is increasing: round down 381 Pi*s to get a lower bound. */ 382 { 383 mpfr_mul (y, p, s, rnd); 384 if (rnd == MPFR_RNDD) 385 mpfr_nextabove (p); /* we will need p rounded above afterwards */ 386 } 387 else /* |sin(Pi*s/2)| is decreasing: round up Pi*s to get a lower bound. */ 388 { 389 if (rnd == MPFR_RNDD) 390 mpfr_nextabove (p); 391 mpfr_mul (y, p, s, MPFR_INVERT_RND(rnd)); 392 } 393 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* exact, rounding mode doesn't matter */ 394 /* The rounding direction of sin depends on its sign. We have: 395 if -4k-2 < s < -4k, then -2k-1 < s/2 < -2k, thus sin(Pi*s/2) < 0; 396 if -4k < s < -4k+2, then -2k < s/2 < -2k+1, thus sin(Pi*s/2) > 0. 397 These cases are distinguished by testing bit 1 of floor(s) as if 398 represented in two's complement (or equivalently, as an unsigned 399 integer mod 4): 400 0: sint = {0,1} mod 4, thus -2k < s/2 < -2k+1 and sin(Pi*s/2) > 0; 401 1: sint = {2,3} mod 4, thus -2k-1 < s/2 < -2k and sin(Pi*s/2) < 0. 402 Let's recall that the comments are for rnd = RNDD. */ 403 if (mpz_tstbit (sint, 1) == 0) /* -2k < s/2 < -2k+1; sin(Pi*s/2) > 0 */ 404 { 405 /* Round sin down to get a lower bound of |sin(Pi*s/2)|. */ 406 mpfr_sin (y, y, rnd); 407 } 408 else /* -2k-1 < s/2 < -2k; sin(Pi*s/2) < 0 */ 409 { 410 /* Round sin up to get a lower bound of |sin(Pi*s/2)|. */ 411 mpfr_sin (y, y, MPFR_INVERT_RND(rnd)); 412 mpfr_abs (y, y, MPFR_RNDN); /* exact, rounding mode doesn't matter */ 413 } 414 mpz_clear (sint); 415 /* now y <= |sin(Pi*s/2)| when rnd=RNDD, y >= |sin(Pi*s/2)| when rnd=RNDU */ 416 mpfr_zeta_pos (z, s1, rnd); /* zeta(1-s) */ 417 mpfr_mul (z, z, y, rnd); 418 /* now z <= |sin(Pi*s/2)|*zeta(1-s) */ 419 mpfr_log (z, z, rnd); 420 /* now z <= log(|sin(Pi*s/2)|*zeta(1-s)) */ 421 mpfr_lngamma (y, s1, rnd); 422 mpfr_add (z, z, y, rnd); 423 /* z <= lngamma(1-s) + log(|sin(Pi*s/2)|*zeta(1-s)) */ 424 /* since s-1 < 0, we want to round log(2*pi) upwards */ 425 mpfr_mul_2ui (y, p, 1, MPFR_INVERT_RND(rnd)); 426 mpfr_log (y, y, MPFR_INVERT_RND(rnd)); 427 mpfr_mul (y, y, s1, MPFR_INVERT_RND(rnd)); 428 mpfr_sub (z, z, y, rnd); 429 mpfr_exp (z, z, rnd); 430 if (rnd == MPFR_RNDD) 431 mpfr_nextbelow (p); /* restore original p */ 432} 433 434int 435mpfr_zeta (mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) 436{ 437 mpfr_t z_pre, s1, y, p; 438 long add; 439 mpfr_prec_t precz, prec1, precs, precs1; 440 int inex; 441 MPFR_GROUP_DECL (group); 442 MPFR_ZIV_DECL (loop); 443 MPFR_SAVE_EXPO_DECL (expo); 444 445 MPFR_LOG_FUNC ( 446 ("s[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode), 447 ("z[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex)); 448 449 /* Zero, Nan or Inf ? */ 450 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s))) 451 { 452 if (MPFR_IS_NAN (s)) 453 { 454 MPFR_SET_NAN (z); 455 MPFR_RET_NAN; 456 } 457 else if (MPFR_IS_INF (s)) 458 { 459 if (MPFR_IS_POS (s)) 460 return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */ 461 MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ 462 MPFR_RET_NAN; 463 } 464 else /* s iz zero */ 465 { 466 MPFR_ASSERTD (MPFR_IS_ZERO (s)); 467 return mpfr_set_si_2exp (z, -1, -1, rnd_mode); 468 } 469 } 470 471 /* s is neither Nan, nor Inf, nor Zero */ 472 473 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0, 474 and for |s| <= 2^(-4), we have |zeta(s) + 1/2| <= |s|. 475 EXP(s) + 1 < -PREC(z) is a sufficient condition to be able to round 476 correctly, for any PREC(z) >= 1 (see algorithms.tex for details). */ 477 if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z)) 478 { 479 int signs = MPFR_SIGN(s); 480 481 MPFR_SAVE_EXPO_MARK (expo); 482 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */ 483 if (rnd_mode == MPFR_RNDA) 484 rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */ 485 if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0) 486 { 487 mpfr_nextabove (z); /* z = -1/2 + epsilon */ 488 inex = 1; 489 } 490 else if (rnd_mode == MPFR_RNDD && signs > 0) 491 { 492 mpfr_nextbelow (z); /* z = -1/2 - epsilon */ 493 inex = -1; 494 } 495 else 496 { 497 if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */ 498 inex = 1; 499 else if (rnd_mode == MPFR_RNDD) 500 inex = -1; /* s < 0: z = -1/2 */ 501 else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */ 502 inex = (signs > 0) ? 1 : -1; 503 } 504 MPFR_SAVE_EXPO_FREE (expo); 505 return mpfr_check_range (z, inex, rnd_mode); 506 } 507 508 /* Check for case s= -2n */ 509 if (MPFR_IS_NEG (s)) 510 { 511 mpfr_t tmp; 512 tmp[0] = *s; 513 MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1; 514 if (mpfr_integer_p (tmp)) 515 { 516 MPFR_SET_ZERO (z); 517 MPFR_SET_POS (z); 518 MPFR_RET (0); 519 } 520 } 521 522 /* Check for case s=1 before changing the exponent range */ 523 if (mpfr_cmp (s, __gmpfr_one) == 0) 524 { 525 MPFR_SET_INF (z); 526 MPFR_SET_POS (z); 527 MPFR_SET_DIVBY0 (); 528 MPFR_RET (0); 529 } 530 531 MPFR_SAVE_EXPO_MARK (expo); 532 533 /* Compute Zeta */ 534 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ 535 inex = mpfr_zeta_pos (z, s, rnd_mode); 536 else /* use reflection formula 537 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */ 538 { 539 int overflow = 0; 540 541 precz = MPFR_PREC (z); 542 precs = MPFR_PREC (s); 543 544 /* Precision precs1 needed to represent 1 - s, and s + 2, 545 without any truncation */ 546 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); 547 /* Precision prec1 is the precision on elementary computations; 548 it ensures a final precision prec1 - add for zeta(s) */ 549 add = compute_add (s, precz); 550 prec1 = precz + add; 551 /* FIXME: To avoid that the working precision (prec1) depends on the 552 input precision, one would need to take into account the error made 553 when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1) 554 below, and also in the case y=Inf (i.e. when gamma(s1) overflows). 555 Make sure that underflows do not occur in intermediate computations. 556 Due to the limited precision, they are probably not possible 557 in practice; add some MPFR_ASSERTN's to be sure that problems 558 do not remain undetected? */ 559 prec1 = MAX (prec1, precs1) + 10; 560 561 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p); 562 MPFR_ZIV_INIT (loop, prec1); 563 for (;;) 564 { 565 mpfr_exp_t ey; 566 mpfr_t z_up; 567 568 mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */ 569 570 mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */ 571 mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */ 572 if (MPFR_IS_INF (y)) /* zeta(s) < 0 for -4k-2 < s < -4k, 573 zeta(s) > 0 for -4k < s < -4k+2 */ 574 { 575 /* FIXME: An overflow in gamma(s1) does not imply that 576 zeta(s) will overflow. A solution: 577 1. Compute 578 log(|zeta(s)|/2) = (s-1)*log(2*pi) + lngamma(1-s) 579 + log(abs(sin(Pi*s/2)) * zeta(1-s)) 580 (possibly sharing computations with the normal case) 581 with a rather good accuracy (see (2)). 582 Memorize the sign of sin(...) for the final sign. 583 2. Take the exponential, ~= |zeta(s)|/2. If there is an 584 overflow, then this means an overflow on the final result 585 (due to the multiplication by 2, which has not been done 586 yet). 587 3. Ziv test. 588 4. Correct the sign from the sign of sin(...). 589 5. Round then multiply by 2. Here, an overflow in either 590 operation means a real overflow. */ 591 mpfr_reflection_overflow (z_pre, s1, s, y, p, MPFR_RNDD); 592 /* z_pre is a lower bound of |zeta(s)|/2, thus if it overflows, 593 or has exponent emax, then |zeta(s)| overflows too. */ 594 if (MPFR_IS_INF (z_pre) || MPFR_GET_EXP(z_pre) == __gmpfr_emax) 595 { /* determine the sign of overflow */ 596 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ 597 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ 598 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1; 599 break; 600 } 601 else /* EXP(z_pre) < __gmpfr_emax */ 602 { 603 int ok = 0; 604 mpfr_t z_down; 605 mpfr_init2 (z_up, mpfr_get_prec (z_pre)); 606 mpfr_reflection_overflow (z_up, s1, s, y, p, MPFR_RNDU); 607 /* if the lower approximation z_pre does not overflow, but 608 z_up does, we need more precision */ 609 if (MPFR_IS_INF (z_up) || MPFR_GET_EXP(z_up) == __gmpfr_emax) 610 goto next_loop; 611 /* check if z_pre and z_up round to the same number */ 612 mpfr_init2 (z_down, precz); 613 mpfr_set (z_down, z_pre, rnd_mode); 614 /* Note: it might be that EXP(z_down) = emax here, in that 615 case we will have overflow below when we multiply by 2 */ 616 mpfr_prec_round (z_up, precz, rnd_mode); 617 ok = mpfr_cmp (z_down, z_up) == 0; 618 mpfr_clear (z_up); 619 mpfr_clear (z_down); 620 if (ok) 621 { 622 /* get correct sign and multiply by 2 */ 623 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ 624 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ 625 if (mpfr_cmp_si_2exp (s1, -1, -1) > 0) 626 mpfr_neg (z_pre, z_pre, rnd_mode); 627 mpfr_mul_2ui (z_pre, z_pre, 1, rnd_mode); 628 break; 629 } 630 else 631 goto next_loop; 632 } 633 } 634 mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ 635 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */ 636 637 /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */ 638 mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */ 639 mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */ 640 mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */ 641 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 642 mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN); 643 644 /* multiply z_pre by sin(Pi*s/2) */ 645 mpfr_mul (y, s, p, MPFR_RNDN); 646 mpfr_div_2ui (p, y, 1, MPFR_RNDN); /* p = s*Pi/2 */ 647 /* FIXME: sinpi will be available, we should replace the mpfr_sin 648 call below by mpfr_sinpi(s/2), where s/2 will be exact. 649 Can mpfr_sin underflow? Moreover, the code below should be 650 improved so that the "if" condition becomes unlikely, e.g. 651 by taking a slightly larger working precision. */ 652 mpfr_sin (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */ 653 ey = MPFR_GET_EXP (y); 654 if (ey < 0) /* take account of cancellation in sin(p) */ 655 { 656 mpfr_t t; 657 658 MPFR_ASSERTN (- ey < MPFR_PREC_MAX - prec1); 659 mpfr_init2 (t, prec1 - ey); 660 mpfr_const_pi (t, MPFR_RNDD); 661 mpfr_mul (t, s, t, MPFR_RNDN); 662 mpfr_div_2ui (t, t, 1, MPFR_RNDN); 663 mpfr_sin (y, t, MPFR_RNDN); 664 mpfr_clear (t); 665 } 666 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 667 668 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz, 669 rnd_mode))) 670 break; 671 672 next_loop: 673 MPFR_ZIV_NEXT (loop, prec1); 674 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p); 675 } 676 MPFR_ZIV_FREE (loop); 677 if (overflow != 0) 678 { 679 inex = mpfr_overflow (z, rnd_mode, overflow); 680 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 681 } 682 else 683 inex = mpfr_set (z, z_pre, rnd_mode); 684 MPFR_GROUP_CLEAR (group); 685 } 686 687 MPFR_SAVE_EXPO_FREE (expo); 688 return mpfr_check_range (z, inex, rnd_mode); 689} 690