1/* mpfr_li2 -- Dilogarithm. 2 3Copyright 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel 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 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/* 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*/ 33static int 34li2_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*/ 134static int 135mpfr_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*/ 180static int 181mpfr_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 */ 217int 218mpfr_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