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           if (MPFR_IS_ZERO(u))
 337             goto next_m;
 338           mpfr_neg (u, u, GMP_RNDN);    /* u = -log(1-x) */
 339           expo_l = MPFR_GET_EXP (u);
 340           k = li2_series (s, u, GMP_RNDU);
 341           err = 1 + MPFR_INT_CEIL_LOG2 (k + 1);
 342
 343           mpfr_sqr (u, u, GMP_RNDU);
 344           mpfr_div_2ui (u, u, 2, GMP_RNDU);     /* u = log^2(1-x) / 4 */
 345           mpfr_sub (s, s, u, GMP_RNDN);
 346
 347           /* error(s) <= (0.5 + 2^(d-EXP(s))
 348              + 2^(3 + MAX(1, - expo_l) - EXP(s))) ulp(s) */
 349           err = MAX (err, MAX (1, - expo_l) - 1) - MPFR_GET_EXP (s);
 350           err = 2 + MAX (-1, err);
 351           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 352             break;
 353
 354         next_m:
 355           MPFR_ZIV_NEXT (loop, m);
 356           mpfr_set_prec (u, m);
 357           mpfr_set_prec (s, m);
 358         }
 359       MPFR_ZIV_FREE (loop);
 360       inexact = mpfr_set (y, s, rnd_mode);
 361
 362       mpfr_clear (u);
 363       mpfr_clear (s);
 364       MPFR_SAVE_EXPO_FREE (expo);
 365       return mpfr_check_range (y, inexact, rnd_mode);
 366     }
 367   else if (!mpfr_cmp_ui (x, 1))
 368     /* Li2(1)= pi^2 / 6 */
 369     {
 370       mpfr_t u;
 371       mpfr_init2 (u, m);
 372
 373       MPFR_ZIV_INIT (loop, m);
 374       for (;;)
 375         {
 376           mpfr_const_pi (u, GMP_RNDU);
 377           mpfr_sqr (u, u, GMP_RNDN);
 378           mpfr_div_ui (u, u, 6, GMP_RNDN);
 379
 380           err = m - 4;          /* error(u) <= 19/2 ulp(u) */
 381           if (MPFR_CAN_ROUND (u, err, yp, rnd_mode))
 382             break;
 383
 384           MPFR_ZIV_NEXT (loop, m);
 385           mpfr_set_prec (u, m);
 386         }
 387       MPFR_ZIV_FREE (loop);
 388       inexact = mpfr_set (y, u, rnd_mode);
 389
 390       mpfr_clear (u);
 391       MPFR_SAVE_EXPO_FREE (expo);
 392       return mpfr_check_range (y, inexact, rnd_mode);
 393     }
 394   else if (mpfr_cmp_ui (x, 2) >= 0)
 395     /* x >= 2: Li2(x) = -S(-log(1-1/x))-log^2(x)/2+log^2(1-1/x)/4+pi^2/3 */
 396     {
 397       int k;
 398       mp_exp_t expo_l;
 399       mpfr_t s, u, xx;
 400
 401       if (mpfr_cmp_ui (x, 38) >= 0)
 402         {
 403           inexact = mpfr_li2_asympt_pos (y, x, rnd_mode);
 404           if (inexact != 0)
 405             goto end_of_case_gt2;
 406         }
 407
 408       mpfr_init2 (u, m);
 409       mpfr_init2 (s, m);
 410       mpfr_init2 (xx, m);
 411
 412       MPFR_ZIV_INIT (loop, m);
 413       for (;;)
 414         {
 415           mpfr_ui_div (xx, 1, x, GMP_RNDN);
 416           mpfr_neg (xx, xx, GMP_RNDN);
 417           mpfr_log1p (u, xx, GMP_RNDD);
 418           mpfr_neg (u, u, GMP_RNDU);    /* u = -log(1-1/x) */
 419           expo_l = MPFR_GET_EXP (u);
 420           k = li2_series (s, u, GMP_RNDN);
 421           mpfr_neg (s, s, GMP_RNDN);
 422           err = MPFR_INT_CEIL_LOG2 (k + 1) + 1; /* error(s) <= 2^err ulp(s) */
 423
 424           mpfr_sqr (u, u, GMP_RNDN);
 425           mpfr_div_2ui (u, u, 2, GMP_RNDN);     /* u= log^2(1-1/x)/4 */
 426           mpfr_add (s, s, u, GMP_RNDN);
 427           err =
 428             MAX (err,
 429                  3 + MAX (1, -expo_l) + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s);
 430           err = 2 + MAX (-1, err);      /* error(s) <= 2^err ulp(s) */
 431           err += MPFR_GET_EXP (s);
 432
 433           mpfr_log (u, x, GMP_RNDU);
 434           mpfr_sqr (u, u, GMP_RNDN);
 435           mpfr_div_2ui (u, u, 1, GMP_RNDN);     /* u = log^2(x)/2 */
 436           mpfr_sub (s, s, u, GMP_RNDN);
 437           err = MAX (err, 3 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s);
 438           err = 2 + MAX (-1, err);      /* error(s) <= 2^err ulp(s) */
 439           err += MPFR_GET_EXP (s);
 440
 441           mpfr_const_pi (u, GMP_RNDU);
 442           mpfr_sqr (u, u, GMP_RNDN);
 443           mpfr_div_ui (u, u, 3, GMP_RNDN);      /* u = pi^2/3 */
 444           mpfr_add (s, s, u, GMP_RNDN);
 445           err = MAX (err, 2) - MPFR_GET_EXP (s);
 446           err = 2 + MAX (-1, err);      /* error(s) <= 2^err ulp(s) */
 447           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 448             break;
 449
 450           MPFR_ZIV_NEXT (loop, m);
 451           mpfr_set_prec (u, m);
 452           mpfr_set_prec (s, m);
 453           mpfr_set_prec (xx, m);
 454         }
 455       MPFR_ZIV_FREE (loop);
 456       inexact = mpfr_set (y, s, rnd_mode);
 457       mpfr_clears (s, u, xx, (mpfr_ptr) 0);
 458
 459     end_of_case_gt2:
 460       MPFR_SAVE_EXPO_FREE (expo);
 461       return mpfr_check_range (y, inexact, rnd_mode);
 462     }
 463   else if (mpfr_cmp_ui (x, 1) > 0)
 464     /* 2 > x > 1: Li2(x) = S(log(x))+log^2(x)/4-log(x)log(x-1)+pi^2/6 */
 465     {
 466       int k;
 467       mp_exp_t e1, e2;
 468       mpfr_t s, u, v, xx;
 469       mpfr_init2 (s, m);
 470       mpfr_init2 (u, m);
 471       mpfr_init2 (v, m);
 472       mpfr_init2 (xx, m);
 473
 474       MPFR_ZIV_INIT (loop, m);
 475       for (;;)
 476         {
 477           mpfr_log (v, x, GMP_RNDU);
 478           k = li2_series (s, v, GMP_RNDN);
 479           e1 = MPFR_GET_EXP (s);
 480
 481           mpfr_sqr (u, v, GMP_RNDN);
 482           mpfr_div_2ui (u, u, 2, GMP_RNDN);     /* u = log^2(x)/4 */
 483           mpfr_add (s, s, u, GMP_RNDN);
 484
 485           mpfr_sub_ui (xx, x, 1, GMP_RNDN);
 486           mpfr_log (u, xx, GMP_RNDU);
 487           e2 = MPFR_GET_EXP (u);
 488           mpfr_mul (u, v, u, GMP_RNDN); /* u = log(x) * log(x-1) */
 489           mpfr_sub (s, s, u, GMP_RNDN);
 490
 491           mpfr_const_pi (u, GMP_RNDU);
 492           mpfr_sqr (u, u, GMP_RNDN);
 493           mpfr_div_ui (u, u, 6, GMP_RNDN);      /* u = pi^2/6 */
 494           mpfr_add (s, s, u, GMP_RNDN);
 495           /* error(s) <= (31 + (k+1) * 2^(1-e1) + 2^(1-e2)) ulp(s)
 496              see algorithms.tex */
 497           err = MAX (MPFR_INT_CEIL_LOG2 (k + 1) + 1 - e1, 1 - e2);
 498           err = 2 + MAX (5, err);
 499           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 500             break;
 501
 502           MPFR_ZIV_NEXT (loop, m);
 503           mpfr_set_prec (s, m);
 504           mpfr_set_prec (u, m);
 505           mpfr_set_prec (v, m);
 506           mpfr_set_prec (xx, m);
 507         }
 508       MPFR_ZIV_FREE (loop);
 509       inexact = mpfr_set (y, s, rnd_mode);
 510
 511       mpfr_clears (s, u, v, xx, (mpfr_ptr) 0);
 512       MPFR_SAVE_EXPO_FREE (expo);
 513       return mpfr_check_range (y, inexact, rnd_mode);
 514     }
 515   else if (mpfr_cmp_ui_2exp (x, 1, -1) > 0) /*  1/2 < x < 1 */
 516     /* 1 > x > 1/2: Li2(x) = -S(-log(x))+log^2(x)/4-log(x)log(1-x)+pi^2/6 */
 517     {
 518       int k;
 519       mpfr_t s, u, v, xx;
 520       mpfr_init2 (s, m);
 521       mpfr_init2 (u, m);
 522       mpfr_init2 (v, m);
 523       mpfr_init2 (xx, m);
 524
 525
 526       MPFR_ZIV_INIT (loop, m);
 527       for (;;)
 528         {
 529           mpfr_log (u, x, GMP_RNDD);
 530           mpfr_neg (u, u, GMP_RNDN);
 531           k = li2_series (s, u, GMP_RNDN);
 532           mpfr_neg (s, s, GMP_RNDN);
 533           err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s);
 534
 535           mpfr_ui_sub (xx, 1, x, GMP_RNDN);
 536           mpfr_log (v, xx, GMP_RNDU);
 537           mpfr_mul (v, v, u, GMP_RNDN); /* v = - log(x) * log(1-x) */
 538           mpfr_add (s, s, v, GMP_RNDN);
 539           err = MAX (err, 1 - MPFR_GET_EXP (v));
 540           err = 2 + MAX (3, err) - MPFR_GET_EXP (s);
 541
 542           mpfr_sqr (u, u, GMP_RNDN);
 543           mpfr_div_2ui (u, u, 2, GMP_RNDN);     /* u = log^2(x)/4 */
 544           mpfr_add (s, s, u, GMP_RNDN);
 545           err = MAX (err, 2 + MPFR_GET_EXP (u)) - MPFR_GET_EXP (s);
 546           err = 2 + MAX (-1, err) + MPFR_GET_EXP (s);
 547
 548           mpfr_const_pi (u, GMP_RNDU);
 549           mpfr_sqr (u, u, GMP_RNDN);
 550           mpfr_div_ui (u, u, 6, GMP_RNDN);      /* u = pi^2/6 */
 551           mpfr_add (s, s, u, GMP_RNDN);
 552           err = MAX (err, 3) - MPFR_GET_EXP (s);
 553           err = 2 + MAX (-1, err);
 554
 555           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 556             break;
 557
 558           MPFR_ZIV_NEXT (loop, m);
 559           mpfr_set_prec (s, m);
 560           mpfr_set_prec (u, m);
 561           mpfr_set_prec (v, m);
 562           mpfr_set_prec (xx, m);
 563         }
 564       MPFR_ZIV_FREE (loop);
 565       inexact = mpfr_set (y, s, rnd_mode);
 566
 567       mpfr_clears (s, u, v, xx, (mpfr_ptr) 0);
 568       MPFR_SAVE_EXPO_FREE (expo);
 569       return mpfr_check_range (y, inexact, rnd_mode);
 570     }
 571   else if (mpfr_cmp_si (x, -1) >= 0)
 572     /* 0 > x >= -1: Li2(x) = -S(log(1-x))-log^2(1-x)/4 */
 573     {
 574       int k;
 575       mp_exp_t expo_l;
 576       mpfr_t s, u, xx;
 577       mpfr_init2 (s, m);
 578       mpfr_init2 (u, m);
 579       mpfr_init2 (xx, m);
 580
 581       MPFR_ZIV_INIT (loop, m);
 582       for (;;)
 583         {
 584           mpfr_neg (xx, x, GMP_RNDN);
 585           mpfr_log1p (u, xx, GMP_RNDN);
 586           k = li2_series (s, u, GMP_RNDN);
 587           mpfr_neg (s, s, GMP_RNDN);
 588           expo_l = MPFR_GET_EXP (u);
 589           err = 1 + MPFR_INT_CEIL_LOG2 (k + 1) - MPFR_GET_EXP (s);
 590
 591           mpfr_sqr (u, u, GMP_RNDN);
 592           mpfr_div_2ui (u, u, 2, GMP_RNDN);     /* u = log^2(1-x)/4 */
 593           mpfr_sub (s, s, u, GMP_RNDN);
 594           err = MAX (err, - expo_l);
 595           err = 2 + MAX (err, 3);
 596           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 597             break;
 598
 599           MPFR_ZIV_NEXT (loop, m);
 600           mpfr_set_prec (s, m);
 601           mpfr_set_prec (u, m);
 602           mpfr_set_prec (xx, m);
 603         }
 604       MPFR_ZIV_FREE (loop);
 605       inexact = mpfr_set (y, s, rnd_mode);
 606
 607       mpfr_clears (s, u, xx, (mpfr_ptr) 0);
 608       MPFR_SAVE_EXPO_FREE (expo);
 609       return mpfr_check_range (y, inexact, rnd_mode);
 610     }
 611   else
 612     /* x < -1: Li2(x)
 613        = S(log(1-1/x))-log^2(-x)/4-log(1-x)log(-x)/2+log^2(1-x)/4-pi^2/6 */
 614     {
 615       int k;
 616       mpfr_t s, u, v, w, xx;
 617
 618       if (mpfr_cmp_si (x, -7) <= 0)
 619         {
 620           inexact = mpfr_li2_asympt_neg (y, x, rnd_mode);
 621           if (inexact != 0)
 622             goto end_of_case_ltm1;
 623         }
 624
 625       mpfr_init2 (s, m);
 626       mpfr_init2 (u, m);
 627       mpfr_init2 (v, m);
 628       mpfr_init2 (w, m);
 629       mpfr_init2 (xx, m);
 630
 631       MPFR_ZIV_INIT (loop, m);
 632       for (;;)
 633         {
 634           mpfr_ui_div (xx, 1, x, GMP_RNDN);
 635           mpfr_neg (xx, xx, GMP_RNDN);
 636           mpfr_log1p (u, xx, GMP_RNDN);
 637           k = li2_series (s, u, GMP_RNDN);
 638
 639           mpfr_ui_sub (xx, 1, x, GMP_RNDN);
 640           mpfr_log (u, xx, GMP_RNDU);
 641           mpfr_neg (xx, x, GMP_RNDN);
 642           mpfr_log (v, xx, GMP_RNDU);
 643           mpfr_mul (w, v, u, GMP_RNDN);
 644           mpfr_div_2ui (w, w, 1, GMP_RNDN);  /* w = log(-x) * log(1-x) / 2 */
 645           mpfr_sub (s, s, w, GMP_RNDN);
 646           err = 1 + MAX (3, MPFR_INT_CEIL_LOG2 (k+1) + 1  - MPFR_GET_EXP (s))
 647             + MPFR_GET_EXP (s);
 648
 649           mpfr_sqr (w, v, GMP_RNDN);
 650           mpfr_div_2ui (w, w, 2, GMP_RNDN);  /* w = log^2(-x) / 4 */
 651           mpfr_sub (s, s, w, GMP_RNDN);
 652           err = MAX (err, 3 + MPFR_GET_EXP(w)) - MPFR_GET_EXP (s);
 653           err = 2 + MAX (-1, err) + MPFR_GET_EXP (s);
 654
 655           mpfr_sqr (w, u, GMP_RNDN);
 656           mpfr_div_2ui (w, w, 2, GMP_RNDN);     /* w = log^2(1-x) / 4 */
 657           mpfr_add (s, s, w, GMP_RNDN);
 658           err = MAX (err, 3 + MPFR_GET_EXP (w)) - MPFR_GET_EXP (s);
 659           err = 2 + MAX (-1, err) + MPFR_GET_EXP (s);
 660
 661           mpfr_const_pi (w, GMP_RNDU);
 662           mpfr_sqr (w, w, GMP_RNDN);
 663           mpfr_div_ui (w, w, 6, GMP_RNDN);      /* w = pi^2 / 6 */
 664           mpfr_sub (s, s, w, GMP_RNDN);
 665           err = MAX (err, 3) - MPFR_GET_EXP (s);
 666           err = 2 + MAX (-1, err) + MPFR_GET_EXP (s);
 667
 668           if (MPFR_CAN_ROUND (s, (mp_exp_t) m - err, yp, rnd_mode))
 669             break;
 670
 671           MPFR_ZIV_NEXT (loop, m);
 672           mpfr_set_prec (s, m);
 673           mpfr_set_prec (u, m);
 674           mpfr_set_prec (v, m);
 675           mpfr_set_prec (w, m);
 676           mpfr_set_prec (xx, m);
 677         }
 678       MPFR_ZIV_FREE (loop);
 679       inexact = mpfr_set (y, s, rnd_mode);
 680       mpfr_clears (s, u, v, w, xx, (mpfr_ptr) 0);
 681
 682     end_of_case_ltm1:
 683       MPFR_SAVE_EXPO_FREE (expo);
 684       return mpfr_check_range (y, inexact, rnd_mode);
 685     }
 686
 687   MPFR_ASSERTN (0);             /* should never reach this point */
 688 }