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