contrib/mpfr/src/li2.c

   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
  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
 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
 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
 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 }