267 lines
8.1 KiB
C++
267 lines
8.1 KiB
C++
/* The functions exp1, expi below are based on translations of the Fortran code
|
|
* by Shanjie Zhang and Jianming Jin from the book
|
|
*
|
|
* Shanjie Zhang, Jianming Jin,
|
|
* Computation of Special Functions,
|
|
* Wiley, 1996,
|
|
* ISBN: 0-471-11963-6,
|
|
* LC: QA351.C45.
|
|
*/
|
|
|
|
#pragma once
|
|
|
|
#include "config.h"
|
|
#include "error.h"
|
|
|
|
#include "cephes/const.h"
|
|
|
|
|
|
namespace xsf {
|
|
|
|
|
|
XSF_HOST_DEVICE inline double exp1(double x) {
|
|
// ============================================
|
|
// Purpose: Compute exponential integral E1(x)
|
|
// Input : x --- Argument of E1(x)
|
|
// Output: E1 --- E1(x) ( x > 0 )
|
|
// ============================================
|
|
int k, m;
|
|
double e1, r, t, t0;
|
|
constexpr double ga = cephes::detail::SCIPY_EULER;
|
|
|
|
if (x == 0.0) {
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
if (x <= 1.0) {
|
|
e1 = 1.0;
|
|
r = 1.0;
|
|
for (k = 1; k < 26; k++) {
|
|
r = -r*k*x/std::pow(k+1.0, 2);
|
|
e1 += r;
|
|
if (std::abs(r) <= std::abs(e1)*1e-15) { break; }
|
|
}
|
|
return -ga - std::log(x) + x*e1;
|
|
}
|
|
m = 20 + (int)(80.0/x);
|
|
t0 = 0.0;
|
|
for (k = m; k > 0; k--) {
|
|
t0 = k / (1.0 + k / (x+t0));
|
|
}
|
|
t = 1.0 / (x + t0);
|
|
return std::exp(-x)*t;
|
|
}
|
|
|
|
XSF_HOST_DEVICE inline float exp1(float x) { return exp1(static_cast<double>(x)); }
|
|
|
|
XSF_HOST_DEVICE inline std::complex<double> exp1(std::complex<double> z) {
|
|
// ====================================================
|
|
// Purpose: Compute complex exponential integral E1(z)
|
|
// Input : z --- Argument of E1(z)
|
|
// Output: CE1 --- E1(z)
|
|
// ====================================================
|
|
constexpr double el = cephes::detail::SCIPY_EULER;
|
|
int k;
|
|
std::complex<double> ce1, cr, zc, zd, zdc;
|
|
double x = z.real();
|
|
double a0 = std::abs(z);
|
|
// Continued fraction converges slowly near negative real axis,
|
|
// so use power series in a wedge around it until radius 40.0
|
|
double xt = -2.0*std::abs(z.imag());
|
|
|
|
if (a0 == 0.0) { return std::numeric_limits<double>::infinity(); }
|
|
if ((a0 < 5.0) || ((x < xt) && (a0 < 40.0))) {
|
|
// Power series
|
|
ce1 = 1.0;
|
|
cr = 1.0;
|
|
for (k = 1; k < 501; k++) {
|
|
cr = -cr*z*static_cast<double>(k / std::pow(k + 1, 2));
|
|
ce1 += cr;
|
|
if (std::abs(cr) < std::abs(ce1)*1e-15) { break; }
|
|
}
|
|
if ((x <= 0.0) && (z.imag() == 0.0)) {
|
|
//Careful on the branch cut -- use the sign of the imaginary part
|
|
// to get the right sign on the factor if pi.
|
|
ce1 = -el - std::log(-z) + z*ce1 - std::copysign(M_PI, z.imag())*std::complex<double>(0.0, 1.0);
|
|
} else {
|
|
ce1 = -el - std::log(z) + z*ce1;
|
|
}
|
|
} else {
|
|
// Continued fraction https://dlmf.nist.gov/6.9
|
|
// 1 1 1 2 2 3 3
|
|
// E1 = exp(-z) * ----- ----- ----- ----- ----- ----- ----- ...
|
|
// Z + 1 + Z + 1 + Z + 1 + Z +
|
|
zc = 0.0;
|
|
zd = static_cast<double>(1) / z;
|
|
zdc = zd;
|
|
zc += zdc;
|
|
for (k = 1; k < 501; k++) {
|
|
zd = static_cast<double>(1) / (zd*static_cast<double>(k) + static_cast<double>(1));
|
|
zdc *= (zd - static_cast<double>(1));
|
|
zc += zdc;
|
|
|
|
zd = static_cast<double>(1) / (zd*static_cast<double>(k) + z);
|
|
zdc *= (z*zd - static_cast<double>(1));
|
|
zc += zdc;
|
|
if ((std::abs(zdc) <= std::abs(zc)*1e-15) && (k > 20)) { break; }
|
|
}
|
|
ce1 = std::exp(-z)*zc;
|
|
if ((x <= 0.0) && (z.imag() == 0.0)) {
|
|
ce1 -= M_PI*std::complex<double>(0.0, 1.0);
|
|
}
|
|
}
|
|
return ce1;
|
|
}
|
|
|
|
XSF_HOST_DEVICE inline std::complex<float> exp1(std::complex<float> z) {
|
|
return static_cast<std::complex<float>>(exp1(static_cast<std::complex<double>>(z)));
|
|
}
|
|
|
|
XSF_HOST_DEVICE inline double expi(double x) {
|
|
// ============================================
|
|
// Purpose: Compute exponential integral Ei(x)
|
|
// Input : x --- Argument of Ei(x)
|
|
// Output: EI --- Ei(x)
|
|
// ============================================
|
|
|
|
constexpr double ga = cephes::detail::SCIPY_EULER;
|
|
double ei, r;
|
|
|
|
if (x == 0.0) {
|
|
ei = -std::numeric_limits<double>::infinity();
|
|
} else if (x < 0) {
|
|
ei = -exp1(-x);
|
|
} else if (std::abs(x) <= 40.0) {
|
|
// Power series around x=0
|
|
ei = 1.0;
|
|
r = 1.0;
|
|
|
|
for (int k = 1; k <= 100; k++) {
|
|
r = r * k * x / ((k + 1.0) * (k + 1.0));
|
|
ei += r;
|
|
if (std::abs(r / ei) <= 1.0e-15) { break; }
|
|
}
|
|
ei = ga + std::log(x) + x * ei;
|
|
} else {
|
|
// Asymptotic expansion (the series is not convergent)
|
|
ei = 1.0;
|
|
r = 1.0;
|
|
for (int k = 1; k <= 20; k++) {
|
|
r = r * k / x;
|
|
ei += r;
|
|
}
|
|
ei = std::exp(x) / x * ei;
|
|
}
|
|
return ei;
|
|
}
|
|
|
|
XSF_HOST_DEVICE inline float expi(float x) { return expi(static_cast<double>(x)); }
|
|
|
|
std::complex<double> expi(std::complex<double> z) {
|
|
// ============================================
|
|
// Purpose: Compute exponential integral Ei(x)
|
|
// Input : x --- Complex argument of Ei(x)
|
|
// Output: EI --- Ei(x)
|
|
// ============================================
|
|
|
|
std::complex<double> cei;
|
|
cei = - exp1(-z);
|
|
if (z.imag() > 0.0) {
|
|
cei += std::complex<double>(0.0, M_PI);
|
|
} else if (z.imag() < 0.0 ) {
|
|
cei -= std::complex<double>(0.0, M_PI);
|
|
} else {
|
|
if (z.real() > 0.0) {
|
|
cei += std::complex<double>(0.0, copysign(M_PI, z.imag()));
|
|
}
|
|
}
|
|
return cei;
|
|
}
|
|
|
|
|
|
XSF_HOST_DEVICE inline std::complex<float> expi(std::complex<float> z) {
|
|
return static_cast<std::complex<float>>(expi(static_cast<std::complex<double>>(z)));
|
|
}
|
|
|
|
namespace detail {
|
|
|
|
//
|
|
// Compute a factor of the exponential integral E1.
|
|
// This is used in scaled_exp1(x) for moderate values of x.
|
|
//
|
|
// The function uses the continued fraction expansion given in equation 5.1.22
|
|
// of Abramowitz & Stegun, "Handbook of Mathematical Functions".
|
|
// For n=1, this is
|
|
//
|
|
// E1(x) = exp(-x)*C(x)
|
|
//
|
|
// where C(x), expressed in the notation used in A&S, is the continued fraction
|
|
//
|
|
// 1 1 1 2 2 3 3
|
|
// C(x) = --- --- --- --- --- --- --- ...
|
|
// x + 1 + x + 1 + x + 1 + x +
|
|
//
|
|
// Here, we pull a factor of 1/z out of C(x), so
|
|
//
|
|
// E1(x) = (exp(-x)/x)*F(x)
|
|
//
|
|
// and a bit of algebra gives the continued fraction expansion of F(x) to be
|
|
//
|
|
// 1 1 1 2 2 3 3
|
|
// F(x) = --- --- --- --- --- --- --- ...
|
|
// 1 + x + 1 + x + 1 + x + 1 +
|
|
//
|
|
XSF_HOST_DEVICE inline double expint1_factor_cont_frac(double x) {
|
|
// The number of terms to use in the truncated continued fraction
|
|
// depends on x. Larger values of x require fewer terms.
|
|
int m = 20 + (int) (80.0 / x);
|
|
double t0 = 0.0;
|
|
for (int k = m; k > 0; --k) {
|
|
t0 = k / (x + k / (1 + t0));
|
|
}
|
|
return 1 / (1 + t0);
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
//
|
|
// Scaled version of the exponential integral E_1(x).
|
|
//
|
|
// Factor E_1(x) as
|
|
//
|
|
// E_1(x) = exp(-x)/x * F(x)
|
|
//
|
|
// This function computes F(x).
|
|
//
|
|
// F(x) has the properties:
|
|
// * F(0) = 0
|
|
// * F is increasing on [0, inf)
|
|
// * lim_{x->inf} F(x) = 1.
|
|
//
|
|
XSF_HOST_DEVICE inline double scaled_exp1(double x) {
|
|
if (x < 0) {
|
|
return std::numeric_limits<double>::quiet_NaN();
|
|
}
|
|
|
|
if (x == 0) {
|
|
return 0.0;
|
|
}
|
|
|
|
if (x <= 1) {
|
|
// For small x, the naive implementation is sufficiently accurate.
|
|
return x * std::exp(x) * exp1(x);
|
|
}
|
|
|
|
if (x <= 1250) {
|
|
// For moderate x, use the continued fraction expansion.
|
|
return detail::expint1_factor_cont_frac(x);
|
|
}
|
|
|
|
// For large x, use the asymptotic expansion. This is equation 5.1.51
|
|
// from Abramowitz & Stegun, "Handbook of Mathematical Functions".
|
|
return 1 + (-1 + (2 + (-6 + (24 - 120 / x) / x) / x) / x) / x;
|
|
}
|
|
|
|
XSF_HOST_DEVICE inline float scaled_exp1(float x) { return scaled_exp1(static_cast<double>(x)); }
|
|
|
|
} // namespace xsf
|