Exponential integral function

[augustt198]

This PR adds the exponential integral En(ν, z) for both complex ν and z. The special case ν = 1 is optimized in function E₁(z), which is based on Steven Johnson's code in #19.

Not all test cases are passing yet, so this is still a draft.

TODO:

• E₁(z)
  • Remove Polynomials.jl dependency
  • Tune Taylor series crossover
  • Use SIMD optimized evalpoly
• En(ν, z)
  • Consistently handle large ν: avoid multiplying complex numbers with Inf
  • Find better crossover boundary between Taylor series & continued fractions
  • Add more test cases with integer arguments: worried about hitting poles in gamma function
  • Better criteria for choosing to use the continued fraction with gamma: En(41 - 23im, 30 + 404im) fails
  • Type stability
• Documentation

[stevengj]

(@augustt198 has been working on this with me as part of UROP project at MIT.)

[stevengj]

You did a bunch of benchmarks against external codes — it would be good to post a few of those here.

[stevengj]

Can we now provide complex support for gamma_inc? How does this code compare to gamma_inc in the cases where they agree?

[stevengj]

Maybe just use ratfn_minimax from https://github.com/simonbyrne/Remez.jl/blob/master/src/Remez.jl?

[stevengj]

ApproxFun claims that a degree-16 polynomial gives machine precision over the 2.15..3.0 interval, for example, and presumably a minimax function could do even better.

[stevengj]

Note that you will want to pass a weighting function of 1/E₁ to ratfn_minimax so that it minimizes the maximum relative error.

[stevengj]

Tests are currently failing, for example:

  Expression: expint(1, -5 + 0.01im) ⩭ -40.18408805785205 - 2.8447696994032254im
   Evaluated: -40.18408805608827 - 2.8447696977664467im ⩭ -40.18408805785205 - 2.8447696994032254im

for example

Test Summary:          | Pass  Fail  Error  Total
expint                 |  925     1      5    931
  exponential integral |  925     1      5    931
    E1 real            |  309                 309
    E1 complex         |  488                 488
    En                 |  128     1      5    134

[stevengj]

test failure:

E1 complex: Test Failed at /Users/travis/build/JuliaMath/SpecialFunctions.jl/test/expint.jl:31
  Expression: expint(x) ≅ y
   Evaluated: 5.733317091369105e-42 + 4.0193199136184414e-39im ≅ 5.7333170913660523e-42 + 4.019319913618439e-39im

The relative error in the real part is 5.3e-13, which is slightly above your tolerance but should be okay — apparently due to a slight compiler difference on another platform?

[stevengj]

real ν > 50 are supported, so we should test some of those cases.

[stevengj]

Need to add expint to https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/docs/src/functions_list.md
 and https://github.com/JuliaMath/SpecialFunctions.jl/blob/master/docs/src/functions_overview.md

[stevengj]

Regarding the failing test, since it seems to just be slightly worse roundoff error in the real part on one machine (and the overall error in z is quite small), I don't think it's worth getting to the bottom of what is probably a compiler quirk — just test it with the lower 5e-12 tolerance from the other tests.

[stevengj]

Note: you can implement the incomplete gamma function in terms of this. To get the equivalent of gamma_inc, do:

gamma_inc2(a,x) = let q = expint(1-a,x) * x^a / gamma(a); 1-q, q; end

(However, this might give an inaccurate first value, the lower incomplete gamma function ratio, if q is close to 1.)

It actually seems to be significantly faster than our current gamma_inc function, however, and has the additional advantage of supporting complex arguments!

[eschnett]

Do you want to make a new minor release for expint?

[stevengj]

Once #263 is fixed, which should happen next week.

[stevengj]

Tagged. Since the api has been stable for a while, and Julia’s package system strongly favors >= 1.0 versions for compatibility tracking, I tagged it as 1.0

[oscardssmith]

Should there be documentation for this?

[stevengj]

It's in the "dev" docs; I'm not sure why it hasn't appeared in the "stable" docs yet?

[giordano]

For some reasons I don't currently understand, documentation is built only on pushes to master, but not creation of tags. At a quick glance, the setup looks correct to me

[giordano]

Oh, I see now: the setup was fixed with #287, which was merged on December 21st, but the latest release is v1.2.1, from December 16th. So the next release should hopefully have the documentation built regularly. See also #286

[stevengj]

Should be fixed by JuliaRegistries/General#30169

[giordano]

Yes, https://juliamath.github.io/SpecialFunctions.jl/stable/ now points to v1.3.0