1 /* mpfr_zeta -- compute the Riemann Zeta function 2 3 Copyright 2003-2023 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 https://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 #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 */ 37 static void 38 mpfr_zeta_part_b (mpfr_ptr 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 */ 102 static void 103 mpfr_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) */ 132 static void 133 mpfr_zeta_part_a (mpfr_ptr 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 */ 160 static int 161 mpfr_zeta_pos (mpfr_ptr 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 /* TODO: Check the error analysis. The following (undocumented?) one 304 does not take into account the replacement of sin(Pi*s/2) by sinpi(s/2) 305 in commit fd5d811d81f6d1839d4099cc1bb2cde705981648, which could have 306 reduced the error bound since the multiplication by Pi is now exact. */ 307 /* return add = 1 + floor(log(c^3*(13+m1))/log(2)) 308 where c = (1+eps)*(1+eps*max(8,m1)), 309 m1 = 1 + max(1/eps,2*sd)*(1+eps), 310 eps = 2^(-precz-14) 311 sd = abs(s-1) 312 */ 313 static long 314 compute_add (mpfr_srcptr s, mpfr_prec_t precz) 315 { 316 mpfr_t t, u, m1; 317 long add; 318 319 mpfr_inits2 (64, t, u, m1, (mpfr_ptr) 0); 320 if (mpfr_cmp_ui (s, 1) >= 0) 321 mpfr_sub_ui (t, s, 1, MPFR_RNDU); 322 else 323 mpfr_ui_sub (t, 1, s, MPFR_RNDU); 324 /* now t = sd = abs(s-1), rounded up */ 325 mpfr_set_ui_2exp (u, 1, - precz - 14, MPFR_RNDU); 326 /* u = eps */ 327 /* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then 328 sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */ 329 if (mpfr_get_exp (t) >= precz + 14) 330 mpfr_mul_2ui (t, t, 1, MPFR_RNDU); 331 else 332 mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU); 333 /* now t = max(1/eps,2*sd) */ 334 mpfr_add_ui (u, u, 1, MPFR_RNDU); /* u = 1+eps, rounded up */ 335 mpfr_mul (t, t, u, MPFR_RNDU); /* t = max(1/eps,2*sd)*(1+eps) */ 336 mpfr_add_ui (m1, t, 1, MPFR_RNDU); 337 if (mpfr_get_exp (m1) <= 3) 338 mpfr_set_ui (t, 8, MPFR_RNDU); 339 else 340 mpfr_set (t, m1, MPFR_RNDU); 341 /* now t = max(8,m1) */ 342 mpfr_div_2ui (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */ 343 mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */ 344 mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */ 345 mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */ 346 mpfr_mul (u, u, t, MPFR_RNDU); /* c*(13+m1) */ 347 mpfr_sqr (t, t, MPFR_RNDU); /* c^2 */ 348 mpfr_mul (u, u, t, MPFR_RNDU); /* c^3*(13+m1) */ 349 add = mpfr_get_exp (u); 350 mpfr_clears (t, u, m1, (mpfr_ptr) 0); 351 return add; 352 } 353 354 /* return in z a lower bound (for rnd = RNDD) or upper bound (for rnd = RNDU) 355 of |zeta(s)|/2, using: 356 log(|zeta(s)|/2) = (s-1)*log(2*Pi) + lngamma(1-s) 357 + log(|sin(Pi*s/2)| * zeta(1-s)). 358 Assumes s < 1/2 and s1 = 1-s exactly, thus s1 > 1/2. 359 y and p are temporary variables. 360 At input, p is Pi rounded down. 361 The comments in the code are for rnd = RNDD. */ 362 static void 363 mpfr_reflection_overflow (mpfr_ptr z, mpfr_ptr s1, mpfr_srcptr s, mpfr_ptr y, 364 mpfr_ptr p, mpfr_rnd_t rnd) 365 { 366 mpz_t sint; 367 368 MPFR_ASSERTD (rnd == MPFR_RNDD || rnd == MPFR_RNDU); 369 370 /* Since log is increasing, we want lower bounds on |sin(Pi*s/2)| and 371 zeta(1-s). */ 372 mpz_init (sint); 373 mpfr_get_z (sint, s, MPFR_RNDD); /* sint = floor(s) */ 374 /* We first compute a lower bound of |sin(Pi*s/2)|, which is a periodic 375 function of period 2. Thus: 376 if 2k < s < 2k+1, then |sin(Pi*s/2)| is increasing; 377 if 2k-1 < s < 2k, then |sin(Pi*s/2)| is decreasing. 378 These cases are distinguished by testing bit 0 of floor(s) as if 379 represented in two's complement (or equivalently, as an unsigned 380 integer mod 2): 381 0: sint = 0 mod 2, thus 2k < s < 2k+1 and |sin(Pi*s/2)| is increasing; 382 1: sint = 1 mod 2, thus 2k-1 < s < 2k and |sin(Pi*s/2)| is decreasing. 383 Let's recall that the comments are for rnd = RNDD. */ 384 if (mpz_tstbit (sint, 0) == 0) /* |sin(Pi*s/2)| is increasing: round down 385 Pi*s to get a lower bound. */ 386 { 387 mpfr_mul (y, p, s, rnd); 388 if (rnd == MPFR_RNDD) 389 mpfr_nextabove (p); /* we will need p rounded above afterwards */ 390 } 391 else /* |sin(Pi*s/2)| is decreasing: round up Pi*s to get a lower bound. */ 392 { 393 if (rnd == MPFR_RNDD) 394 mpfr_nextabove (p); 395 mpfr_mul (y, p, s, MPFR_INVERT_RND(rnd)); 396 } 397 mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* exact, rounding mode doesn't matter */ 398 /* The rounding direction of sin depends on its sign. We have: 399 if -4k-2 < s < -4k, then -2k-1 < s/2 < -2k, thus sin(Pi*s/2) < 0; 400 if -4k < s < -4k+2, then -2k < s/2 < -2k+1, thus sin(Pi*s/2) > 0. 401 These cases are distinguished by testing bit 1 of floor(s) as if 402 represented in two's complement (or equivalently, as an unsigned 403 integer mod 4): 404 0: sint = {0,1} mod 4, thus -2k < s/2 < -2k+1 and sin(Pi*s/2) > 0; 405 1: sint = {2,3} mod 4, thus -2k-1 < s/2 < -2k and sin(Pi*s/2) < 0. 406 Let's recall that the comments are for rnd = RNDD. */ 407 if (mpz_tstbit (sint, 1) == 0) /* -2k < s/2 < -2k+1; sin(Pi*s/2) > 0 */ 408 { 409 /* Round sin down to get a lower bound of |sin(Pi*s/2)|. */ 410 mpfr_sin (y, y, rnd); 411 } 412 else /* -2k-1 < s/2 < -2k; sin(Pi*s/2) < 0 */ 413 { 414 /* Round sin up to get a lower bound of |sin(Pi*s/2)|. */ 415 mpfr_sin (y, y, MPFR_INVERT_RND(rnd)); 416 mpfr_abs (y, y, MPFR_RNDN); /* exact, rounding mode doesn't matter */ 417 } 418 mpz_clear (sint); 419 /* now y <= |sin(Pi*s/2)| when rnd=RNDD, y >= |sin(Pi*s/2)| when rnd=RNDU */ 420 mpfr_zeta_pos (z, s1, rnd); /* zeta(1-s) */ 421 mpfr_mul (z, z, y, rnd); 422 /* now z <= |sin(Pi*s/2)|*zeta(1-s) */ 423 mpfr_log (z, z, rnd); 424 /* now z <= log(|sin(Pi*s/2)|*zeta(1-s)) */ 425 mpfr_lngamma (y, s1, rnd); 426 mpfr_add (z, z, y, rnd); 427 /* z <= lngamma(1-s) + log(|sin(Pi*s/2)|*zeta(1-s)) */ 428 /* since s-1 < 0, we want to round log(2*pi) upwards */ 429 mpfr_mul_2ui (y, p, 1, MPFR_INVERT_RND(rnd)); 430 mpfr_log (y, y, MPFR_INVERT_RND(rnd)); 431 mpfr_mul (y, y, s1, MPFR_INVERT_RND(rnd)); 432 mpfr_sub (z, z, y, rnd); 433 mpfr_exp (z, z, rnd); 434 if (rnd == MPFR_RNDD) 435 mpfr_nextbelow (p); /* restore original p */ 436 } 437 438 int 439 mpfr_zeta (mpfr_ptr z, mpfr_srcptr s, mpfr_rnd_t rnd_mode) 440 { 441 mpfr_t z_pre, s1, y, p; 442 long add; 443 mpfr_prec_t precz, prec1, precs, precs1; 444 int inex; 445 MPFR_GROUP_DECL (group); 446 MPFR_ZIV_DECL (loop); 447 MPFR_SAVE_EXPO_DECL (expo); 448 449 MPFR_LOG_FUNC ( 450 ("s[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode), 451 ("z[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex)); 452 453 /* Zero, Nan or Inf ? */ 454 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s))) 455 { 456 if (MPFR_IS_NAN (s)) 457 { 458 MPFR_SET_NAN (z); 459 MPFR_RET_NAN; 460 } 461 else if (MPFR_IS_INF (s)) 462 { 463 if (MPFR_IS_POS (s)) 464 return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */ 465 MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */ 466 MPFR_RET_NAN; 467 } 468 else /* s iz zero */ 469 { 470 MPFR_ASSERTD (MPFR_IS_ZERO (s)); 471 return mpfr_set_si_2exp (z, -1, -1, rnd_mode); 472 } 473 } 474 475 /* s is neither Nan, nor Inf, nor Zero */ 476 477 /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0, 478 and for |s| <= 2^(-4), we have |zeta(s) + 1/2| <= |s|. 479 EXP(s) + 1 < -PREC(z) is a sufficient condition to be able to round 480 correctly, for any PREC(z) >= 1 (see algorithms.tex for details). */ 481 if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z)) 482 { 483 int signs = MPFR_SIGN(s); 484 485 MPFR_SAVE_EXPO_MARK (expo); 486 mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */ 487 if (rnd_mode == MPFR_RNDA) 488 rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */ 489 if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0) 490 { 491 mpfr_nextabove (z); /* z = -1/2 + epsilon */ 492 inex = 1; 493 } 494 else if (rnd_mode == MPFR_RNDD && signs > 0) 495 { 496 mpfr_nextbelow (z); /* z = -1/2 - epsilon */ 497 inex = -1; 498 } 499 else 500 { 501 if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */ 502 inex = 1; 503 else if (rnd_mode == MPFR_RNDD) 504 inex = -1; /* s < 0: z = -1/2 */ 505 else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */ 506 inex = (signs > 0) ? 1 : -1; 507 } 508 MPFR_SAVE_EXPO_FREE (expo); 509 return mpfr_check_range (z, inex, rnd_mode); 510 } 511 512 /* Check for case s= -2n */ 513 if (MPFR_IS_NEG (s)) 514 { 515 mpfr_t tmp; 516 tmp[0] = *s; 517 MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1; 518 if (mpfr_integer_p (tmp)) 519 { 520 MPFR_SET_ZERO (z); 521 MPFR_SET_POS (z); 522 MPFR_RET (0); 523 } 524 } 525 526 /* Check for case s=1 before changing the exponent range */ 527 if (mpfr_equal_p (s, __gmpfr_one)) 528 { 529 MPFR_SET_INF (z); 530 MPFR_SET_POS (z); 531 MPFR_SET_DIVBY0 (); 532 MPFR_RET (0); 533 } 534 535 MPFR_SAVE_EXPO_MARK (expo); 536 537 /* Compute Zeta */ 538 if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */ 539 inex = mpfr_zeta_pos (z, s, rnd_mode); 540 else /* use reflection formula 541 zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */ 542 { 543 int overflow = 0; 544 545 precz = MPFR_PREC (z); 546 precs = MPFR_PREC (s); 547 548 /* Precision precs1 needed to represent 1 - s, and s + 2, 549 without any truncation */ 550 precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s)); 551 /* Precision prec1 is the precision on elementary computations; 552 it ensures a final precision prec1 - add for zeta(s) */ 553 add = compute_add (s, precz); 554 prec1 = precz + add; 555 /* FIXME: To avoid that the working precision (prec1) depends on the 556 input precision, one would need to take into account the error made 557 when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1) 558 below, and also in the case y=Inf (i.e. when gamma(s1) overflows). 559 Make sure that underflows do not occur in intermediate computations. 560 Due to the limited precision, they are probably not possible 561 in practice; add some MPFR_ASSERTN's to be sure that problems 562 do not remain undetected? */ 563 prec1 = MAX (prec1, precs1) + 10; 564 565 MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p); 566 MPFR_ZIV_INIT (loop, prec1); 567 for (;;) 568 { 569 mpfr_t z_up; 570 571 mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */ 572 573 mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */ 574 mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */ 575 if (MPFR_IS_INF (y)) /* zeta(s) < 0 for -4k-2 < s < -4k, 576 zeta(s) > 0 for -4k < s < -4k+2 */ 577 { 578 /* FIXME: An overflow in gamma(s1) does not imply that 579 zeta(s) will overflow. A solution: 580 1. Compute 581 log(|zeta(s)|/2) = (s-1)*log(2*pi) + lngamma(1-s) 582 + log(abs(sin(Pi*s/2)) * zeta(1-s)) 583 (possibly sharing computations with the normal case) 584 with a rather good accuracy (see (2)). 585 Memorize the sign of sin(...) for the final sign. 586 2. Take the exponential, ~= |zeta(s)|/2. If there is an 587 overflow, then this means an overflow on the final result 588 (due to the multiplication by 2, which has not been done 589 yet). 590 3. Ziv test. 591 4. Correct the sign from the sign of sin(...). 592 5. Round then multiply by 2. Here, an overflow in either 593 operation means a real overflow. */ 594 mpfr_reflection_overflow (z_pre, s1, s, y, p, MPFR_RNDD); 595 /* z_pre is a lower bound of |zeta(s)|/2, thus if it overflows, 596 or has exponent emax, then |zeta(s)| overflows too. */ 597 if (MPFR_IS_INF (z_pre) || MPFR_GET_EXP(z_pre) == __gmpfr_emax) 598 { /* determine the sign of overflow */ 599 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ 600 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ 601 overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1; 602 break; 603 } 604 else /* EXP(z_pre) < __gmpfr_emax */ 605 { 606 int ok = 0; 607 mpfr_t z_down; 608 mpfr_init2 (z_up, mpfr_get_prec (z_pre)); 609 mpfr_reflection_overflow (z_up, s1, s, y, p, MPFR_RNDU); 610 /* if the lower approximation z_pre does not overflow, but 611 z_up does, we need more precision */ 612 if (MPFR_IS_INF (z_up) || MPFR_GET_EXP(z_up) == __gmpfr_emax) 613 goto next_loop; 614 /* check if z_pre and z_up round to the same number */ 615 mpfr_init2 (z_down, precz); 616 mpfr_set (z_down, z_pre, rnd_mode); 617 /* Note: it might be that EXP(z_down) = emax here, in that 618 case we will have overflow below when we multiply by 2 */ 619 mpfr_prec_round (z_up, precz, rnd_mode); 620 ok = mpfr_equal_p (z_down, z_up); 621 mpfr_clear (z_up); 622 mpfr_clear (z_down); 623 if (ok) 624 { 625 /* get correct sign and multiply by 2 */ 626 mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */ 627 mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */ 628 if (mpfr_cmp_si_2exp (s1, -1, -1) > 0) 629 mpfr_neg (z_pre, z_pre, rnd_mode); 630 mpfr_mul_2ui (z_pre, z_pre, 1, rnd_mode); 631 break; 632 } 633 else 634 goto next_loop; 635 } 636 } 637 mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */ 638 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */ 639 640 /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */ 641 mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */ 642 mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */ 643 mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */ 644 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 645 mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN); 646 647 /* multiply z_pre by sin(Pi*s/2) */ 648 mpfr_div_2ui (p, s, 1, MPFR_RNDN); /* p = s/2 */ 649 /* Can mpfr_sinpi underflow? While with mpfr_sin, we could not 650 answer in any precision without a theoretical study (though 651 an underflow would have been very unlikely as we have a 652 huge exponent range), with mpfr_sinpi, an underflow could 653 occur only in a huge, unsupported precision. Indeed, if 654 mpfr_sinpi underflows, this means that 0 < |sinpi(s/2)| < m, 655 where m is the minimum representable positive number, and in 656 this case, r being the reduced argument such that |r| <= 1/2, 657 one has |sinpi(r)| > |2r|, so that |2r| < m; this can be 658 possible only if |s/2| > 1/2 (otherwise |s| = |2r| < m and 659 s would not be representable as an MPFR number) and s has 660 non-zero bits of exponent less than the minimum exponent 661 (s/2 - r being an integer), i.e. the precision is at least 662 MPFR_EMAX_MAX + 2. With such a huge precision, there would 663 probably be failures before reaching this point. */ 664 mpfr_sinpi (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */ 665 mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); 666 667 if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz, 668 rnd_mode))) 669 break; 670 671 next_loop: 672 MPFR_ZIV_NEXT (loop, prec1); 673 MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p); 674 } 675 MPFR_ZIV_FREE (loop); 676 if (overflow != 0) 677 { 678 inex = mpfr_overflow (z, rnd_mode, overflow); 679 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); 680 } 681 else 682 inex = mpfr_set (z, z_pre, rnd_mode); 683 MPFR_GROUP_CLEAR (group); 684 } 685 686 MPFR_SAVE_EXPO_FREE (expo); 687 return mpfr_check_range (z, inex, rnd_mode); 688 } 689
Indexes created Sat Sep 05 00:25:07 UTC 2026