gamma(a, x) returns incorrect values for a <= ~1/2

[jgwillingham]

I am using Julia v1.5.3 and SpecialFunctions v1.3.0.

The Problem

When evaluating the upper incomplete gamma function as gamma(1/2, x), the output is simply incorrect for certain values of x. Here is an example:

using SpecialFunctions
using Plots

X = range(0.0, stop=5.0, length=1000)
Y = gamma.(1/2, X)
plot(X,Y)

This produces the plot: That bump is clearly not right. If I zoom in on it, it seems like the problem begins at 3.0 and ends at pi. From the REPL:

julia> gamma(1/2, 3.0)
1.6882545542259388

julia> gamma(1/2, 3.0-1e-15)
0.02535650932346355

julia> gamma(1/2, pi)
0.0216042311666874

julia> gamma(1/2, pi-1e-14)
1.698504537077344

I have no idea why that interval might be special here. When I run gamma(a,x) for a < 1/2, I find similar problems which begin at x=3.0 but do not seem to end at x=pi. But I don't see any problem for a > 1/2.

EDIT: Upon a closer look, the problem above actually persists just above a=1/2 but it rapidly goes away. Increasing a by steps of 0.001, for length=100000, the interval of wrong values shrinks until it cannot be seen at a=0.508.

[stevengj]

Looks like the problem stems from a bug in the expint function:

julia> expint(0.5, 2.99) # correct
0.014831255784108155

julia> expint(0.5, 3.0) # incorrect
0.9747142213429575

cc @augustt198

[stevengj]

The crossover at x=3.0 comes from this line, which switches over from Taylor series to continued fractions at that point.