contrib/mpfr/eint.c

   1 /* mpfr_eint, mpfr_eint1 -- the exponential integral
   2
   3 Copyright 2005, 2006, 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 /* 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 mp_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   mp_exp_t e, sizeinbase;
  42   mp_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) <= - (mp_exp_t) w)
  52     {
  53       mpfr_set (y, x, GMP_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_exp (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, GMP_RNDN); /* eps[0] = 0 */
  78   mpfr_set_ui (errs, 0, GMP_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, GMP_RNDU);
  87       mpfr_add_z (eps, eps, t, GMP_RNDU);
  88       MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
  89       mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, GMP_RNDU);
  90       mpfr_div_ui (eps, eps, k, GMP_RNDU);
  91       mpfr_add_ui (eps, eps, 1, GMP_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, GMP_RNDU);
 102       mpfr_add_ui (erru, erru, 1, GMP_RNDU);
 103       /* and that on s is the sum of all errors on u */
 104       mpfr_add (errs, errs, erru, GMP_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, GMP_RNDU);
 117   mpfr_div_ui (eps, eps, k, GMP_RNDU);
 118   mpfr_abs (erru, x, GMP_RNDU); /* |x| */
 119   mpfr_mul (eps, eps, erru, GMP_RNDU);
 120   mpfr_ui_sub (erru, k, erru, GMP_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, GMP_RNDU);
 129       mpfr_add (errs, errs, eps, GMP_RNDU);
 130       mpfr_set_z (y, s, GMP_RNDN);
 131       mpfr_div_2ui (y, y, w, GMP_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 mp_exp_t
 152 mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x)
 153 {
 154   mp_prec_t p = MPFR_PREC(y);
 155   mpfr_t invx, t, err;
 156   unsigned long k;
 157   mp_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, GMP_RNDN); /* invx = 1/x*(1+u) with |u|<=2^(1-p) */
 163   mpfr_set_ui (t, 1, GMP_RNDN); /* exact */
 164   mpfr_set (y, t, GMP_RNDN);
 165   mpfr_set_ui (err, 0, GMP_RNDN);
 166   for (k = 1; MPFR_GET_EXP(t) + (mp_exp_t) p > MPFR_GET_EXP(y); k++)
 167     {
 168       mpfr_mul (t, t, invx, GMP_RNDN); /* 2 more roundings */
 169       mpfr_mul_ui (t, t, k, GMP_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), GMP_RNDU);
 175       mpfr_add_ui (err, err, 6 * k, GMP_RNDU);
 176       /* transform back in terms of ulp(y) */
 177       mpfr_div_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), GMP_RNDU);
 178       mpfr_add (y, y, t, GMP_RNDN);
 179     }
 180   /* add the truncation error bounded by ulp(y): 1 ulp */
 181   mpfr_mul (y, y, invx, GMP_RNDN); /* err <= 2*err + 3/2 */
 182   mpfr_exp (t, x, GMP_RNDN); /* err(t) <= 1/2*ulp(t) */
 183   mpfr_mul (y, y, t, GMP_RNDN); /* again: err <= 2*err + 3/2 */
 184   mpfr_mul_2ui (err, err, 2, GMP_RNDU);
 185   mpfr_add_ui (err, err, 8, GMP_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, mp_rnd_t rnd)
 195 {
 196   int inex;
 197   mpfr_t tmp, ump;
 198   mp_exp_t err, te;
 199   mp_prec_t prec;
 200   MPFR_SAVE_EXPO_DECL (expo);
 201   MPFR_ZIV_DECL (loop);
 202
 203   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
 204                  ("y[%#R]=%R inexact=%d", y, y, inex));
 205
 206   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 207     {
 208       /* exp(NaN) = exp(-Inf) = NaN */
 209       if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x)))
 210         {
 211           MPFR_SET_NAN (y);
 212           MPFR_RET_NAN;
 213         }
 214       /* eint(+inf) = +inf */
 215       else if (MPFR_IS_INF (x))
 216         {
 217           MPFR_SET_INF(y);
 218           MPFR_SET_POS(y);
 219           MPFR_RET(0);
 220         }
 221       else /* eint(+/-0) = -Inf */
 222         {
 223           MPFR_SET_INF(y);
 224           MPFR_SET_NEG(y);
 225           MPFR_RET(0);
 226         }
 227     }
 228
 229   /* eint(x) = NaN for x < 0 */
 230   if (MPFR_IS_NEG(x))
 231     {
 232       MPFR_SET_NAN (y);
 233       MPFR_RET_NAN;
 234     }
 235
 236   MPFR_SAVE_EXPO_MARK (expo);
 237
 238   /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
 239      Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
 240      then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
 241   mpfr_init2 (tmp, 64);
 242   mpfr_init2 (ump, 64);
 243   mpfr_log (tmp, x, GMP_RNDU);
 244   mpfr_sub (ump, x, tmp, GMP_RNDD);
 245   mpfr_const_log2 (tmp, GMP_RNDU);
 246   mpfr_div (ump, ump, tmp, GMP_RNDD);
 247   /* FIXME: We really need mpfr_set_exp_t and mpfr_cmp_exp_t functions. */
 248   MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
 249   if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
 250     {
 251       mpfr_clear (tmp);
 252       mpfr_clear (ump);
 253       MPFR_SAVE_EXPO_FREE (expo);
 254       return mpfr_overflow (y, rnd, 1);
 255     }
 256
 257   /* Init stuff */
 258   prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;
 259
 260   /* eint() has a root 0.37250741078136663446..., so if x is near,
 261      already take more bits */
 262   if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
 263     {
 264       double d;
 265       d = mpfr_get_d (x, GMP_RNDN) - 0.37250741078136663;
 266       d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d);
 267       prec += -d;
 268     }
 269
 270   mpfr_set_prec (tmp, prec);
 271   mpfr_set_prec (ump, prec);
 272
 273   MPFR_ZIV_INIT (loop, prec);            /* Initialize the ZivLoop controler */
 274   for (;;)                               /* Infinite loop */
 275     {
 276       /* We need that the smallest value of k!/x^k is smaller than 2^(-p).
 277          The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x
 278          for x>=1. */
 279       if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec +
 280                             0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
 281         err = mpfr_eint_asympt (tmp, x);
 282       else
 283         {
 284           err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
 285           te = MPFR_GET_EXP(tmp);
 286           mpfr_const_euler (ump, GMP_RNDN); /* 0.577 -> EXP(ump)=0 */
 287           mpfr_add (tmp, tmp, ump, GMP_RNDN);
 288           /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
 289              <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
 290              <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */
 291           err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp);
 292           err = MAX(0, err);
 293           te = MPFR_GET_EXP(tmp);
 294           mpfr_log (ump, x, GMP_RNDN);
 295           mpfr_add (tmp, tmp, ump, GMP_RNDN);
 296           /* same formula as above, except now EXP(ump) is not 0 */
 297           err += te + 1;
 298           if (MPFR_LIKELY (!MPFR_IS_ZERO (ump)))
 299             err = MAX (MPFR_GET_EXP (ump), err);
 300           err = MAX(0, err - MPFR_GET_EXP (tmp));
 301         }
 302       if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
 303         break;
 304       MPFR_ZIV_NEXT (loop, prec);        /* Increase used precision */
 305       mpfr_set_prec (tmp, prec);
 306       mpfr_set_prec (ump, prec);
 307     }
 308   MPFR_ZIV_FREE (loop);                  /* Free the ZivLoop Controler */
 309
 310   inex = mpfr_set (y, tmp, rnd);    /* Set y to the computed value */
 311   mpfr_clear (tmp);
 312   mpfr_clear (ump);
 313
 314   MPFR_SAVE_EXPO_FREE (expo);
 315   return mpfr_check_range (y, inex, rnd);
 316 }