1 /* mpfr_li2 -- Dilogarithm.
2
3 Copyright 2007-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 #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_ptr 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 an additional check that z
44 is (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[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
225 ("y[%Pd]=%.*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 }