contrib/mpfr/li2.c

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