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
li2_series(mpfr_t sum,mpfr_srcptr z,mpfr_rnd_t rnd_mode)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
mpfr_li2_asympt_pos(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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
mpfr_li2_asympt_neg(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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
mpfr_li2(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)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