--- /dev/null
+/* mpfr_eint, mpfr_eint1 -- the exponential integral
+
+Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#define MPFR_NEED_LONGLONG_H
+#include "mpfr-impl.h"
+
+/* eint1(x) = -gamma - log(x) - sum((-1)^k*z^k/k/k!, k=1..infinity) for x > 0
+ = - eint(-x) for x < 0
+ where
+ eint (x) = gamma + log(x) + sum(z^k/k/k!, k=1..infinity) for x > 0
+ eint (x) is undefined for x < 0.
+*/
+
+/* compute in y an approximation of sum(x^k/k/k!, k=1..infinity),
+ and return e such that the absolute error is bound by 2^e ulp(y) */
+static mp_exp_t
+mpfr_eint_aux (mpfr_t y, mpfr_srcptr x)
+{
+ mpfr_t eps; /* dynamic (absolute) error bound on t */
+ mpfr_t erru, errs;
+ mpz_t m, s, t, u;
+ mp_exp_t e, sizeinbase;
+ mp_prec_t w = MPFR_PREC(y);
+ unsigned long k;
+ MPFR_GROUP_DECL (group);
+
+ /* for |x| <= 1, we have S := sum(x^k/k/k!, k=1..infinity) = x + R(x)
+ where |R(x)| <= (x/2)^2/(1-x/2) <= 2*(x/2)^2
+ thus |R(x)/x| <= |x|/2
+ thus if |x| <= 2^(-PREC(y)) we have |S - o(x)| <= ulp(y) */
+
+ if (MPFR_GET_EXP(x) <= - (mp_exp_t) w)
+ {
+ mpfr_set (y, x, GMP_RNDN);
+ return 0;
+ }
+
+ mpz_init (s); /* initializes to 0 */
+ mpz_init (t);
+ mpz_init (u);
+ mpz_init (m);
+ MPFR_GROUP_INIT_3 (group, 31, eps, erru, errs);
+ e = mpfr_get_z_exp (m, x); /* x = m * 2^e */
+ MPFR_ASSERTD (mpz_sizeinbase (m, 2) == MPFR_PREC (x));
+ if (MPFR_PREC (x) > w)
+ {
+ e += MPFR_PREC (x) - w;
+ mpz_tdiv_q_2exp (m, m, MPFR_PREC (x) - w);
+ }
+ /* remove trailing zeroes from m: this will speed up much cases where
+ x is a small integer divided by a power of 2 */
+ k = mpz_scan1 (m, 0);
+ mpz_tdiv_q_2exp (m, m, k);
+ e += k;
+ /* initialize t to 2^w */
+ mpz_set_ui (t, 1);
+ mpz_mul_2exp (t, t, w);
+ mpfr_set_ui (eps, 0, GMP_RNDN); /* eps[0] = 0 */
+ mpfr_set_ui (errs, 0, GMP_RNDN);
+ for (k = 1;; k++)
+ {
+ /* let eps[k] be the absolute error on t[k]:
+ since t[k] = trunc(t[k-1]*m*2^e/k), we have
+ eps[k+1] <= 1 + eps[k-1]*m*2^e/k + t[k-1]*m*2^(1-w)*2^e/k
+ = 1 + (eps[k-1] + t[k-1]*2^(1-w))*m*2^e/k
+ = 1 + (eps[k-1]*2^(w-1) + t[k-1])*2^(1-w)*m*2^e/k */
+ mpfr_mul_2ui (eps, eps, w - 1, GMP_RNDU);
+ mpfr_add_z (eps, eps, t, GMP_RNDU);
+ MPFR_MPZ_SIZEINBASE2 (sizeinbase, m);
+ mpfr_mul_2si (eps, eps, sizeinbase - (w - 1) + e, GMP_RNDU);
+ mpfr_div_ui (eps, eps, k, GMP_RNDU);
+ mpfr_add_ui (eps, eps, 1, GMP_RNDU);
+ mpz_mul (t, t, m);
+ if (e < 0)
+ mpz_tdiv_q_2exp (t, t, -e);
+ else
+ mpz_mul_2exp (t, t, e);
+ mpz_tdiv_q_ui (t, t, k);
+ mpz_tdiv_q_ui (u, t, k);
+ mpz_add (s, s, u);
+ /* the absolute error on u is <= 1 + eps[k]/k */
+ mpfr_div_ui (erru, eps, k, GMP_RNDU);
+ mpfr_add_ui (erru, erru, 1, GMP_RNDU);
+ /* and that on s is the sum of all errors on u */
+ mpfr_add (errs, errs, erru, GMP_RNDU);
+ /* we are done when t is smaller than errs */
+ if (mpz_sgn (t) == 0)
+ sizeinbase = 0;
+ else
+ MPFR_MPZ_SIZEINBASE2 (sizeinbase, t);
+ if (sizeinbase < MPFR_GET_EXP (errs))
+ break;
+ }
+ /* the truncation error is bounded by (|t|+eps)/k*(|x|/k + |x|^2/k^2 + ...)
+ <= (|t|+eps)/k*|x|/(k-|x|) */
+ mpz_abs (t, t);
+ mpfr_add_z (eps, eps, t, GMP_RNDU);
+ mpfr_div_ui (eps, eps, k, GMP_RNDU);
+ mpfr_abs (erru, x, GMP_RNDU); /* |x| */
+ mpfr_mul (eps, eps, erru, GMP_RNDU);
+ mpfr_ui_sub (erru, k, erru, GMP_RNDD);
+ if (MPFR_IS_NEG (erru))
+ {
+ /* the truncated series does not converge, return fail */
+ e = w;
+ }
+ else
+ {
+ mpfr_div (eps, eps, erru, GMP_RNDU);
+ mpfr_add (errs, errs, eps, GMP_RNDU);
+ mpfr_set_z (y, s, GMP_RNDN);
+ mpfr_div_2ui (y, y, w, GMP_RNDN);
+ /* errs was an absolute error bound on s. We must convert it to an error
+ in terms of ulp(y). Since ulp(y) = 2^(EXP(y)-PREC(y)), we must
+ divide the error by 2^(EXP(y)-PREC(y)), but since we divided also
+ y by 2^w = 2^PREC(y), we must simply divide by 2^EXP(y). */
+ e = MPFR_GET_EXP (errs) - MPFR_GET_EXP (y);
+ }
+ MPFR_GROUP_CLEAR (group);
+ mpz_clear (s);
+ mpz_clear (t);
+ mpz_clear (u);
+ mpz_clear (m);
+ return e;
+}
+
+/* Return in y an approximation of Ei(x) using the asymptotic expansion:
+ Ei(x) = exp(x)/x * (1 + 1/x + 2/x^2 + ... + k!/x^k + ...)
+ Assumes x >= PREC(y) * log(2).
+ Returns the error bound in terms of ulp(y).
+*/
+static mp_exp_t
+mpfr_eint_asympt (mpfr_ptr y, mpfr_srcptr x)
+{
+ mp_prec_t p = MPFR_PREC(y);
+ mpfr_t invx, t, err;
+ unsigned long k;
+ mp_exp_t err_exp;
+
+ mpfr_init2 (t, p);
+ mpfr_init2 (invx, p);
+ mpfr_init2 (err, 31); /* error in ulps on y */
+ mpfr_ui_div (invx, 1, x, GMP_RNDN); /* invx = 1/x*(1+u) with |u|<=2^(1-p) */
+ mpfr_set_ui (t, 1, GMP_RNDN); /* exact */
+ mpfr_set (y, t, GMP_RNDN);
+ mpfr_set_ui (err, 0, GMP_RNDN);
+ for (k = 1; MPFR_GET_EXP(t) + (mp_exp_t) p > MPFR_GET_EXP(y); k++)
+ {
+ mpfr_mul (t, t, invx, GMP_RNDN); /* 2 more roundings */
+ mpfr_mul_ui (t, t, k, GMP_RNDN); /* 1 more rounding: t = k!/x^k*(1+u)^e
+ with u=2^{-p} and |e| <= 3*k */
+ /* we use the fact that |(1+u)^n-1| <= 2*|n*u| for |n*u| <= 1, thus
+ the error on t is less than 6*k*2^{-p}*t <= 6*k*ulp(t) */
+ /* err is in terms of ulp(y): transform it in terms of ulp(t) */
+ mpfr_mul_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), GMP_RNDU);
+ mpfr_add_ui (err, err, 6 * k, GMP_RNDU);
+ /* transform back in terms of ulp(y) */
+ mpfr_div_2si (err, err, MPFR_GET_EXP(y) - MPFR_GET_EXP(t), GMP_RNDU);
+ mpfr_add (y, y, t, GMP_RNDN);
+ }
+ /* add the truncation error bounded by ulp(y): 1 ulp */
+ mpfr_mul (y, y, invx, GMP_RNDN); /* err <= 2*err + 3/2 */
+ mpfr_exp (t, x, GMP_RNDN); /* err(t) <= 1/2*ulp(t) */
+ mpfr_mul (y, y, t, GMP_RNDN); /* again: err <= 2*err + 3/2 */
+ mpfr_mul_2ui (err, err, 2, GMP_RNDU);
+ mpfr_add_ui (err, err, 8, GMP_RNDU);
+ err_exp = MPFR_GET_EXP(err);
+ mpfr_clear (t);
+ mpfr_clear (invx);
+ mpfr_clear (err);
+ return err_exp;
+}
+
+int
+mpfr_eint (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd)
+{
+ int inex;
+ mpfr_t tmp, ump;
+ mp_exp_t err, te;
+ mp_prec_t prec;
+ MPFR_SAVE_EXPO_DECL (expo);
+ MPFR_ZIV_DECL (loop);
+
+ MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd),
+ ("y[%#R]=%R inexact=%d", y, y, inex));
+
+ if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
+ {
+ /* exp(NaN) = exp(-Inf) = NaN */
+ if (MPFR_IS_NAN (x) || (MPFR_IS_INF (x) && MPFR_IS_NEG(x)))
+ {
+ MPFR_SET_NAN (y);
+ MPFR_RET_NAN;
+ }
+ /* eint(+inf) = +inf */
+ else if (MPFR_IS_INF (x))
+ {
+ MPFR_SET_INF(y);
+ MPFR_SET_POS(y);
+ MPFR_RET(0);
+ }
+ else /* eint(+/-0) = -Inf */
+ {
+ MPFR_SET_INF(y);
+ MPFR_SET_NEG(y);
+ MPFR_RET(0);
+ }
+ }
+
+ /* eint(x) = NaN for x < 0 */
+ if (MPFR_IS_NEG(x))
+ {
+ MPFR_SET_NAN (y);
+ MPFR_RET_NAN;
+ }
+
+ MPFR_SAVE_EXPO_MARK (expo);
+
+ /* Since eint(x) >= exp(x)/x, we have log2(eint(x)) >= (x-log(x))/log(2).
+ Let's compute k <= (x-log(x))/log(2) in a low precision. If k >= emax,
+ then log2(eint(x)) >= emax, and eint(x) >= 2^emax, i.e. it overflows. */
+ mpfr_init2 (tmp, 64);
+ mpfr_init2 (ump, 64);
+ mpfr_log (tmp, x, GMP_RNDU);
+ mpfr_sub (ump, x, tmp, GMP_RNDD);
+ mpfr_const_log2 (tmp, GMP_RNDU);
+ mpfr_div (ump, ump, tmp, GMP_RNDD);
+ /* FIXME: We really need mpfr_set_exp_t and mpfr_cmp_exp_t functions. */
+ MPFR_ASSERTN (MPFR_EMAX_MAX <= LONG_MAX);
+ if (mpfr_cmp_ui (ump, __gmpfr_emax) >= 0)
+ {
+ mpfr_clear (tmp);
+ mpfr_clear (ump);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_overflow (y, rnd, 1);
+ }
+
+ /* Init stuff */
+ prec = MPFR_PREC (y) + 2 * MPFR_INT_CEIL_LOG2 (MPFR_PREC (y)) + 6;
+
+ /* eint() has a root 0.37250741078136663446..., so if x is near,
+ already take more bits */
+ if (MPFR_GET_EXP(x) == -1) /* 1/4 <= x < 1/2 */
+ {
+ double d;
+ d = mpfr_get_d (x, GMP_RNDN) - 0.37250741078136663;
+ d = (d == 0.0) ? -53 : __gmpfr_ceil_log2 (d);
+ prec += -d;
+ }
+
+ mpfr_set_prec (tmp, prec);
+ mpfr_set_prec (ump, prec);
+
+ MPFR_ZIV_INIT (loop, prec); /* Initialize the ZivLoop controler */
+ for (;;) /* Infinite loop */
+ {
+ /* We need that the smallest value of k!/x^k is smaller than 2^(-p).
+ The minimum is obtained for x=k, and it is smaller than e*sqrt(x)/e^x
+ for x>=1. */
+ if (MPFR_GET_EXP (x) > 0 && mpfr_cmp_d (x, ((double) prec +
+ 0.5 * (double) MPFR_GET_EXP (x)) * LOG2 + 1.0) > 0)
+ err = mpfr_eint_asympt (tmp, x);
+ else
+ {
+ err = mpfr_eint_aux (tmp, x); /* error <= 2^err ulp(tmp) */
+ te = MPFR_GET_EXP(tmp);
+ mpfr_const_euler (ump, GMP_RNDN); /* 0.577 -> EXP(ump)=0 */
+ mpfr_add (tmp, tmp, ump, GMP_RNDN);
+ /* error <= 1/2 + 1/2*2^(EXP(ump)-EXP(tmp)) + 2^(te-EXP(tmp)+err)
+ <= 1/2 + 2^(MAX(EXP(ump), te+err+1) - EXP(tmp))
+ <= 2^(MAX(0, 1 + MAX(EXP(ump), te+err+1) - EXP(tmp))) */
+ err = MAX(1, te + err + 2) - MPFR_GET_EXP(tmp);
+ err = MAX(0, err);
+ te = MPFR_GET_EXP(tmp);
+ mpfr_log (ump, x, GMP_RNDN);
+ mpfr_add (tmp, tmp, ump, GMP_RNDN);
+ /* same formula as above, except now EXP(ump) is not 0 */
+ err += te + 1;
+ if (MPFR_LIKELY (!MPFR_IS_ZERO (ump)))
+ err = MAX (MPFR_GET_EXP (ump), err);
+ err = MAX(0, err - MPFR_GET_EXP (tmp));
+ }
+ if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - err, MPFR_PREC (y), rnd)))
+ break;
+ MPFR_ZIV_NEXT (loop, prec); /* Increase used precision */
+ mpfr_set_prec (tmp, prec);
+ mpfr_set_prec (ump, prec);
+ }
+ MPFR_ZIV_FREE (loop); /* Free the ZivLoop Controler */
+
+ inex = mpfr_set (y, tmp, rnd); /* Set y to the computed value */
+ mpfr_clear (tmp);
+ mpfr_clear (ump);
+
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (y, inex, rnd);
+}
projects / msp430-gcc.git / blobdiff