search

JuliaMath / SpecialFunctions.jl

Special mathematical functions in Julia
https://specialfunctions.juliamath.org/stable/
Other
352 stars 99 forks source link

exponential integral (Ei, E₁, Eₙ...) function #19

Closed stevengj closed 3 years ago

stevengj commented 7 years ago

This is a common special function that it would be nice to include. It already is supplied by MPFR, which gives us a BigFloat version. (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.)

ararslan commented 7 years ago

You could always provide extra credit to students who submit PRs for that here 😉

elaineVRC commented 7 years ago

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 E1CFJ

Edit by ararslan: Code formatting for readability

stevengj commented 7 years ago

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
stevengj commented 7 years ago

(Unfortunately, the forwards recurrence for Eₙ used in expint(n::Integer, z) in my post above is inaccurate for large z due to cancellation errors, so it will need to switch to an explicit continued fraction for Eₙ.)

mschauer commented 7 years ago

For the time being, your code snippet above is under the MIT licence, I assume?

stevengj commented 7 years ago

Yes.

mschauer commented 6 years ago

I believe (take with grain) there is something not right in the way @E₁_cf64 z 15 is chosen for negative x with x² ≥ 7.84:

julia> z
-2.893143983832364 + 0.0im

julia> expint(z)
26031.00572721477 - 0.0im
stevengj commented 6 years ago

@mschauer, I'm not surprised, I mostly only looked at the right-half complex plane so far.

mschauer commented 6 years ago

Ah, I was maybe mislead by if (x² ≥ 546121) & (real(z) > 0)....

fabienlefloch commented 6 years ago

Another useful reference is the CERN library CEXPINT function, which is open source Kölbig, K. S. (1990). Exponential Integral for Complex Argument (No. CERNLIB-C338). CERN.

His rational and padé expansions are quite good.

pdenapo commented 5 years ago

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 Stacktrace: [1] top-level scope at none:0 [2] include at ./boot.jl:317 [inlined] [3] include_relative(::Module, ::String) at ./loading.jl:1041 [4] include(::Module, ::String) at ./sysimg.jl:29 [5] include(::String) at ./client.jl:388 [6] top-level scope at none:0 in expression starting at /home/pablo/Dropbox/pablo-ariel/exponencial_intergal.jl:8

mschauer commented 5 years ago

I do not know what you did, but in any case Expr(:macrocall, Symbol("@evalpoly"), :t, c...) will not work anymore because the IR for macro calls got a line number node and the :t is not hygienic (used in evalpoly)

mschauer commented 5 years ago

I fixed it on Bridge master https://github.com/mschauer/Bridge.jl/pull/37

feanor12 commented 5 years ago

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
end
using 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)
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     453.523 ns (0.00% GC)
  median time:      464.629 ns (0.00% GC)
  mean time:        505.681 ns (2.46% GC)
  maximum time:     124.984 μs (99.34% GC)
  --------------
  samples:          10000
  evals/sample:     197
using BenchmarkTools
@benchmark pyspecial.expi(-100)
BenchmarkTools.Trial: 
  memory estimate:  208 bytes
  allocs estimate:  6
  --------------
  minimum time:     3.191 μs (0.00% GC)
  median time:      3.373 μs (0.00% GC)
  mean time:        3.576 μs (0.00% GC)
  maximum time:     46.222 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     8
lnacquaroli commented 4 years ago

Hello, I was just wondering if this will be merged at some point. The implementation from @mschauer (https://github.com/JuliaMath/SpecialFunctions.jl/issues/19#issuecomment-435017198) seems to work for me, but I don't feel like installing Bridge.jl just for the expint function.

stevengj commented 3 years ago

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).