1 /* mpfr_li2 -- Dilogarithm.
3 Copyright 2007-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
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
30 Return the maximum index of the truncature which is useful
31 for determinating the relative error.
34 li2_series (mpfr_t sum
, mpfr_srcptr z
, mpfr_rnd_t rnd_mode
)
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
46 MPFR_ASSERTD (MPFR_IS_STRICTPOS (z
));
47 MPFR_ASSERTD (mpfr_cmp_d (z
, 0.6953125) <= 0);
49 sump
= MPFR_PREC (sum
); /* target precision */
50 p
= sump
+ MPFR_INT_CEIL_LOG2 (sump
) + 4; /* the working precision */
56 B
= mpfr_bernoulli_internal ((mpz_t
*) 0, 0);
59 MPFR_ZIV_INIT (loop
, p
);
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
);
71 B
= mpfr_bernoulli_internal (B
, Bmax
++); /* B_2i*(2i+1)!, exact */
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 */
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) */
84 mpfr_add (s
, s
, w
, MPFR_RNDN
);
86 err
= MAX (err
+ se
, 5 * i
+ 8 + MPFR_GET_EXP (w
))
88 err
= 2 + MAX (-1, err
);
89 se
= MPFR_GET_EXP (s
);
90 if (MPFR_GET_EXP (w
) <= se
- (mpfr_exp_t
) p
)
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
))
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
);
107 MPFR_ZIV_FREE (loop
);
108 mpfr_set (sum
, s
, rnd_mode
);
113 (*__gmp_free_func
) (B
, Bmax
* sizeof (mpz_t
));
114 mpfr_clears (s
, u
, v
, w
, (mpfr_ptr
) 0);
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).
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.
135 mpfr_li2_asympt_pos (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
138 mpfr_prec_t w
= MPFR_PREC (y
) + 20;
141 MPFR_ASSERTN (mpfr_cmp_ui (x
, 38) >= 0);
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|
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).
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
);
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.
181 mpfr_li2_asympt_neg (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
184 mpfr_prec_t w
= MPFR_PREC (y
) + 20;
187 MPFR_ASSERTN (mpfr_cmp_si (x
, -7) <= 0);
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|
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).
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
);
215 /* Compute the real part of the dilogarithm defined by
216 Li2(x) = -\Int_{t=0}^x log(1-t)/t dt */
218 mpfr_li2 (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
223 MPFR_ZIV_DECL (loop
);
224 MPFR_SAVE_EXPO_DECL (expo
);
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
));
231 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
238 else if (MPFR_IS_INF (x
))
246 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
247 MPFR_SET_SAME_SIGN (y
, x
);
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) */
257 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y
, x
, -MPFR_GET_EXP (x
), 1, 1, rnd_mode
,
260 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y
, x
, -MPFR_GET_EXP (x
), 2, 0, rnd_mode
,
263 MPFR_SAVE_EXPO_MARK (expo
);
265 m
= yp
+ MPFR_INT_CEIL_LOG2 (yp
) + 13;
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 */
277 MPFR_ZIV_INIT (loop
, m
);
280 mpfr_ui_sub (u
, 1, x
, MPFR_RNDN
);
281 mpfr_log (u
, u
, MPFR_RNDU
);
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);
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
);
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
))
301 MPFR_ZIV_NEXT (loop
, m
);
302 mpfr_set_prec (u
, m
);
303 mpfr_set_prec (s
, m
);
305 MPFR_ZIV_FREE (loop
);
306 inexact
= mpfr_set (y
, s
, rnd_mode
);
310 MPFR_SAVE_EXPO_FREE (expo
);
311 return mpfr_check_range (y
, inexact
, rnd_mode
);
313 else if (!mpfr_cmp_ui (x
, 1))
314 /* Li2(1)= pi^2 / 6 */
319 MPFR_ZIV_INIT (loop
, m
);
322 mpfr_const_pi (u
, MPFR_RNDU
);
323 mpfr_sqr (u
, u
, MPFR_RNDN
);
324 mpfr_div_ui (u
, u
, 6, MPFR_RNDN
);
326 err
= m
- 4; /* error(u) <= 19/2 ulp(u) */
327 if (MPFR_CAN_ROUND (u
, err
, yp
, rnd_mode
))
330 MPFR_ZIV_NEXT (loop
, m
);
331 mpfr_set_prec (u
, m
);
333 MPFR_ZIV_FREE (loop
);
334 inexact
= mpfr_set (y
, u
, rnd_mode
);
337 MPFR_SAVE_EXPO_FREE (expo
);
338 return mpfr_check_range (y
, inexact
, rnd_mode
);
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 */
347 if (mpfr_cmp_ui (x
, 38) >= 0)
349 inexact
= mpfr_li2_asympt_pos (y
, x
, rnd_mode
);
351 goto end_of_case_gt2
;
358 MPFR_ZIV_INIT (loop
, m
);
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) */
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
);
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
);
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
);
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
))
396 MPFR_ZIV_NEXT (loop
, m
);
397 mpfr_set_prec (u
, m
);
398 mpfr_set_prec (s
, m
);
399 mpfr_set_prec (xx
, m
);
401 MPFR_ZIV_FREE (loop
);
402 inexact
= mpfr_set (y
, s
, rnd_mode
);
403 mpfr_clears (s
, u
, xx
, (mpfr_ptr
) 0);
406 MPFR_SAVE_EXPO_FREE (expo
);
407 return mpfr_check_range (y
, inexact
, rnd_mode
);
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 */
420 MPFR_ZIV_INIT (loop
, m
);
423 mpfr_log (v
, x
, MPFR_RNDU
);
424 k
= li2_series (s
, v
, MPFR_RNDN
);
425 e1
= MPFR_GET_EXP (s
);
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
);
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
);
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
))
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
);
454 MPFR_ZIV_FREE (loop
);
455 inexact
= mpfr_set (y
, s
, rnd_mode
);
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
);
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 */
472 MPFR_ZIV_INIT (loop
, m
);
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
);
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
);
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
);
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
);
501 if (MPFR_CAN_ROUND (s
, (mpfr_exp_t
) m
- err
, yp
, rnd_mode
))
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
);
510 MPFR_ZIV_FREE (loop
);
511 inexact
= mpfr_set (y
, s
, rnd_mode
);
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
);
517 else if (mpfr_cmp_si (x
, -1) >= 0)
518 /* 0 > x >= -1: Li2(x) = -S(log(1-x))-log^2(1-x)/4 */
527 MPFR_ZIV_INIT (loop
, m
);
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
);
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
))
545 MPFR_ZIV_NEXT (loop
, m
);
546 mpfr_set_prec (s
, m
);
547 mpfr_set_prec (u
, m
);
548 mpfr_set_prec (xx
, m
);
550 MPFR_ZIV_FREE (loop
);
551 inexact
= mpfr_set (y
, s
, rnd_mode
);
553 mpfr_clears (s
, u
, xx
, (mpfr_ptr
) 0);
554 MPFR_SAVE_EXPO_FREE (expo
);
555 return mpfr_check_range (y
, inexact
, rnd_mode
);
559 = S(log(1-1/x))-log^2(-x)/4-log(1-x)log(-x)/2+log^2(1-x)/4-pi^2/6 */
562 mpfr_t s
, u
, v
, w
, xx
;
564 if (mpfr_cmp_si (x
, -7) <= 0)
566 inexact
= mpfr_li2_asympt_neg (y
, x
, rnd_mode
);
568 goto end_of_case_ltm1
;
577 MPFR_ZIV_INIT (loop
, m
);
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
);
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
))
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
);
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
);
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
);
614 if (MPFR_CAN_ROUND (s
, (mpfr_exp_t
) m
- err
, yp
, rnd_mode
))
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
);
624 MPFR_ZIV_FREE (loop
);
625 inexact
= mpfr_set (y
, s
, rnd_mode
);
626 mpfr_clears (s
, u
, v
, w
, xx
, (mpfr_ptr
) 0);
629 MPFR_SAVE_EXPO_FREE (expo
);
630 return mpfr_check_range (y
, inexact
, rnd_mode
);
633 MPFR_RET_NEVER_GO_HERE ();