contrib/mpfr/src/eint.c

   1 /* mpfr_eint, mpfr_eint1 -- the exponential integral
   2
   3 Copyright 2005, 2006, 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 /* eint1(x) = -gamma - log(x) - sum((-1)^k*z^k/k/k!, k=1..infinity) for x > 0
  27             = - eint(-x) for x < 0
  28    where
  29    eint (x) = gamma + log(x) + sum(z^k/k/k!, k=1..infinity) for x > 0
  30    eint (x) is undefined for x < 0.
  31 */
  32
  33 /* compute in y an approximation of sum(x^k/k/k!, k=1..infinity),
  34    and return e such that the absolute error is bound by 2^e ulp(y) */
  35 static mpfr_exp_t
  36 mpfr_eint_aux (mpfr_t y, mpfr_srcptr x)
  37 {
  38   mpfr_t eps; /* dynamic (absolute) error bound on t */
  39   mpfr_t erru, errs;
  40   mpz_t m, s, t, u;
  41   mpfr_exp_t e, sizeinbase;
  42   mpfr_prec_t w = MPFR_PREC(y);
  43   unsigned long k;
  44   MPFR_GROUP_DECL (group);
  45
  46   /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x)
  47      where |R(x)| <= (x/2)^2/(1-x/2) <= 2*(x/2)^2
  48      thus |R(x)/x| <= |x|/2
  49      thus if |x| <= 2^(-PREC(y)) we have |S - o(x)| <= ulp(y) */
  50
  51   if (MPFR_GET_EXP(x) <= - (mpfr_exp_t) w)
  52     {
  53       mpfr_set (y, x, MPFR_RNDN);
  54       return 0;
  55     }
  56
  57   mpz_init (s); /* initializes to 0 */
  58   mpz_init (t);
  59   mpz_init (u);
  60   mpz_init (m);
  61   MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs);
  62   e = mpfr_get_z_2exp (m, x); /* x = m * 2^e */
  63   MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x));
  64   if (MPFR_PREC (x) > w)
  65     {
  66       e += MPFR_PREC (x) - w;
  67       mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w);
  68     }
  69   /* remove trailing zeroes from m: this will speed up much cases where
  70      x is a small integer divided by a power of 2 */
  71   k = mpz_scan1 (m, 0);
  72   mpz_tdiv_q_2exp (m, m, k);
  73   e += k;
  74   /* initialize t to 2^w */
  75   mpz_set_ui (t, 1);
  76   mpz_mul_2exp (t, t, w);
  77   mpfr_set_ui (eps, 0, MPFR_RNDN); /* eps[0] = 0 */
  78   mpfr_set_ui (errs, 0, MPFR_RNDN);
  79   for (k = 1;; k++)
  80     {
  81       /* let eps[k] be the absolute error on t[k]:
  82          since t[k] = trunc(t[k-1]*m*2^e/k), we have
  83          eps[k+1] <= 1 + eps[k-1]*m*2^e/k + t[k-1]*m*2^(1-w)*2^e/k
  84                   =  1 + (eps[k-1] + t[k-1]*2^(1-w))*m*2^e/k
  85                   = 1 + (eps[k-1]*2^(w-1) + t[k-1])*2^(1-w)*m*2^e/k */
  86       mpfr_mul_2ui (eps, eps, w - 1, MPFR_RNDU);
  87       mpfr_add_z (eps, eps, t, MPFR_RNDU);
  88       MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
  89       mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, MPFR_RNDU);
  90       mpfr_div_ui (eps, eps, k, MPFR_RNDU);
  91       mpfr_add_ui (eps, eps, 1, MPFR_RNDU);
  92       mpz_mul (t, t, m);
  93       if (e < 0)
  94         mpz_tdiv_q_2exp (t, t, -e);
  95       else
  96         mpz_mul_2exp (t, t, e);
  97       mpz_tdiv_q_ui (t, t, k);
  98       mpz_tdiv_q_ui (u, t, k);
  99       mpz_add (s, s, u);
 100       /* the absolute error on u is <= 1 + eps[k]/k */
 101       mpfr_div_ui (erru, eps, k, MPFR_RNDU);
 102       mpfr_add_ui (erru, erru, 1, MPFR_RNDU);
 103       /* and that on s is the sum of all errors on u */
 104       mpfr_add (errs, errs, erru, MPFR_RNDU);
 105       /* we are done when t is smaller than errs */
 106       if (mpz_sgn (t) == 0)
 107         sizeinbase = 0;
 108       else
 109         MPFR_MPZ_SIZEINBASE2 (sizeinbase, t);
 110       if (sizeinbase < MPFR_GET_EXP (errs))
 111         break;
 112     }
 113   /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
 114      <= (|t|+eps)/k*|x|/(k-|x|) */
 115   mpz_abs (t, t);
 116   mpfr_add_z (eps, eps, t, MPFR_RNDU);
 117   mpfr_div_ui (eps, eps, k, MPFR_RNDU);
 118   mpfr_abs (erru, x, MPFR_RNDU); /* |x| */
 119   mpfr_mul (eps, eps, erru, MPFR_RNDU);
 120   mpfr_ui_sub (erru, k, erru, MPFR_RNDD);
 121   if (MPFR_IS_NEG (erru))
 122     {
 123       /* the truncated series does not converge, return fail */
 124       e = w;
 125     }
 126   else
 127     {
 128       mpfr_div (eps, eps, erru, MPFR_RNDU);
 129       mpfr_add (errs, errs, eps, MPFR_RNDU);
 130       mpfr_set_z (y, s, MPFR_RNDN);
 131       mpfr_div_2ui (y, y, w, MPFR_RNDN);
 132       /* errs was an absolute error bound on s. We must convert it to an error
 133          in terms of ulp(y). Since ulp(y) = 2^(EXP(y)-PREC(y)), we must
 134          divide the error by 2^(EXP(y)-PREC(y)), but since we divided also
 135          y by 2^w = 2^PREC(y), we must simply divide by 2^EXP(y). */
 136       e = MPFR_GET_EXP (errs) - MPFR_GET_EXP (y);
 137     }
 138   MPFR_GROUP_CLEAR (group);
 139   mpz_clear (s);
 140   mpz_clear (t);
 141   mpz_clear (u);
 142   mpz_clear (m);
 143   return e;
 144 }
 145
 146 /* Return in y an approximation of Ei(x) using the asymptotic expansion:
 147    Ei(x) = exp(x)/x * (1 + 1/x + 2/x^2 + ... + k!/x^k + ...)
 148    Assumes x >= PREC(y) * log(2).
 149    Returns the error bound in terms of ulp(y).
 150 */
 151 static mpfr_exp_t
 152 mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x)
 153 {
 154   mpfr_prec_t p = MPFR_PREC(y);
 155   mpfr_t invx, t, err;
 156   unsigned long k;
 157   mpfr_exp_t err_exp;
 158
 159   mpfr_init2 (t, p);
 160   mpfr_init2 (invx, p);
 161   mpfr_init2 (err, 31); /* error in ulps on y */
 162   mpfr_ui_div (invx, 1, x, MPFR_RNDN); /* invx = 1/x*(1+u) with |u|<=2^(1-p) */
 163   mpfr_set_ui (t, 1, MPFR_RNDN); /* exact */
 164   mpfr_set (y, t, MPFR_RNDN);
 165   mpfr_set_ui (err, 0, MPFR_RNDN);
 166   for (k = 1; MPFR_GET_EXP(t) + (mpfr_exp_t) p > MPFR_GET_EXP(y); k++)
 167     {
 168       mpfr_mul (t, t, invx, MPFR_RNDN); /* 2 more roundings */
 169       mpfr_mul_ui (t, t, k, MPFR_RNDN); /* 1 more rounding: t = k!/x^k*(1+u)^e
 170                                           with u=2^{-p} and |e| <= 3*k */
 171       /* we use the fact that |(1+u)^n-1| <= 2*|n*u| for |n*u| <= 1, thus
 172          the error on t is less than 6*k*2^{-p}*t <= 6*k*ulp(t) */
 173       /* err is in terms of ulp(y): transform it in terms of ulp(t) */
 174       mpfr_mul_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
 175       mpfr_add_ui (err, err, 6 * k, MPFR_RNDU);
 176       /* transform back in terms of ulp(y) */
 177       mpfr_div_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), MPFR_RNDU);
 178       mpfr_add (y, y, t, MPFR_RNDN);
 179     }
 180   /* add the truncation error bounded by ulp(y): 1 ulp */
 181   mpfr_mul (y, y, invx, MPFR_RNDN); /* err <= 2*err + 3/2 */
 182   mpfr_exp (t, x, MPFR_RNDN); /* err(t) <= 1/2*ulp(t) */
 183   mpfr_mul (y, y, t, MPFR_RNDN); /* again: err <= 2*err + 3/2 */
 184   mpfr_mul_2ui (err, err, 2, MPFR_RNDU);
 185   mpfr_add_ui (err, err, 8, MPFR_RNDU);
 186   err_exp = MPFR_GET_EXP(err);
 187   mpfr_clear (t);
 188   mpfr_clear (invx);
 189   mpfr_clear (err);
 190   return err_exp;
 191 }
 192
 193 int
 194 mpfr_eint (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
 195 {
 196   int inex;
 197   mpfr_t tmp, ump;
 198   mpfr_exp_t err, te;
 199   mpfr_prec_t prec;
 200   MPFR_SAVE_EXPO_DECL (expo);
 201   MPFR_ZIV_DECL (loop);
 202
 203   MPFR_LOG_FUNC (
 204     ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd),
 205     ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, inex));
 206
 207   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 208     {
 209       /* exp(NaN) = exp(-Inf) = NaN */
 210       if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x)))
 211         {
 212           MPFR_SET_NAN (y);
 213           MPFR_RET_NAN;
 214         }
 215       /* eint(+inf) = +inf */
 216       else if (MPFR_IS_INF (x))
 217         {
 218           MPFR_SET_INF(y);
 219           MPFR_SET_POS(y);
 220           MPFR_RET(0);
 221         }
 222       else /* eint(+/-0) = -Inf */
 223         {
 224           MPFR_SET_INF(y);
 225           MPFR_SET_NEG(y);
 226           mpfr_set_divby0 ();
 227           MPFR_RET(0);
 228         }
 229     }
 230
 231   /* eint(x) = NaN for x < 0 */
 232   if (MPFR_IS_NEG(x))
 233     {
 234       MPFR_SET_NAN (y);
 235       MPFR_RET_NAN;
 236     }
 237
 238   MPFR_SAVE_EXPO_MARK (expo);
 239
 240   /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
 241      Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
 242      then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
 243   mpfr_init2 (tmp, 64);
 244   mpfr_init2 (ump, 64);
 245   mpfr_log (tmp, x, MPFR_RNDU);
 246   mpfr_sub (ump, x, tmp, MPFR_RNDD);
 247   mpfr_const_log2 (tmp, MPFR_RNDU);
 248   mpfr_div (ump, ump, tmp, MPFR_RNDD);
 249   /* FIXME: We really need mpfr_set_exp_t and mpfr_cmpfr_exp_t functions. */
 250   MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
 251   if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
 252     {
 253       mpfr_clear (tmp);
 254       mpfr_clear (ump);
 255       MPFR_SAVE_EXPO_FREE (expo);
 256       return mpfr_overflow (y, rnd, 1);
 257     }
 258
 259   /* Init stuff */
 260   prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;
 261
 262   /* eint() has a root 0.37250741078136663446..., so if x is near,
 263      already take more bits */
 264   /* FIXME: do not use native floating-point here. */
 265   if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
 266     {
 267       double d;
 268       d = mpfr_get_d (x, MPFR_RNDN) - 0.37250741078136663;
 269       d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d);
 270       prec += -d;
 271     }
 272
 273   mpfr_set_prec (tmp, prec);
 274   mpfr_set_prec (ump, prec);
 275
 276   MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
 277   for (;;)                               /* Infinite loop */
 278     {
 279       /* We need that the smallest value of k!/x^k is smaller than 2^(-p).
 280          The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x
 281          for x>=1. */
 282       if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec +
 283                             0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
 284         err = mpfr_eint_asympt (tmp, x);
 285       else
 286         {
 287           err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
 288           te = MPFR_GET_EXP(tmp);
 289           mpfr_const_euler (ump, MPFR_RNDN); /* 0.577 -> EXP(ump)=0 */
 290           mpfr_add (tmp, tmp, ump, MPFR_RNDN);
 291           /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
 292              <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
 293              <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */
 294           err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp);
 295           err = MAX(0, err);
 296           te = MPFR_GET_EXP(tmp);
 297           mpfr_log (ump, x, MPFR_RNDN);
 298           mpfr_add (tmp, tmp, ump, MPFR_RNDN);
 299           /* same formula as above, except now EXP(ump) is not 0 */
 300           err += te + 1;
 301           if (MPFR_LIKELY (!MPFR_IS_ZERO (ump)))
 302             err = MAX (MPFR_GET_EXP (ump), err);
 303           err = MAX(0, err - MPFR_GET_EXP (tmp));
 304         }
 305       if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
 306         break;
 307       MPFR_ZIV_NEXT (loop, prec);        /* Increase used precision */
 308       mpfr_set_prec (tmp, prec);
 309       mpfr_set_prec (ump, prec);
 310     }
 311   MPFR_ZIV_FREE (loop);                  /* Free the ZivLoop Controler */
 312
 313   inex = mpfr_set (y, tmp, rnd);    /* Set y to the computed value */
 314   mpfr_clear (tmp);
 315   mpfr_clear (ump);
 316
 317   MPFR_SAVE_EXPO_FREE (expo);
 318   return mpfr_check_range (y, inex, rnd);
 319 }