1 /* mpfr_li2 -- Dilogarithm. 2 3 Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC 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 /* Compute the alternating series 27 s = S(z) = \sum_{k=0}^infty B_{2k} (z))^{2k+1} / (2k+1)! 28 with 0 < z <= log(2) to the precision of s rounded in the direction 29 rnd_mode. 30 Return the maximum index of the truncature which is useful 31 for determinating the relative error. 32 */ 33 static int 34 li2_series (mpfr_t sum, mpfr_srcptr z, mpfr_rnd_t rnd_mode) 35 { 36 int i, Bm, Bmax; 37 mpfr_t s, u, v, w; 38 mpfr_prec_t sump, p; 39 mpfr_exp_t se, err; 40 mpz_t *B; 41 MPFR_ZIV_DECL (loop); 42 43 /* The series converges for |z| < 2 pi, but in mpfr_li2 the argument is 44 reduced so that 0 < z <= log(2). Here is additionnal check that z is 45 (nearly) correct */ 46 MPFR_ASSERTD (MPFR_IS_STRICTPOS (z)); 47 MPFR_ASSERTD (mpfr_cmp_d (z, 0.6953125) <= 0); 48 49 sump = MPFR_PREC (sum); /* target precision */ 50 p = sump + MPFR_INT_CEIL_LOG2 (sump) + 4; /* the working precision */ 51 mpfr_init2 (s, p); 52 mpfr_init2 (u, p); 53 mpfr_init2 (v, p); 54 mpfr_init2 (w, p); 55 56 B = mpfr_bernoulli_internal ((mpz_t *) 0, 0); 57 Bm = Bmax = 1; 58 59 MPFR_ZIV_INIT (loop, p); 60 for (;;) 61 { 62 mpfr_sqr (u, z, MPFR_RNDU); 63 mpfr_set (v, z, MPFR_RNDU); 64 mpfr_set (s, z, MPFR_RNDU); 65 se = MPFR_GET_EXP (s); 66 err = 0; 67 68 for (i = 1;; i++) 69 { 70 if (i >= Bmax) 71 B = mpfr_bernoulli_internal (B, Bmax++); /* B_2i*(2i+1)!, exact */ 72 73 mpfr_mul (v, u, v, MPFR_RNDU); 74 mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); 75 mpfr_div_ui (v, v, 2 * i, MPFR_RNDU); 76 mpfr_div_ui (v, v, 2 * i + 1, MPFR_RNDU); 77 mpfr_div_ui (v, v, 2 * i + 1, MPFR_RNDU); 78 /* here, v_2i = v_{2i-2} / (2i * (2i+1))^2 */ 79 80 mpfr_mul_z (w, v, B[i], MPFR_RNDN); 81 /* here, w_2i = v_2i * B_2i * (2i+1)! with 82 error(w_2i) < 2^(5 * i + 8) ulp(w_2i) (see algorithms.tex) */ 83 84 mpfr_add (s, s, w, MPFR_RNDN); 85 86 err = MAX (err + se, 5 * i + 8 + MPFR_GET_EXP (w)) 87 - MPFR_GET_EXP (s); 88 err = 2 + MAX (-1, err); 89 se = MPFR_GET_EXP (s); 90 if (MPFR_GET_EXP (w) <= se - (mpfr_exp_t) p) 91 break; 92 } 93 94 /* the previous value of err is the rounding error, 95 the truncation error is less than EXP(z) - 6 * i - 5 96 (see algorithms.tex) */ 97 err = MAX (err, MPFR_GET_EXP (z) - 6 * i - 5) + 1; 98 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) p - err, sump, rnd_mode)) 99 break; 100 101 MPFR_ZIV_NEXT (loop, p); 102 mpfr_set_prec (s, p); 103 mpfr_set_prec (u, p); 104 mpfr_set_prec (v, p); 105 mpfr_set_prec (w, p); 106 } 107 MPFR_ZIV_FREE (loop); 108 mpfr_set (sum, s, rnd_mode); 109 110 Bm = Bmax; 111 while (Bm--) 112 mpz_clear (B[Bm]); 113 (*__gmp_free_func) (B, Bmax * sizeof (mpz_t)); 114 mpfr_clears (s, u, v, w, (mpfr_ptr) 0); 115 116 /* Let K be the returned value. 117 1. As we compute an alternating series, the truncation error has the same 118 sign as the next term w_{K+2} which is positive iff K%4 == 0. 119 2. Assume that error(z) <= (1+t) z', where z' is the actual value, then 120 error(s) <= 2 * (K+1) * t (see algorithms.tex). 121 */ 122 return 2 * i; 123 } 124 125 /* try asymptotic expansion when x is large and positive: 126 Li2(x) = -log(x)^2/2 + Pi^2/3 - 1/x + O(1/x^2). 127 More precisely for x >= 2 we have for g(x) = -log(x)^2/2 + Pi^2/3: 128 -2 <= x * (Li2(x) - g(x)) <= -1 129 thus |Li2(x) - g(x)| <= 2/x. 130 Assumes x >= 38, which ensures log(x)^2/2 >= 2*Pi^2/3, and g(x) <= -3.3. 131 Return 0 if asymptotic expansion failed (unable to round), otherwise 132 returns correct ternary value. 133 */ 134 static int 135 mpfr_li2_asympt_pos (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 136 { 137 mpfr_t g, h; 138 mpfr_prec_t w = MPFR_PREC (y) + 20; 139 int inex = 0; 140 141 MPFR_ASSERTN (mpfr_cmp_ui (x, 38) >= 0); 142 143 mpfr_init2 (g, w); 144 mpfr_init2 (h, w); 145 mpfr_log (g, x, MPFR_RNDN); /* rel. error <= |(1 + theta) - 1| */ 146 mpfr_sqr (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */ 147 mpfr_div_2ui (g, g, 1, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 148 mpfr_const_pi (h, MPFR_RNDN); /* error <= 2^(1-w) */ 149 mpfr_sqr (h, h, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 150 mpfr_div_ui (h, h, 3, MPFR_RNDN); /* rel. error <= |(1 + theta)^4 - 1| 151 <= 5 * 2^(-w) */ 152 /* since x is chosen such that log(x)^2/2 >= 2 * (Pi^2/3), we should have 153 g >= 2*h, thus |g-h| >= |h|, and the relative error on g is at most 154 multiplied by 2 in the difference, and that by h is unchanged. */ 155 MPFR_ASSERTN (MPFR_EXP (g) > MPFR_EXP (h)); 156 mpfr_sub (g, h, g, MPFR_RNDN); /* err <= ulp(g)/2 + g*2^(3-w) + g*5*2^(-w) 157 <= ulp(g) * (1/2 + 8 + 5) < 14 ulp(g). 158 159 If in addition 2/x <= 2 ulp(g), i.e., 160 1/x <= ulp(g), then the total error is 161 bounded by 16 ulp(g). */ 162 if ((MPFR_EXP (x) >= (mpfr_exp_t) w - MPFR_EXP (g)) && 163 MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) 164 inex = mpfr_set (y, g, rnd_mode); 165 166 mpfr_clear (g); 167 mpfr_clear (h); 168 169 return inex; 170 } 171 172 /* try asymptotic expansion when x is large and negative: 173 Li2(x) = -log(-x)^2/2 - Pi^2/6 - 1/x + O(1/x^2). 174 More precisely for x <= -2 we have for g(x) = -log(-x)^2/2 - Pi^2/6: 175 |Li2(x) - g(x)| <= 1/|x|. 176 Assumes x <= -7, which ensures |log(-x)^2/2| >= Pi^2/6, and g(x) <= -3.5. 177 Return 0 if asymptotic expansion failed (unable to round), otherwise 178 returns correct ternary value. 179 */ 180 static int 181 mpfr_li2_asympt_neg (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 182 { 183 mpfr_t g, h; 184 mpfr_prec_t w = MPFR_PREC (y) + 20; 185 int inex = 0; 186 187 MPFR_ASSERTN (mpfr_cmp_si (x, -7) <= 0); 188 189 mpfr_init2 (g, w); 190 mpfr_init2 (h, w); 191 mpfr_neg (g, x, MPFR_RNDN); 192 mpfr_log (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta) - 1| */ 193 mpfr_sqr (g, g, MPFR_RNDN); /* rel. error <= |(1 + theta)^3 - 1| <= 2^(2-w) */ 194 mpfr_div_2ui (g, g, 1, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 195 mpfr_const_pi (h, MPFR_RNDN); /* error <= 2^(1-w) */ 196 mpfr_sqr (h, h, MPFR_RNDN); /* rel. error <= 2^(2-w) */ 197 mpfr_div_ui (h, h, 6, MPFR_RNDN); /* rel. error <= |(1 + theta)^4 - 1| 198 <= 5 * 2^(-w) */ 199 MPFR_ASSERTN (MPFR_EXP (g) >= MPFR_EXP (h)); 200 mpfr_add (g, g, h, MPFR_RNDN); /* err <= ulp(g)/2 + g*2^(2-w) + g*5*2^(-w) 201 <= ulp(g) * (1/2 + 4 + 5) < 10 ulp(g). 202 203 If in addition |1/x| <= 4 ulp(g), then the 204 total error is bounded by 16 ulp(g). */ 205 if ((MPFR_EXP (x) >= (mpfr_exp_t) (w - 2) - MPFR_EXP (g)) && 206 MPFR_CAN_ROUND (g, w - 4, MPFR_PREC (y), rnd_mode)) 207 inex = mpfr_neg (y, g, rnd_mode); 208 209 mpfr_clear (g); 210 mpfr_clear (h); 211 212 return inex; 213 } 214 215 /* Compute the real part of the dilogarithm defined by 216 Li2(x) = -\Int_{t=0}^x log(1-t)/t dt */ 217 int 218 mpfr_li2 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 219 { 220 int inexact; 221 mpfr_exp_t err; 222 mpfr_prec_t yp, m; 223 MPFR_ZIV_DECL (loop); 224 MPFR_SAVE_EXPO_DECL (expo); 225 226 MPFR_LOG_FUNC 227 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 228 ("y[%Pu]=%.*Rg inexact=%d", 229 mpfr_get_prec (y), mpfr_log_prec, y, inexact)); 230 231 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 232 { 233 if (MPFR_IS_NAN (x)) 234 { 235 MPFR_SET_NAN (y); 236 MPFR_RET_NAN; 237 } 238 else if (MPFR_IS_INF (x)) 239 { 240 MPFR_SET_NEG (y); 241 MPFR_SET_INF (y); 242 MPFR_RET (0); 243 } 244 else /* x is zero */ 245 { 246 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 247 MPFR_SET_SAME_SIGN (y, x); 248 MPFR_SET_ZERO (y); 249 MPFR_RET (0); 250 } 251 } 252 253 /* Li2(x) = x + x^2/4 + x^3/9 + ..., more precisely for 0 < x <= 1/2 254 we have |Li2(x) - x| < x^2/2 <= 2^(2EXP(x)-1) and for -1/2 <= x < 0 255 we have |Li2(x) - x| < x^2/4 <= 2^(2EXP(x)-2) */ 256 if (MPFR_IS_POS (x)) 257 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x), 1, 1, rnd_mode, 258 {}); 259 else 260 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -MPFR_GET_EXP (x), 2, 0, rnd_mode, 261 {}); 262 263 MPFR_SAVE_EXPO_MARK (expo); 264 yp = MPFR_PREC (y); 265 m = yp + MPFR_INT_CEIL_LOG2 (yp) + 13; 266 267 if (MPFR_LIKELY ((mpfr_cmp_ui (x, 0) > 0) && (mpfr_cmp_d (x, 0.5) <= 0))) 268 /* 0 < x <= 1/2: Li2(x) = S(-log(1-x))-log^2(1-x)/4 */ 269 { 270 mpfr_t s, u; 271 mpfr_exp_t expo_l; 272 int k; 273 274 mpfr_init2 (u, m); 275 mpfr_init2 (s, m); 276 277 MPFR_ZIV_INIT (loop, m); 278 for (;;) 279 { 280 mpfr_ui_sub (u, 1, x, MPFR_RNDN); 281 mpfr_log (u, u, MPFR_RNDU); 282 if (MPFR_IS_ZERO(u)) 283 goto next_m; 284 mpfr_neg (u, u, MPFR_RNDN); /* u = -log(1-x) */ 285 expo_l = MPFR_GET_EXP (u); 286 k = li2_series (s, u, MPFR_RNDU); 287 err = 1 + MPFR_INT_CEIL_LOG2 (k + 1); 288 289 mpfr_sqr (u, u, MPFR_RNDU); 290 mpfr_div_2ui (u, u, 2, MPFR_RNDU); /* u = log^2(1-x) / 4 */ 291 mpfr_sub (s, s, u, MPFR_RNDN); 292 293 /* error(s) <= (0.5 + 2^(d-EXP(s)) 294 + 2^(3 + MAX(1, - expo_l) - EXP(s))) ulp(s) */ 295 err = MAX (err, MAX (1, - expo_l) - 1) - MPFR_GET_EXP (s); 296 err = 2 + MAX (-1, err); 297 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 298 break; 299 300 next_m: 301 MPFR_ZIV_NEXT (loop, m); 302 mpfr_set_prec (u, m); 303 mpfr_set_prec (s, m); 304 } 305 MPFR_ZIV_FREE (loop); 306 inexact = mpfr_set (y, s, rnd_mode); 307 308 mpfr_clear (u); 309 mpfr_clear (s); 310 MPFR_SAVE_EXPO_FREE (expo); 311 return mpfr_check_range (y, inexact, rnd_mode); 312 } 313 else if (!mpfr_cmp_ui (x, 1)) 314 /* Li2(1)= pi^2 / 6 */ 315 { 316 mpfr_t u; 317 mpfr_init2 (u, m); 318 319 MPFR_ZIV_INIT (loop, m); 320 for (;;) 321 { 322 mpfr_const_pi (u, MPFR_RNDU); 323 mpfr_sqr (u, u, MPFR_RNDN); 324 mpfr_div_ui (u, u, 6, MPFR_RNDN); 325 326 err = m - 4; /* error(u) <= 19/2 ulp(u) */ 327 if (MPFR_CAN_ROUND (u, err, yp, rnd_mode)) 328 break; 329 330 MPFR_ZIV_NEXT (loop, m); 331 mpfr_set_prec (u, m); 332 } 333 MPFR_ZIV_FREE (loop); 334 inexact = mpfr_set (y, u, rnd_mode); 335 336 mpfr_clear (u); 337 MPFR_SAVE_EXPO_FREE (expo); 338 return mpfr_check_range (y, inexact, rnd_mode); 339 } 340 else if (mpfr_cmp_ui (x, 2) >= 0) 341 /* x >= 2: Li2(x) = -S(-log(1-1/x))-log^2(x)/2+log^2(1-1/x)/4+pi^2/3 */ 342 { 343 int k; 344 mpfr_exp_t expo_l; 345 mpfr_t s, u, xx; 346 347 if (mpfr_cmp_ui (x, 38) >= 0) 348 { 349 inexact = mpfr_li2_asympt_pos (y, x, rnd_mode); 350 if (inexact != 0) 351 goto end_of_case_gt2; 352 } 353 354 mpfr_init2 (u, m); 355 mpfr_init2 (s, m); 356 mpfr_init2 (xx, m); 357 358 MPFR_ZIV_INIT (loop, m); 359 for (;;) 360 { 361 mpfr_ui_div (xx, 1, x, MPFR_RNDN); 362 mpfr_neg (xx, xx, MPFR_RNDN); 363 mpfr_log1p (u, xx, MPFR_RNDD); 364 mpfr_neg (u, u, MPFR_RNDU); /* u = -log(1-1/x) */ 365 expo_l = MPFR_GET_EXP (u); 366 k = li2_series (s, u, MPFR_RNDN); 367 mpfr_neg (s, s, MPFR_RNDN); 368 err = MPFR_INT_CEIL_LOG2 (k + 1) + 1; /* error(s) <= 2^err ulp(s) */ 369 370 mpfr_sqr (u, u, MPFR_RNDN); 371 mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u= log^2(1-1/x)/4 */ 372 mpfr_add (s, s, u, MPFR_RNDN); 373 err = 374 MAX (err, 375 3 + MAX (1, -expo_l) + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 376 err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 377 err += MPFR_GET_EXP (s); 378 379 mpfr_log (u, x, MPFR_RNDU); 380 mpfr_sqr (u, u, MPFR_RNDN); 381 mpfr_div_2ui (u, u, 1, MPFR_RNDN); /* u = log^2(x)/2 */ 382 mpfr_sub (s, s, u, MPFR_RNDN); 383 err = MAX (err, 3 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 384 err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 385 err += MPFR_GET_EXP (s); 386 387 mpfr_const_pi (u, MPFR_RNDU); 388 mpfr_sqr (u, u, MPFR_RNDN); 389 mpfr_div_ui (u, u, 3, MPFR_RNDN); /* u = pi^2/3 */ 390 mpfr_add (s, s, u, MPFR_RNDN); 391 err = MAX (err, 2) - MPFR_GET_EXP (s); 392 err = 2 + MAX (-1, err); /* error(s) <= 2^err ulp(s) */ 393 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 394 break; 395 396 MPFR_ZIV_NEXT (loop, m); 397 mpfr_set_prec (u, m); 398 mpfr_set_prec (s, m); 399 mpfr_set_prec (xx, m); 400 } 401 MPFR_ZIV_FREE (loop); 402 inexact = mpfr_set (y, s, rnd_mode); 403 mpfr_clears (s, u, xx, (mpfr_ptr) 0); 404 405 end_of_case_gt2: 406 MPFR_SAVE_EXPO_FREE (expo); 407 return mpfr_check_range (y, inexact, rnd_mode); 408 } 409 else if (mpfr_cmp_ui (x, 1) > 0) 410 /* 2 > x > 1: Li2(x) = S(log(x))+log^2(x)/4-log(x)log(x-1)+pi^2/6 */ 411 { 412 int k; 413 mpfr_exp_t e1, e2; 414 mpfr_t s, u, v, xx; 415 mpfr_init2 (s, m); 416 mpfr_init2 (u, m); 417 mpfr_init2 (v, m); 418 mpfr_init2 (xx, m); 419 420 MPFR_ZIV_INIT (loop, m); 421 for (;;) 422 { 423 mpfr_log (v, x, MPFR_RNDU); 424 k = li2_series (s, v, MPFR_RNDN); 425 e1 = MPFR_GET_EXP (s); 426 427 mpfr_sqr (u, v, MPFR_RNDN); 428 mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(x)/4 */ 429 mpfr_add (s, s, u, MPFR_RNDN); 430 431 mpfr_sub_ui (xx, x, 1, MPFR_RNDN); 432 mpfr_log (u, xx, MPFR_RNDU); 433 e2 = MPFR_GET_EXP (u); 434 mpfr_mul (u, v, u, MPFR_RNDN); /* u = log(x) * log(x-1) */ 435 mpfr_sub (s, s, u, MPFR_RNDN); 436 437 mpfr_const_pi (u, MPFR_RNDU); 438 mpfr_sqr (u, u, MPFR_RNDN); 439 mpfr_div_ui (u, u, 6, MPFR_RNDN); /* u = pi^2/6 */ 440 mpfr_add (s, s, u, MPFR_RNDN); 441 /* error(s) <= (31 + (k+1) * 2^(1-e1) + 2^(1-e2)) ulp(s) 442 see algorithms.tex */ 443 err = MAX (MPFR_INT_CEIL_LOG2 (k + 1) + 1 - e1, 1 - e2); 444 err = 2 + MAX (5, err); 445 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 446 break; 447 448 MPFR_ZIV_NEXT (loop, m); 449 mpfr_set_prec (s, m); 450 mpfr_set_prec (u, m); 451 mpfr_set_prec (v, m); 452 mpfr_set_prec (xx, m); 453 } 454 MPFR_ZIV_FREE (loop); 455 inexact = mpfr_set (y, s, rnd_mode); 456 457 mpfr_clears (s, u, v, xx, (mpfr_ptr) 0); 458 MPFR_SAVE_EXPO_FREE (expo); 459 return mpfr_check_range (y, inexact, rnd_mode); 460 } 461 else if (mpfr_cmp_ui_2exp (x, 1, -1) > 0) /* 1/2 < x < 1 */ 462 /* 1 > x > 1/2: Li2(x) = -S(-log(x))+log^2(x)/4-log(x)log(1-x)+pi^2/6 */ 463 { 464 int k; 465 mpfr_t s, u, v, xx; 466 mpfr_init2 (s, m); 467 mpfr_init2 (u, m); 468 mpfr_init2 (v, m); 469 mpfr_init2 (xx, m); 470 471 472 MPFR_ZIV_INIT (loop, m); 473 for (;;) 474 { 475 mpfr_log (u, x, MPFR_RNDD); 476 mpfr_neg (u, u, MPFR_RNDN); 477 k = li2_series (s, u, MPFR_RNDN); 478 mpfr_neg (s, s, MPFR_RNDN); 479 err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s); 480 481 mpfr_ui_sub (xx, 1, x, MPFR_RNDN); 482 mpfr_log (v, xx, MPFR_RNDU); 483 mpfr_mul (v, v, u, MPFR_RNDN); /* v = - log(x) * log(1-x) */ 484 mpfr_add (s, s, v, MPFR_RNDN); 485 err = MAX (err, 1 - MPFR_GET_EXP (v)); 486 err = 2 + MAX (3, err) - MPFR_GET_EXP (s); 487 488 mpfr_sqr (u, u, MPFR_RNDN); 489 mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(x)/4 */ 490 mpfr_add (s, s, u, MPFR_RNDN); 491 err = MAX (err, 2 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s); 492 err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 493 494 mpfr_const_pi (u, MPFR_RNDU); 495 mpfr_sqr (u, u, MPFR_RNDN); 496 mpfr_div_ui (u, u, 6, MPFR_RNDN); /* u = pi^2/6 */ 497 mpfr_add (s, s, u, MPFR_RNDN); 498 err = MAX (err, 3) - MPFR_GET_EXP (s); 499 err = 2 + MAX (-1, err); 500 501 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 502 break; 503 504 MPFR_ZIV_NEXT (loop, m); 505 mpfr_set_prec (s, m); 506 mpfr_set_prec (u, m); 507 mpfr_set_prec (v, m); 508 mpfr_set_prec (xx, m); 509 } 510 MPFR_ZIV_FREE (loop); 511 inexact = mpfr_set (y, s, rnd_mode); 512 513 mpfr_clears (s, u, v, xx, (mpfr_ptr) 0); 514 MPFR_SAVE_EXPO_FREE (expo); 515 return mpfr_check_range (y, inexact, rnd_mode); 516 } 517 else if (mpfr_cmp_si (x, -1) >= 0) 518 /* 0 > x >= -1: Li2(x) = -S(log(1-x))-log^2(1-x)/4 */ 519 { 520 int k; 521 mpfr_exp_t expo_l; 522 mpfr_t s, u, xx; 523 mpfr_init2 (s, m); 524 mpfr_init2 (u, m); 525 mpfr_init2 (xx, m); 526 527 MPFR_ZIV_INIT (loop, m); 528 for (;;) 529 { 530 mpfr_neg (xx, x, MPFR_RNDN); 531 mpfr_log1p (u, xx, MPFR_RNDN); 532 k = li2_series (s, u, MPFR_RNDN); 533 mpfr_neg (s, s, MPFR_RNDN); 534 expo_l = MPFR_GET_EXP (u); 535 err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s); 536 537 mpfr_sqr (u, u, MPFR_RNDN); 538 mpfr_div_2ui (u, u, 2, MPFR_RNDN); /* u = log^2(1-x)/4 */ 539 mpfr_sub (s, s, u, MPFR_RNDN); 540 err = MAX (err, - expo_l); 541 err = 2 + MAX (err, 3); 542 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 543 break; 544 545 MPFR_ZIV_NEXT (loop, m); 546 mpfr_set_prec (s, m); 547 mpfr_set_prec (u, m); 548 mpfr_set_prec (xx, m); 549 } 550 MPFR_ZIV_FREE (loop); 551 inexact = mpfr_set (y, s, rnd_mode); 552 553 mpfr_clears (s, u, xx, (mpfr_ptr) 0); 554 MPFR_SAVE_EXPO_FREE (expo); 555 return mpfr_check_range (y, inexact, rnd_mode); 556 } 557 else 558 /* x < -1: Li2(x) 559 = S(log(1-1/x))-log^2(-x)/4-log(1-x)log(-x)/2+log^2(1-x)/4-pi^2/6 */ 560 { 561 int k; 562 mpfr_t s, u, v, w, xx; 563 564 if (mpfr_cmp_si (x, -7) <= 0) 565 { 566 inexact = mpfr_li2_asympt_neg (y, x, rnd_mode); 567 if (inexact != 0) 568 goto end_of_case_ltm1; 569 } 570 571 mpfr_init2 (s, m); 572 mpfr_init2 (u, m); 573 mpfr_init2 (v, m); 574 mpfr_init2 (w, m); 575 mpfr_init2 (xx, m); 576 577 MPFR_ZIV_INIT (loop, m); 578 for (;;) 579 { 580 mpfr_ui_div (xx, 1, x, MPFR_RNDN); 581 mpfr_neg (xx, xx, MPFR_RNDN); 582 mpfr_log1p (u, xx, MPFR_RNDN); 583 k = li2_series (s, u, MPFR_RNDN); 584 585 mpfr_ui_sub (xx, 1, x, MPFR_RNDN); 586 mpfr_log (u, xx, MPFR_RNDU); 587 mpfr_neg (xx, x, MPFR_RNDN); 588 mpfr_log (v, xx, MPFR_RNDU); 589 mpfr_mul (w, v, u, MPFR_RNDN); 590 mpfr_div_2ui (w, w, 1, MPFR_RNDN); /* w = log(-x) * log(1-x) / 2 */ 591 mpfr_sub (s, s, w, MPFR_RNDN); 592 err = 1 + MAX (3, MPFR_INT_CEIL_LOG2 (k+1) + 1 - MPFR_GET_EXP (s)) 593 + MPFR_GET_EXP (s); 594 595 mpfr_sqr (w, v, MPFR_RNDN); 596 mpfr_div_2ui (w, w, 2, MPFR_RNDN); /* w = log^2(-x) / 4 */ 597 mpfr_sub (s, s, w, MPFR_RNDN); 598 err = MAX (err, 3 + MPFR_GET_EXP(w)) - MPFR_GET_EXP (s); 599 err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 600 601 mpfr_sqr (w, u, MPFR_RNDN); 602 mpfr_div_2ui (w, w, 2, MPFR_RNDN); /* w = log^2(1-x) / 4 */ 603 mpfr_add (s, s, w, MPFR_RNDN); 604 err = MAX (err, 3 + MPFR_GET_EXP (w)) - MPFR_GET_EXP (s); 605 err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 606 607 mpfr_const_pi (w, MPFR_RNDU); 608 mpfr_sqr (w, w, MPFR_RNDN); 609 mpfr_div_ui (w, w, 6, MPFR_RNDN); /* w = pi^2 / 6 */ 610 mpfr_sub (s, s, w, MPFR_RNDN); 611 err = MAX (err, 3) - MPFR_GET_EXP (s); 612 err = 2 + MAX (-1, err) + MPFR_GET_EXP (s); 613 614 if (MPFR_CAN_ROUND (s, (mpfr_exp_t) m - err, yp, rnd_mode)) 615 break; 616 617 MPFR_ZIV_NEXT (loop, m); 618 mpfr_set_prec (s, m); 619 mpfr_set_prec (u, m); 620 mpfr_set_prec (v, m); 621 mpfr_set_prec (w, m); 622 mpfr_set_prec (xx, m); 623 } 624 MPFR_ZIV_FREE (loop); 625 inexact = mpfr_set (y, s, rnd_mode); 626 mpfr_clears (s, u, v, w, xx, (mpfr_ptr) 0); 627 628 end_of_case_ltm1: 629 MPFR_SAVE_EXPO_FREE (expo); 630 return mpfr_check_range (y, inexact, rnd_mode); 631 } 632 633 MPFR_ASSERTN (0); /* should never reach this point */ 634 } 635