Exponential integral
The exponential integral is a special function arising in mathematical analysis and applied mathematics, defined principally as E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt, where the integration path in the complex plane avoids the negative real axis and the origin to ensure the principal value.[1] This function is related to the incomplete gamma function via E_1(z) = \Gamma(0, z), and it exhibits a branch cut along the negative real axis, making it multi-valued for arguments on that ray.[1] Common notations include \Ei(x) for the real-variable case, defined as the Cauchy principal value \Ei(x) = -\int_{-x}^\infty \frac{e^{-t}}{t} \, dt = \int_{-\infty}^x \frac{e^t}{t} \, dt for x > 0, with the relation \Ei(-x) = -E_1(x).[1] Another variant is the complementary exponential integral \Ein(z) = \int_0^z \frac{1 - e^{-t}}{t} \, dt, an entire function connected to E_1(z) by E_1(z) = \Ein(z) - \ln z - \gamma, where \gamma \approx 0.57721 is the Euler-Mascheroni constant.[1] These functions satisfy recurrence relations, such as E_{n+1}(z) = \frac{1}{n} [e^{-z} - z E_n(z)] for the generalized form E_n(z) = \int_1^\infty \frac{e^{-zt}}{t^n} \, dt, and they appear in series expansions involving the Euler constant and logarithms.[1] The exponential integral plays a crucial role in diverse fields, including solutions to differential equations in quantum mechanics and thermodynamics, where it models decay processes and energy radiation from oscillators. In number theory, the related logarithmic integral \li(x) = \Ei(\ln x) provides asymptotic approximations for the prime-counting function \pi(x), with \li(x) - \pi(x) bounded by terms on the order of \sqrt{x} \ln x. Additional applications encompass quantum field theory, Gibbs phenomena in Fourier analysis, and Laplace equation solutions in semiconductor physics, underscoring its utility in handling integrals of the form \int \frac{e^{at}}{t} \, dt.[2]Definitions
Principal definitions
The exponential integral \operatorname{Ei}(x) is defined for real x > 0 as the Cauchy principal value \operatorname{Ei}(x) = \pvint_{-\infty}^{x} \frac{e^{t}}{t} \, \mathrm{d}t, where the principal value accounts for the singularity at t = 0; for x < 0, it is given by \operatorname{Ei}(x) = \int_{-\infty}^{x} \frac{e^{t}}{t} \, \mathrm{d}t, with the function exhibiting a discontinuity at x = 0.[1][2][3] The generalized exponential integral E_n(x) for positive integer n \geq 1 and \operatorname{Re}(x) > 0 is defined as E_n(x) = \int_1^\infty \frac{e^{-xt}}{t^n} \, \mathrm{d}t, with the special case E_1(x) serving as the primary form related to \operatorname{Ei}(x) via \operatorname{Ei}(x) = -E_1(-x) for x > 0.[1][2] For complex arguments z, the function E_1(z) is analytically continued with a branch cut along the negative real axis (-\infty, 0], where the integral path avoids crossing this cut or passing through the origin, ensuring the principal value is well-defined.[1][2] The notation for the exponential integral was introduced by J. W. L. Glaisher in 1870, with subsequent standardization appearing in comprehensive mathematical tables such as those by Abramowitz and Stegun in 1964.[2][3]Alternative integral representations
The exponential integral E_1(z) is represented by the line integral E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt, where the path of integration is straight from z to \infty in the direction of increasing real part, valid for |\arg z| < \pi. This form provides the principal branch and facilitates analytic continuation within the sector. The generalized exponential integral E_n(x) for integer n \ge 2 and x > 0 admits a double integral representation E_n(x) = \frac{e^{-x}}{(n-1)!} \int_0^\infty e^{-u} \frac{u^{n-1}}{x + u} \, du. This expression arises from interchanging the order of integration in the standard definition and the gamma function representation of powers, enabling efficient computation for large n or small x.[4]Properties
Series expansions
The exponential integral E_1(z) admits a convergent power series expansion that incorporates a logarithmic term to account for its branch point singularity at z = 0. Specifically, for z in the complex plane cut along the negative real axis (i.e., |\arg z| < \pi), E_1(z) = -\gamma - \ln z + \sum_{n=1}^\infty \frac{(-1)^{n+1} z^n}{n \cdot n!}, where \gamma is the Euler-Mascheroni constant.[4] This expansion is valid for all finite |z| in the cut plane, with the power series portion having an infinite radius of convergence.[4] The logarithmic term \ln z (principal branch) highlights the non-analytic behavior at the origin, where E_1(z) diverges as -\ln z for small |z|, while the series provides the regular (Taylor-expandable) part of the function around z = 0.[4] An alternative form expresses the series using the digamma function \psi(n+1): E_1(z) = -\ln z + e^{-z} \sum_{n=0}^\infty \frac{z^n}{n!} \psi(n+1), which also converges for all finite |z| in the cut plane and underscores the connection to the incomplete gamma function via \psi.[4] For the related exponential integral \operatorname{Ei}(x) with x > 0, the series is \operatorname{Ei}(x) = \gamma + \ln x + \sum_{n=1}^\infty \frac{x^n}{n \cdot n!}, converging for all finite x > 0.[4] This follows from the interrelation \operatorname{Ei}(x) = -E_1(-x) (principal value), where the analytic continuation of E_1(-x) across the branch cut yields the real-valued \operatorname{Ei}(x).[1] Near x = 0^+, \operatorname{Ei}(x) exhibits a logarithmic divergence similar to E_1(z), with the series capturing the smooth contributions.[4]Asymptotic expansions
The asymptotic expansion of the exponential integral E_1(z) for large |z| in the sector |\arg z| < 3\pi/2 is a divergent series derived from repeated integration by parts of the integral representation: E_1(z) \sim \frac{e^{-z}}{z} \sum_{m=0}^{\infty} \frac{(-1)^m m!}{z^m}, \quad |z| \to \infty, \quad |\arg z| \leq \frac{3\pi}{2} - \delta, where \delta > 0 is arbitrary but fixed. This expansion captures the leading behavior e^{-z}/z with successive corrections, but the factorial growth of the coefficients m! renders the series divergent for all finite z; optimal truncation at the term where the coefficients begin to increase minimizes the error, which is bounded by the first neglected term.[5] The sector of validity is determined by the analytic continuation of E_1(z), with Stokes lines occurring at \arg z = \pm \pi, where the dominant contribution switches and the subdominant exponential terms become significant. Across these lines, the asymptotic form acquires an additional exponentially small term, \pm 2\pi i e^{\mp i \pi} z^{-1}, ensuring smooth transition in overlapping sectors such as \pi/2 + \delta \leq \pm \arg z \leq 3\pi/2 - \delta. This Stokes phenomenon highlights the sectorial nature of the expansion, limiting its uniform applicability in the complex plane.[5] Beyond the standard truncation, asymptotics to all orders incorporate exponentially improved expansions, where the remainder after n terms is re-expanded using complementary error functions or generalized exponential integrals, yielding higher precision in transition regions near the Stokes lines. Additionally, the divergent series admits Borel summability: the Borel transform of \sum_{m=0}^{\infty} (-1)^m m! / z^{m+1} is \sum_{m=0}^{\infty} (t/z)^{m} / m! = e^{t/z}, and integrating \int_0^{\infty} e^{-s} e^{s/z} ds / z = E_1(z) e^{-z}/z recovers the exact function, providing a rigorous resummation method for the asymptotic series. These techniques reveal exponential corrections that are negligible in the primary sector but essential for global approximations.[5][6] For real positive arguments x > 0, the expansion implies simple bracketing inequalities that bound E_1(x) without truncation: \frac{e^{-x}}{x+1} < E_1(x) < \frac{e^{-x}}{x}. These follow from integral comparisons or continued fraction representations and sharpen for x > 1, where the relative error in the leading term e^{-x}/x is O(1/x).[7]Functional equations and identities
The exponential integral E_n(x) for positive integer n and x > 0 satisfies a fundamental recurrence relation obtained via integration by parts on its integral definition: E_{n+1}(x) = \frac{e^{-x} - x E_n(x)}{n}. This identity allows computation of higher-order exponential integrals from lower-order ones and holds for n \geq 1. A related functional equation arises from differentiating the defining integral, yielding \frac{d}{dx} E_{n+1}(x) = -E_n(x), or equivalently, x \frac{d}{dx} E_{n+1}(x) + n E_{n+1}(x) = e^{-x}, for n \geq 1. This linear first-order differential equation connects the derivative of E_{n+1}(x) to the function itself. For addition formulas, the exponential integral E_1(x) admits representations that relate E_1(x+y) to integrals involving shifted arguments, such as E_1(x+y) = e^{-y} \int_x^\infty \frac{e^{-u}}{u+y} \, du, \quad x > 0, \, y > 0. This form expresses the sum in the argument through an integral transform, without a simple closed-form addition theorem in terms of E_1(x) and E_1(y) alone. More general integral identities link E_1(x+y) to auxiliary functions, but they remain representational rather than algebraic. The exponential integral of imaginary argument, \operatorname{Ei}(x), exhibits symmetry and reflection properties across its branch cut along the negative real axis. For x > 0, \operatorname{Ei}(-x) = -E_1(x). This relates the principal value to the standard E_1. Additionally, the function satisfies a reflection formula reflecting the monodromy around the origin: \operatorname{Ei}(x e^{2m\pi i}) = \operatorname{Ei}(x) + 2m\pi i, \quad m \in \mathbb{Z}, for x > 0, capturing the discontinuity across the branch cut, where the jump is $2\pi i. These identities highlight the multivalued nature of \operatorname{Ei}(x) for complex arguments.Derivatives and integrals
The first derivative of the exponential integral E_1(z) is \frac{d}{dz} E_1(z) = -\frac{e^{-z}}{z}, valid for \operatorname{[Re](/page/Re)}(z) > 0, and follows from differentiating the integral representation E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt using the Leibniz rule for differentiating under the integral sign with a variable lower limit. Higher-order derivatives can be obtained recursively using the relation \frac{d}{dz} E_n(z) = -E_{n-1}(z), where E_n(z) denotes the generalized exponential integral for positive integer n, with E_0(z) = \frac{e^{-z}}{z}, yielding \frac{d^n}{dz^n} E_1(z) = (-1)^n \int_1^\infty t^{n-1} e^{-z t} \, dt for positive integer n and z > 0.[8] This integral form connects to the broader family of exponential integrals through recurrence relations, such as E_{n+1}(z) = \frac{1}{n} \left[ e^{-z} - z E_n(z) \right], which tie into functional identities for transforming arguments. The indefinite integral of E_1(z) is \int E_1(z) \, dz = z E_1(z) - e^{-z} + C, for \operatorname{Re}(z) > 0 and constant C, obtained by integration by parts with u = E_1(z) and dv = dz, using the first derivative formula. Differentiating this result recovers E_1(z), confirming its validity. Applications of the Leibniz rule arise in differentiating products involving exponential integrals, such as z E_1(z), where \frac{d}{dz} \left[ z E_1(z) \right] = E_1(z) + z \left( -\frac{e^{-z}}{z} \right) = E_1(z) - e^{-z}, linking directly to the recurrence E_2(z) = e^{-z} - z E_1(z). Similarly, for the product e^{z} E_n(z) with n \geq 1, \frac{d}{dz} \left[ e^{z} E_n(z) \right] = e^{z} E_{n-1}(z), providing a useful transformation for solving differential equations or evaluating related integrals. Higher-order Leibniz applications extend this to expressions like \frac{d^k}{dz^k} \left[ e^{z} E_n(z) \right], yielding sums involving lower-order E_m(z) terms.Related functions
Exponential integral of imaginary argument
The exponential integral of imaginary argument arises when the argument z is purely imaginary, z = i y with y real, and is particularly important in contexts involving oscillatory integrals, such as in quantum mechanics and signal processing. This variant exhibits complex values with real and imaginary parts directly tied to the sine and cosine integrals, highlighting its role in connecting exponential and trigonometric special functions. The function Ei(i y) is given by the relation \Ei(i y) = \Ci(y) + i \left( \Si(y) + \frac{\pi}{2} \right), where \Si(y) = \int_0^y \frac{\sin t}{t} , dt and \Ci(y) = -\int_y^\infty \frac{\cos t}{t} , dt for y > 0, with the cosine integral defined via its principal value. This expression follows from the analytic continuation of the exponential integral across the complex plane, substituting the imaginary argument into the integral representation. Likewise, the generalized exponential integral E_1(i y) satisfies E_1(i y) = -\Ci(y) + i \left( \Si(y) - \frac{\pi}{2} \right) for y > 0. These relations allow the imaginary-argument case to be computed using well-established algorithms for \Si(y) and \Ci(y), which are real-valued and exhibit bounded oscillatory behavior as y increases. For y < 0, the formulas hold with |y| and appropriate sign adjustments in the phase due to the even/odd properties of the underlying trigonometric functions. The principal branch of both Ei(z) and E_1(z) is defined with a branch cut along the negative real axis, ensuring the function is single-valued in the cut plane. Along the imaginary axis, there is no branch cut, and the functions are analytic; the principal value is inherently satisfied since paths from i y to infinity can avoid the cut entirely by staying in the right half-plane. This behavior ensures that Ei(i y) and E_1(i y) are continuous and well-defined for all real y ≠ 0, with a removable singularity at y = 0 approachable from either side. For large |y|, the asymptotic behavior of E_1(i y) is dominated by the leading term E_1(i y) \sim \frac{e^{-i y}}{i y}, followed by higher-order oscillatory corrections from the divergent asymptotic series \sum_{n=0}^\infty (-1)^n \frac{n! , e^{-i y}}{(i y)^{n+1}}. This expansion captures the decaying oscillatory nature, with amplitude scaling as 1/|y| and phase determined by the exponential, making it suitable for approximating integrals in high-frequency regimes. The series is asymptotic and diverges, but truncating before the minimal term provides accurate approximations for moderately large |y|.[9]Relations to other special functions
The exponential integral E_1(z) is directly related to the upper incomplete gamma function \Gamma(s, z) through the identity E_1(z) = \Gamma(0, z) for \Re(z) > 0, where \Gamma(s, z) = \int_z^\infty t^{s-1} e^{-t} \, dt. This connection arises because both functions share the integral form involving an exponential decay divided by the variable, and it allows the use of properties of the incomplete gamma for analyzing E_1(z).[10] More generally, the nth-order exponential integral E_n(z) is expressed as E_n(z) = z^{n-1} \Gamma(1 - n, z). The logarithmic integral \mathrm{li}(x), defined as the principal value integral \mathrm{pv} \int_0^x \frac{dt}{\ln t} for x > 1, is connected to the exponential integral by \mathrm{li}(x) = \mathrm{Ei}(\ln x), where \mathrm{Ei}(z) = -\int_{-z}^\infty \frac{e^{-t}}{t} \, dt is the principal branch. This relation is particularly useful in number theory, as \mathrm{li}(x) approximates the prime-counting function, and the exponential integral provides an analytic continuation.[11] The generalized exponential integral E_\nu(z) can be represented using the confluent hypergeometric function of the second kind U(a, b, z) as E_\nu(z) = z^{\nu-1} U(\nu, \nu, z) for \Re(\nu) > 0 and \Re(z) > 0.[12] This expression leverages the series and asymptotic properties of U to derive expansions for E_\nu(z), highlighting its role as a special case of confluent hypergeometric functions. For imaginary arguments, the exponential integral relates to the error function \mathrm{erf}(z) = \frac{2}{\sqrt{\pi}} \int_0^z e^{-t^2} \, dt through connections involving the imaginary error function \mathrm{erfi}(z) = -i \, \mathrm{erf}(i z). The Dawson function D(x) = e^{-x^2} \int_0^x e^{t^2} \, dt = \frac{\sqrt{\pi}}{2} e^{-x^2} \mathrm{erfi}(x) provides an example of this tie, with indirect connections to the exponential integral via complex parameter extensions and auxiliary functions like the Faddeeva function.[13] This extends to the sine integral \mathrm{Si}(y), which for imaginary arguments involves Gaussian integrals akin to those defining \mathrm{erf}.[14]Generalizations
Generalized exponential integrals
The generalized exponential integral extends the standard exponential integral E_1(z) to higher integer orders n > 1, defined for \Re z > 0 and positive integer n as E_n(z) = \int_1^\infty \frac{e^{-z t}}{t^n} \, dt. This form arises naturally in applications requiring repeated integration by parts of the base exponential integral, providing a unified framework for computing integrals involving powers in the denominator.[15] A key property facilitating computation and analysis is the reduction formula, which relates E_n(z) to the previous order: E_n(z) = \frac{1}{n-1} \left( e^{-z} - z E_{n-1}(z) \right), \quad n > 1. This recurrence allows recursive evaluation starting from E_1(z), with stability considerations for numerical implementation when n is large. In the more general form for real p > 1, p E_{p+1}(z) + z E_p(z) = e^{-z}, which rearranges to the analogous reduction E_{p+1}(z) = \frac{1}{p} \left( e^{-z} - z E_p(z) \right). These relations highlight the hierarchical structure of the family.[15] For non-integer orders, the exponential integral is generalized to E_\nu(z) where \nu \in \mathbb{R} with \Re \nu > 0 and \Re z > 0, defined via the integral representation E_\nu(z) = \int_1^\infty \frac{e^{-z t}}{t^\nu} \, dt = z^{\nu-1} \Gamma(1 - \nu, z), with \Gamma(s, z) denoting the upper incomplete gamma function \Gamma(s, z) = \int_z^\infty t^{s-1} e^{-t} \, dt. This extension leverages the analytic continuation of the incomplete gamma function, enabling the function to be defined for fractional \nu while preserving continuity with the integer case. The relation to the gamma function ensures that E_\nu(z) inherits meromorphic properties, with a branch point at z = 0 unless \nu is a nonpositive integer.[15] Further generalizations, such as those explored in the context of fractional calculus, incorporate the Wright function to handle non-integer parameters in differential equations involving exponential kernels. The Wright function {}_p \Psi_q \left( z \mid \begin{smallmatrix} (a_i, A_i)_{1,p} \\ (b_j, B_j)_{1,q} \end{smallmatrix} \right), introduced by E. M. Wright, provides a series representation that subsumes exponential integrals for specific parameter choices, particularly when modeling anomalous diffusion where fractional orders appear. This form is particularly useful for asymptotic analysis in non-integer regimes, though it requires careful parameter selection to recover the standard E_\nu(z).Exponential integrals with complex parameters
The exponential integral E_1(z) for complex argument z \neq 0 is defined via the principal branch by the contour integral E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, dt, where the path of integration starts at z and proceeds to infinity in a direction such that \operatorname{Re} t \to +\infty, without crossing the branch cut along the negative real axis or passing through the origin. This representation ensures analyticity in the complex plane cut along (-\infty, 0]. The function possesses a branch point at z = 0, rendering it multi-valued in the complex plane. Analytic continuation around the origin in the counterclockwise direction yields a monodromy relation E_1(z e^{2\pi i}) = E_1(z) - 2\pi i, arising from the logarithmic term in its Laurent series expansion near z = 0: E_1(z) = -\gamma - \ln z + \sum_{n=1}^\infty \frac{(-1)^{n+1} z^n}{n \cdot n!}, valid for |z| < 2\pi and |\arg(-z)| < \pi, where \gamma is the Euler-Mascheroni constant and the principal logarithm is used. The Riemann surface of E_1(z) is thus a logarithmic surface with infinitely many sheets, connected via the monodromy operator that shifts the function by multiples of -2\pi i. At infinity, E_1(z) exhibits an essential singularity, as its behavior cannot be captured by a Laurent series but requires an asymptotic expansion that diverges. The leading asymptotic form is E_1(z) \sim \frac{e^{-z}}{z} \sum_{k=0}^\infty \frac{(-1)^k k!}{z^k}, \quad |z| \to \infty, valid in the wide sector |\arg z| \leq \frac{3\pi}{2} - \delta for arbitrary fixed \delta > 0. Remainder estimates depend on the subsector: for |\arg z| \leq \frac{\pi}{2}, the error is bounded by the first omitted term, while in \frac{\pi}{2} < |\arg z| < \pi, it involves a factor of \csc(|\arg z|) times the first omitted term. To extend asymptotics across the full complex plane, sector decomposition is employed, dividing the plane into regions separated by Stokes lines (where \arg(-z) = (2m+1)\pi, m \in \mathbb{Z}) and matching expansions via connection formulas that account for subdominant exponential contributions. For the generalized exponential integral E_p(z) with complex order p \in \mathbb{C} and argument z \in \mathbb{C} \setminus \{0\}, the definition extends to E_p(z) = z^{p-1} \int_1^\infty e^{-z t} t^{-p} \, dt, assuming the path avoids the negative real axis; analytic continuation across the cut is achieved by deforming the integration path while respecting the branch of z^{p-1}. Unless p is a nonpositive integer, a branch point exists at z = 0, and each branch of E_p(z) is an entire function of p for fixed z \neq 0. The monodromy relation generalizes to E_p(z e^{2 m \pi i}) = \frac{2\pi i \, e^{m p \pi i}}{\Gamma(p)} \frac{\sin(m p \pi)}{\sin(p \pi)} z^{p-1} + E_p(z), for integer m, reflecting the interplay between the branching of z^{p-1} and the integral. Asymptotic expansions in complex sectors follow similar sectorial validity, with adjustments for the parameter p.Inverse exponential integral
Definition
The inverse exponential integral, denoted \mathrm{Ei}^{-1}(y), is defined as the function that satisfies the equation \mathrm{Ei}(x) = y, where x lies in the principal branch of the exponential integral \mathrm{Ei}(x). This inverse serves as the unique solution to the transcendental equation relating the input y to the output of the principal branch of \mathrm{Ei}. The principal branch of \mathrm{Ei}(x) for real x > 0 is strictly increasing and maps (0, \infty) bijectively onto (-\infty, \infty). Thus, the real inverse exists uniquely for all real y. In particular, \mathrm{Ei}(\ln \mu) = 0, where \mu \approx 1.45137 is the Ramanujan–Soldner constant.[1] Graphically, \mathrm{Ei}^{-1}(y) inverts the curve of \mathrm{Ei}(x), providing a tool to solve transcendental equations of the form \mathrm{Ei}(x) = y numerically or analytically where direct inversion is challenging.Properties and series
The inverse exponential integral \Ei^{-1}(y), defined as the unique positive solution to \Ei(x) = y for y > -\infty, is strictly increasing on its principal branch because the derivative \Ei'(x) = e^x / x > 0 for x > 0, ensuring monotonicity and uniqueness for each y in the range of \Ei.[1] Near y = 0, the function admits a Taylor series expansion \Ei^{-1}(y) = \sum_{n=0}^\infty \frac{y^n}{n!} P_n(\ln \mu) \mu^n, where \mu \approx 1.451 is the Ramanujan–Soldner constant, defined as the unique positive solution to li(μ) = 0 with li(x) = Ei(ln x) the logarithmic integral, and the polynomials P_n(z) are defined recursively by P_0(z) = z, P_{n+1}(z) = z (P_n'(z) - n P_n(z)) for |y| < \mu \ln \mu. This series provides a local representation around its zero x_0 = ln μ ≈ 0.37276, where \Ei(x_0) = 0. Higher-order terms in alternative expansions of the inverse can involve polylogarithms when inverting the series form of \Ei(x).[16] For large positive y, the asymptotic behavior is \Ei^{-1}(y) \sim \ln y + \ln \ln y + o(1), derived from the leading asymptotic of \Ei(x) \sim e^x / x for large x, which leads to the equation x e^{-x} = 1/y. Solving this yields x = -W_{-1}(-1/y), where W_{-1} is the -1 branch of the Lambert W function, whose asymptotic expansion for small negative arguments gives the stated form. This relation to the Lambert W function provides approximate solutions in limits where the exponential integral appears in transcendental equations, such as in certain physical models involving radiative transfer or population dynamics.Numerical evaluation
Approximation methods
For small values of z near the origin, where the exponential integral E_1(z) exhibits a logarithmic singularity, Padé approximants provide effective rational approximations by matching the Taylor series expansion of the regular part up to a specified order. These approximants are ratios of polynomials that converge more rapidly than the power series alone, particularly useful for hand calculations in the region |z| \lesssim 1. For instance, a [5/5] Padé approximant matches the series terms up to order 10, yielding relative errors on the order of $10^{-8} or better in the principal branch. A prominent continued fraction representation for E_1(z) facilitates approximations across a wide range of z, especially for moderate to large |z| away from the branch cut. Specifically, E_1(z) = \frac{e^{-z}}{z + \cfrac{1}{1 + \cfrac{1}{z + \cfrac{2}{1 + \cfrac{2}{z + \cfrac{3}{1 + \cdots}}}}} with convergence for |\ph z| < \pi. This form arises from integral representations and allows successive convergents to approximate E_1(z) by truncating the fraction, with rapid convergence for \Re z > 0. The partial denominators follow the pattern where even levels are 1 and odd levels are z + k for integer k. The convergents of this continued fraction yield inequality-based bounds for error estimation, as they alternate in sign around the true value, providing rigorous upper and lower enclosures. For example, the first convergent gives E_1(z) > e^{-z}/z, while the second provides E_1(z) < e^{-z} (z+1)/(z(z+1) + 1), with the error bounded by the difference between consecutive convergents. This property stems from the general theory of continued fractions for Stieltjes transforms, ensuring monotonic improvement and explicit error control without additional computation.[17] For large real arguments x > 2, simple rational approximations derived from the asymptotic expansion offer quick estimates suitable for initial hand calculations. One such form is E_1(x) \approx e^{-x} \left( \frac{1}{x} - \frac{1}{x^2} + \frac{2}{x^3} \right), which truncates the divergent series after three terms, achieving relative accuracy better than 1% for x \gtrsim 3 and improving as x increases. This rational expression in $1/x, multiplied by e^{-x}, captures the leading decay behavior while remaining computationally straightforward.Computational algorithms
For small values of the argument z, computational algorithms typically employ series summation of the Taylor expansion around zero, combined with argument reduction techniques such as scaling by powers of e to prevent overflow in intermediate computations. This approach ensures stability and convergence for |z| < 1, where the series terms decrease rapidly. For large |z|, asymptotic expansions are used, truncated at the optimal point to minimize the error, with the remainder estimated via the Euler-Maclaurin formula to achieve high accuracy without diverging terms. These methods are particularly effective in regions where |z| > 10, balancing computational cost and precision. Handling complex arguments requires careful branch selection based on the argument of z to respect the principal branch cut along the negative real axis, often partitioning the complex plane into sectors and applying power series expansions within each sector for uniform accuracy. Backward recurrence relations on three-term functional equations can further refine computations away from the branch cut.[18] Implementations in numerical libraries leverage these strategies for robust evaluation. SciPy'sexp1 function, based on the Cephes library, uses power series for small arguments and continued fraction approximations for larger ones, supporting complex inputs with relative errors below machine epsilon.[19][20] The GNU Scientific Library (GSL) provides gsl_sf_expint for real and integer-order cases, drawing from series and asymptotic methods in Abramowitz and Stegun.[21] For arbitrary precision, the Arb library computes exponential integrals using ball arithmetic to enclose results with certified error bounds, with updates through its 2023 merger into FLINT improving complex support and efficiency for high-precision regimes.[22]
Error analysis in these algorithms targets relative precision of approximately $10^{-15} or better for double-precision floating-point arithmetic, achieved through rigorous bounding of truncation and rounding errors in series and asymptotic terms, ensuring reliable results across the complex plane except near the branch cut.[23]