Skip to content

Commit 48f22a7

Browse files
authored
fix spurious expint domain errors (#273)
* fix spurious expint domain errors * uncomment test
1 parent eff906e commit 48f22a7

File tree

2 files changed

+8
-1
lines changed

2 files changed

+8
-1
lines changed

src/expint.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,7 @@ end
197197

198198
# Calculate Γ(1 - ν) * z^(ν-1) safely
199199
En_safe_gamma_term::Number, z::Number) = exp((ν - 1)*log(z) + loggamma(1 - oftype(z, ν)))
200+
En_safe_gamma_term::Integer, z::Real) = (z 0 || isodd(ν) ? 1 : -1) * exp((ν - 1)*log(abs(z)) + loggamma(1 - oftype(z, ν)))
200201

201202
# continued fraction for En(ν, z) that uses the gamma function:
202203
# https://functions.wolfram.com/GammaBetaErf/ExpIntegralE/10/0005/
@@ -369,7 +370,7 @@ function expint(ν::Number, z::Number, niter::Int=1000)
369370
return oftype(z, NaN) * z
370371
end
371372

372-
if z isa Real && (isinteger(ν) ? (z < 0 && ν > 0) : (z < 0 || ν < 0))
373+
if z isa Real && (isinteger(ν) ? (z < 0 && ν > 0) : z < 0)
373374
throw(DomainError(z, "En will only return a complex result if called with a complex argument"))
374375
end
375376

test/expint.jl

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,6 +121,12 @@ using Base.MathConstants
121121
@test expint(-4, big"5.7") big"0.0013051784816041557674219259305064008704427556690521012788554962045721241678128274" rtol=1e-76
122122
@test expint(+4, big"5.7") big"0.00035851619983733810063957822090566023258844652925986297947736355731994111784745604432" rtol=1e-74
123123
end
124+
125+
# issue #272: more cases with real results
126+
@test expint(-2, 2.2) 0.11696351427428933590118175681379776 rtol=1e-14
127+
@test expint(-2, -2.2) -2.0680909972407264331884989 rtol=1e-14
128+
@test expint(-2.2, 3.2) 0.024950173497409329191241353395358 rtol=1e-14
129+
#@test expint(+2.2, 3.2) ≈ 0.008044603700773423319087602010 rtol=1e-14
124130
end
125131

126132
expinti_real(x) = invoke(expinti, Tuple{Real}, x)

0 commit comments

Comments
 (0)