exponential integral (Ei, E₁, Eₙ...) function #19
Comments
|
You could always provide extra credit to students who submit PRs for that here 😉 |
|
This is not a finished E1; this student has a lot to learn. My main concern was accuracy. I used series expansion and continued fraction type J from Cuyt (Handbook of Continued Fractions for Special Functions); also refer to AS (Abramowitz and Stegun) and NR (Numerical Recipes). I need to speed up the calculations. In [1]:
function eint(x::BigFloat)
z = BigFloat()
ccall((:mpfr_eint, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Int32), &z, &x, Base.MPFR.ROUNDING_MODE[end])
return z
end
eint(big(3.)),exp(-3.)*log(1.+1./3)#upper bound of E1(3.)
#lower bound is (exp(-3.)/2)*log(1 + 2./3)
#(9.93383257062541655800833601921676526299065302298758211969303
#4065242496228332754e+00 with 256 bits of precision,0.014322847009365602)
#this from mpfr is really Ei(x)
Out[1]:
(9.933832570625416558008336019216765262990653022987582119693034065242496228332754,0.014322847009365602)
In [ ]:
In [2]:
using PyCall
@pyimport scipy.special as s
s.exp1(3.0)
#0.013048381094197039
#reboot and call again !!
#PyCall not found
#at In[2]:1
###################################
#ipython notebook --profile julia
###################################
Out[2]:
0.013048381094197039
In [3]:
#-digamma(BigFloat("1"))#Euler-Mascheroni gamma
#5.77215664901532860606512090082402431042159335939923598805767
#2348848677267776685e-01 with 256 bits of precision
#WARNING: BigFloat(s::AbstractString) is deprecated,
# use parse(BigFloat,s) instead.
-digamma(parse(BigFloat,"1"))
Out[3]:
5.772156649015328606065120900824024310421593359399235988057672348848677267776685e-01
In [4]:
#series E1 accurate for xk in range > 0 to somewhere between (1. ,2.5)
function E1series(xk)
errfloat=1e-16;#make this smaller later
egamma= .5772156649015328606065120900;#change later
yk= -egamma -log(xk);#uses 5.1.22 AS
j=1;
pterm=xk;
term=xk;
while(abs(term) > errfloat)
yk=yk + term;
j=j+1;
pterm=-xk*pterm/j;
term=pterm/j;
end
return yk
end
#println(E1series(3))
#0.013048381094196789
0.013048381094197233
In [5]:
#E1CFs continued fraction Stype for x > somewhere in (1.,2.5)
function E1CFS(xk)
errfloat=1e-16;#make this smaller later
egamma= .5772156649015328606065120900;#change later
n = 1; #Eq 5.1.22 AS , 14.1.16 Cuyt S type continued fraction
#solving by Wallis method of forward recursions
# we're calculating E1(x)
# initialization
am2 = 0;
bm2 = 1;
am1 = 1;
bm1 = xk;
f = am1 / bm1;
oldf = Inf;
j = 2;
while (abs(f-oldf) > (100*errfloat.*abs(f)))
# calculate the coefficients of the recursion formulas for j even
alpha = n-1+(j/2); # note: beta= 1
#calculate A(j), B(j), and f(j)
a = am1 + alpha * am2;
b = bm1 + alpha * bm2;
# save new normalized variables for next pass through the loop
# note: normalization to avoid overflow or underflow
am2 = am1 / b;
bm2 = bm1 / b;
am1 = a / b;
bm1 = 1;
f = am1;
j = j+1;
# calculate the coefficients for j odd
alpha = (j-1)/2;
beta = xk;
a = beta * am1 + alpha * am2;
b = beta * bm1 + alpha * bm2;
am2 = am1 / b;
bm2 = bm1 / b;
am1 = a / b;
bm1 = 1;
oldf = f;
f = am1;
j = j+1;
end # end of the while loop
yk= exp(-xk) * f ;
return yk
end
# println( E1CFS(3))
#0.013048381094197106
0.013048381094197106
In [16]:
#J type continued fraction Cuyt eq. 14.1.23
#(NR 6.3.5)even part of S type converges faster
#using Lentz's algorithm to evaluate
#(sec 5.2 NR p 171)p 224 NR variation of recursion
function E1CFJ(x)
n=1
nm1=n-1
FPMIN=1e-30 #smallest number reset later
errfloat=1e-16
b=x+n
c=1./FPMIN
d=1./b
h=d
MAXIT=100
for i=1:MAXIT
a=-i*(nm1+i)
b=b+ 2.
d=1./(a*d + b)
c=b + a/c
del=c*d
h=h*del
if (abs(del -1.) < errfloat)
ans=h*exp(-x)
return ans
break
end# if
end#for
end#function
#println(E1CFJ(3.))
#0.013048381094197026
0.013048381094197026
In [78]:
using PyPlot
numx=1000;#was 100 or 200
#numx=200
x=linspace(2.45,2.5,numx);y=ones(numx);# 21.1 to 4.1 and 1.1 to 2.1
tic()
for xi=1:numx
xt=x[xi];
y[xi]=E1series(xt);
end
toc()
println ("time for $numx E1series")
println()
plot(x,y,color="b",marker=".",linestyle="None",label="E1series");
#plot(x, y, color="b", linewidth=2.0, linestyle="--",label="E1series")
#title("E1series")
hold(true)
#x=linspace(2.1,2.5,numx);y=ones(numx);
tic()
for xi=1:numx
xt=x[xi];
y[xi]=E1CFS(xt);
end
toc()
println("time for $numx E1CFS")
println()
plot(x, y, color="g",marker=".", linestyle="None",label="E1CFS");
hold(true)
#x=linspace(2.1,2.5,numx);y=ones(numx);
tic()
for xi=1:numx
xt=x[xi];
y[xi]=E1CFJ(xt);
end
toc()
println("time for $numx E1CFJ")
plot(x, y, color="m",marker=".", linestyle="None",label="E1CFJ");
elapsed time: 0.006186058 seconds
time for 1000 E1series
elapsed time:
Out[78]:
PyObject <matplotlib.legend.Legend object at 0x00000000411D09E8>
0.024557845 seconds
time for 1000 E1CFS
elapsed time: 0.007128253 seconds
time for 1000 E1CFJEdit by ararslan: Code formatting for readability |
|
For reference, here is my implementation. In a quick benchmark, it is 3–50 times faster than the above, presumably because of the inlining of all the polynomial and continued-fraction evaluations: ########################################################################
# Inlined, optimized code for the exponential integral E₁ in double precison.
# For more explanations, see course notes by Steven G. Johnson at
# https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
# n coefficients of the Taylor series of E₁(z) + log(z), in type T:
function E₁_taylor_coefficients{T<:Number}(::Type{T}, n::Integer)
n < 0 && throw(ArgumentError("$n ≥ 0 is required"))
n == 0 && return T[]
n == 1 && return T[-eulergamma]
# iteratively compute the terms in the series, starting with k=1
term::T = 1
terms = T[-eulergamma, term]
for k=2:n
term = -term * (k-1) / (k * k)
push!(terms, term)
end
return terms
end
# inline the Taylor expansion for a given order n, in double precision
macro E₁_taylor64(z, n::Integer)
c = E₁_taylor_coefficients(Float64, n)
taylor = Expr(:macrocall, Symbol("@evalpoly"), :t, c...)
quote
let t = $(esc(z))
$taylor - log(t)
end
end
end
# for numeric-literal coefficients: simplify to a ratio of two polynomials:
import Polynomials
# return (p,q): the polynomials p(x) / q(x) corresponding to E₁_cf(x, a...),
# but without the exp(-x) term
function E₁_cfpoly{T<:Real}(n::Integer, ::Type{T}=BigInt)
q = Polynomials.Poly(T[1])
p = x = Polynomials.Poly(T[0,1])
for i = n:-1:1
p, q = x*p+(1+i)*q, p # from cf = x + (1+i)/cf = x + (1+i)*q/p
p, q = p + i*q, p # from cf = 1 + i/cf = 1 + i*q/p
end
# do final 1/(x + inv(cf)) = 1/(x + q/p) = p/(x*p + q)
return p, x*p + q
end
macro E₁_cf64(x, n::Integer)
p,q = E₁_cfpoly(n, BigInt)
evalpoly = Symbol("@evalpoly")
num_expr = Expr(:macrocall, evalpoly, :t, Float64.(Polynomials.coeffs(p))...)
den_expr = Expr(:macrocall, evalpoly, :t, Float64.(Polynomials.coeffs(q))...)
quote
let t = $(esc(x))
exp(-t) * $num_expr / $den_expr
end
end
end
# exponential integral function E₁(z)
function expint(z::Union{Float64,Complex{Float64}})
x² = real(z)^2
y² = imag(z)^2
if x² + 0.233*y² ≥ 7.84 # use cf expansion, ≤ 30 terms
if (x² ≥ 546121) & (real(z) > 0) # underflow
return zero(z)
elseif x² + 0.401*y² ≥ 58.0 # ≤ 15 terms
if x² + 0.649*y² ≥ 540.0 # ≤ 8 terms
x² + y² ≥ 4e4 && return @E₁_cf64 z 4
return @E₁_cf64 z 8
end
return @E₁_cf64 z 15
end
return @E₁_cf64 z 30
else # use Taylor expansion, ≤ 37 terms
r² = x² + y²
return r² ≤ 0.36 ? (r² ≤ 2.8e-3 ? (r² ≤ 2e-7 ? @E₁_taylor64(z,4) :
@E₁_taylor64(z,8)) :
@E₁_taylor64(z,15)) :
@E₁_taylor64(z,37)
end
end
expint{T<:Integer}(z::Union{T,Complex{T},Rational{T},Complex{Rational{T}}}) = expint(float(z))
######################################################################
# exponential integral Eₙ(z)
function expint(n::Integer, z)
if n == 1
return expint(z)
elseif n < 1
# backwards recurrence from E₀ = e⁻ᶻ/z
zinv = inv(z)
e⁻ᶻ = exp(-z)
Eᵢ = zinv * e⁻ᶻ
for i = 1:-n
Eᵢ = zinv * (e⁻ᶻ + i * Eᵢ)
end
return Eᵢ
elseif n > 1
# forwards recurrence from E₁
e⁻ᶻ = exp(-z)
Eᵢ = expint(z)
for i = 2:n
Eᵢ = (e⁻ᶻ - z*Eᵢ) / (i - 1)
end
return Eᵢ
end
end |
|
(Unfortunately, the forwards recurrence for Eₙ used in |
|
For the time being, your code snippet above is under the MIT licence, I assume? |
|
Yes. |
|
I believe (take with grain) there is something not right in the way |
|
@mschauer, I'm not surprised, I mostly only looked at the right-half complex plane so far. |
|
Ah, I was maybe mislead by |
|
Another useful reference is the CERN library CEXPINT function, which is open source His rational and padé expansions are quite good. |
|
I have tried the code above by Steven Johnson on Julia 1.0.1 and got the following error message. May be something has changed in the language since then? Or some package is missing? ERROR: LoadError: UndefVarError: T not defined |
|
I do not know what you did, but in any case |
|
I fixed it on Bridge master |
|
I implemented the function EIX from https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f in julia. I'm not sure if the copyright is a problem. Also I think it is not type stable yet. using PyCall
pyspecial=pyimport("scipy.special")function _e1xb(n)
if n==0
return Inf
elseif n<=0
# power series around x=0
E1=1.0
R=1.0
for k in 1:25
R=-R*k*n/(k+1.0)^2
E1+=R
if abs(R)<=abs(E1)*1e-15
break
end
end
GA = 0.5772156649015328
return -GA-log(n)+n*E1
else
M=20+div(80,n)
T0=0.0
for k in M:-1:1
T0 = k/(1.0+k/(n+T0))
end
T0=1.0/(n+T0)
return exp(-n)*T0
end
end
function _expi(n::AbstractFloat)
if n==0
return -Inf
elseif n < 0
return -_e1xb(-n)
elseif n <= 40
#power series around x=0
E1=1.0
R=1.0
for k in 1:100
R=R*k*n/(k+1.0)^2
E1+=R
if abs(R/E1)<=1e-15
break
end
end
GA = 0.5772156649015328
return GA+log(n)+n*E1
else
# Asymptotic expansion (the series is not convergent)
E1=1.0
R=1.0
for k in 1:20
R*=k/n
E1+=R
end
return exp(n)/n*E1
end
endusing Test
@test _expi(-12.0)≈pyspecial.expi(-12.0)
@test _expi(12.0)≈pyspecial.expi(12.0)
@test _expi(0.0)≈pyspecial.expi(0.0)
@test _expi(-100.0)≈pyspecial.expi(-100.0)
@test _expi(100.0)≈pyspecial.expi(100.0)Test Passed using BenchmarkTools
@benchmark _expi(-100.0)using BenchmarkTools
@benchmark pyspecial.expi(-100) |
|
Hello, I was just wondering if this will be merged at some point. The implementation from @mschauer (#19 (comment)) seems to work for me, but I don't feel like installing Bridge.jl just for the |
* Added expint from Bridge.jl and JuliaMath/SpecialFunctions.jl#19 * Renamed some files * Renamed some functions * Updated the style of the code * Added anisotropy options to RI functions * Added new RI functions * Added multioscillators options to RI functions
|
Closed by #236 — the implementation by @augustt198 seems to be the fastest available code compared to those we've benchmarked, and it is also the most general (supporting arbitrary complex arguments and complex order). |
This is a common special function that it would be nice to include. It already is supplied by MPFR, which gives us a
BigFloatversion. (Issue copied/moved from JuliaLang/julia#7089.)See also the julia-users discussion on exponential integrals.
Some potentially useful references:
Of course, the usual copyright caveats apply: don't even look at the source code of any of these, just the equations. Especially for the ACM TOMS journal, the most evil journal in numerical analysis.
I assigned the problem of implementing
E₁to double-precision accuracy in the right-half complex plain (real(z) ≥ 0) in problem set 3 of our course 18.S096 this January. In the solutions, I showed how to implement it using a combination of Taylor series and continued fractions, with some custom macros to do inlining, and got performance 5–6 times faster than the Fortran code used in SciPy.This should be helpful as a starting point for more full-featured implementations. (The macros to simplify and inline continued-fraction expansions may also be useful in other special functions.)
The text was updated successfully, but these errors were encountered: