Skip to content

Commit d6cabd8

Browse files
andreasnoackandreasnoack
authored
Merge pull request #394 from JuliaMath/an/backport
More bugfixes to backport
2 parents 2ec1157 + f80715a commit d6cabd8

File tree

5 files changed

+120
-21
lines changed

5 files changed

+120
-21
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
1313
ChainRulesCore = "0.9.44, 0.10, 1"
1414
ChainRulesTestUtils = "0.6.8, 0.7, 1"
1515
IrrationalConstants = "0.1"
16-
LogExpFunctions = "0.3"
16+
LogExpFunctions = "0.3.2"
1717
OpenLibm_jll = "0.7, 0.8"
1818
OpenSpecFun_jll = "0.5"
1919
julia = "1.3"

‎src/beta_inc.jl

Lines changed: 13 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -371,19 +371,21 @@ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::F
371371
end
372372
z = -nu*lnx
373373
if b*z == 0.0
374-
return error("expansion can't be computed")
374+
@debug("underflow: expansion can't be computed")
375+
return w
375376
end
376377

377378
# COMPUTATION OF THE EXPANSION
378379
#SET R = EXP(-Z)*Z**B/GAMMA(B)
379380
r = b*(1.0 + rgamma1pm1(b))*exp(b*log(z))
380381
r *= exp(a*lnx)*exp(0.5*bm1*lnx)
381-
u = loggammadiv(b,a) + b*log(nu)
382+
u = loggammadiv(b, a) + b*log(nu)
382383
u = r*exp(-u)
383384
if u == 0.0
384-
return error("expansion can't be computed")
385+
@debug("underflow: expansion can't be computed")
386+
return w
385387
end
386-
(p, q) = gamma_inc(b,z,0)
388+
p, q = gamma_inc(b, z, 0)
387389
v = inv(nu)^2/4
388390
t2 = lnx^2/4
389391
l = w/u
@@ -412,7 +414,8 @@ function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::F
412414
dj = d[n] * j
413415
sm += dj
414416
if sm <= 0.0
415-
return error("expansion can't be computed")
417+
@debug("underflow: expansion can't be computed")
418+
return w
416419
end
417420
if abs(dj) <= epps*(sm+l)
418421
break
@@ -837,7 +840,7 @@ function _beta_inc(a::Float64, b::Float64, x::Float64, y::Float64=1-x)
837840
elseif a0 < min(epps, epps*b0) && b0*x0 <= 1.0
838841
q = beta_inc_power_series1(a0, b0, x0, epps)
839842
p = 1.0 - q
840-
elseif max(a0,b0) > 1.0
843+
elseif max(a0, b0) > 1.0
841844
if b0 <= 1.0
842845
p = beta_inc_power_series(a0, b0, x0, epps)
843846
q = 1.0 - p
@@ -1001,7 +1004,10 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
10011004
while true
10021005
p_approx = beta_inc(a, b, x)[1]
10031006
xin = x
1004-
p_approx = (p_approx - p)*min(floatmax(), exp(lb + r*log(xin) + t*log1p(-xin)))
1007+
p_approx = (p_approx - p)*min(
1008+
floatmax(),
1009+
exp(lb + LogExpFunctions.xlogy(r, xin) + LogExpFunctions.xlog1py(t, -xin))
1010+
)
10051011

10061012
if p_approx * p_approx_prev <= 0.0
10071013
prev = max(sq, fpu)

‎src/gamma_inc.jl

Lines changed: 90 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -226,7 +226,15 @@ Refer Eqn (3.12) in the paper
226226
"""
227227
function coeff1(eta::Float64)
228228
if abs(eta) < 1.0
229-
coeff1 = @horner(eta, -3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2, -4.923635739372e-3, -4.293658292782e-5) / @horner(eta, 1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1, 3.048648397436e-2, 1.605037988091e-3)
229+
coeff1 = @horner(
230+
eta,
231+
-3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2,
232+
-4.923635739372e-3, -4.293658292782e-5
233+
) / @horner(
234+
eta,
235+
1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1,
236+
3.048648397436e-2, 1.605037988091e-3
237+
)
230238
else
231239
la = lambdaeta(eta)
232240
coeff1 = log(eta/(la - 1.0))/eta
@@ -243,16 +251,39 @@ Refer Eqn (3.12) in the paper
243251
function coeff2(eta::Float64)
244252

245253
if eta < -5.0
246-
x=eta*eta
254+
x = eta*eta
247255
lnmeta = log(-eta)
248256
coeff2 = (12.0 - x - 6.0*lnmeta*lnmeta)/(12.0*x*eta)
249257
elseif eta < -2.0
250-
coeff2 = @horner(eta, -1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3, -6.06834887760e-4, -6.14830384279e-6) / @horner(eta, 1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1, 5.79490176079e-2, 5.74558524851e-3)
258+
coeff2 = @horner(
259+
eta,
260+
-1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3,
261+
-6.06834887760e-4, -6.14830384279e-6
262+
) / @horner(
263+
eta,
264+
1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1,
265+
5.79490176079e-2, 5.74558524851e-3)
251266
elseif eta < 2.0
252-
coeff2 = @horner(eta, -1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3, -3.91032032692e-4, 2.49634036069e-6) / @horner(eta, 1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1, 4.43843438769e-2, 4.24073217211e-3)
267+
coeff2 = @horner(
268+
eta,
269+
-1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3,
270+
-3.91032032692e-4, 2.49634036069e-6
271+
) / @horner(
272+
eta,
273+
1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1,
274+
4.43843438769e-2, 4.24073217211e-3
275+
)
253276
elseif eta < 1000.0
254277
x = 1.0/eta
255-
coeff2 = @horner(x, 9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2, 7.31901559577e+2, 4.55174411671e+1) / @horner(x, 1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2, 3.11993802124e+3, 3.97003311219e+3)
278+
coeff2 = @horner(
279+
x,
280+
9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2,
281+
7.31901559577e+2, 4.55174411671e+1
282+
) / @horner(
283+
x,
284+
1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2,
285+
3.11993802124e+3, 3.97003311219e+3
286+
)/(-12.0*eta)
256287
else
257288
coeff2 = -1.0/(12.0*eta)
258289
end
@@ -269,19 +300,65 @@ function coeff3(eta::Float64)
269300
if eta < -8.0
270301
x=eta*eta
271302
y = log(-eta)/eta
272-
coeff3=(-30.0+eta*y*(6.0*x*y*y-12.0+x))/(12.0*eta*x*x)
303+
coeff3=(-30.0 + eta*y*(6.0*x*y*y - 12.0 + x))/(12.0*eta*x*x)
273304
elseif eta < -4.0
274-
coeff3 = (@horner(eta, 4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3, 5.12634846317e-4, -2.01411722031e-5) / @horner(eta, 1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1, 4.64854522477e-2, 4.03751193496e-3))/(eta*eta)
305+
coeff3 = (
306+
@horner(
307+
eta,
308+
4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3,
309+
5.12634846317e-4, -2.01411722031e-5
310+
) / @horner(
311+
eta,
312+
1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1,
313+
4.64854522477e-2, 4.03751193496e-3
314+
)
315+
)/(eta*eta)
275316
elseif eta < -2.0
276-
coeff3 = @horner(eta, 4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5, -5.04476066942e-5, -5.35770949796e-6) / @horner(eta, 1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1, 9.01638932349e-2, 9.48935714996e-3)
317+
coeff3 = @horner(
318+
eta,
319+
4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5,
320+
-5.04476066942e-5, -5.35770949796e-6
321+
) / @horner(
322+
eta,
323+
1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1,
324+
9.01638932349e-2, 9.48935714996e-3
325+
)
277326
elseif eta < 2.0
278-
coeff3 = @horner(eta, 4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4, 5.29110969589e-6, 1.57166771750e-7) / @horner(eta, 1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1, 7.03527806143e-2, 8.06110846078e-3)
327+
coeff3 = @horner(
328+
eta,
329+
4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4,
330+
5.29110969589e-6, 1.57166771750e-7
331+
) / @horner(
332+
eta,
333+
1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1,
334+
7.03527806143e-2, 8.06110846078e-3
335+
)
279336
elseif eta < 10.0
280-
x=1.0/eta
281-
coeff3 = (@horner(x, -1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0, -2.18472031183e-1, 7.30002451555e-2) / @horner(x, 1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1, 2.18938950816e+2, 2.77067027185e+2))/(eta*eta)
337+
x = 1.0/eta
338+
coeff3 = (
339+
@horner(
340+
x,
341+
-1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0,
342+
-2.18472031183e-1, 7.30002451555e-2
343+
) / @horner(
344+
x,
345+
1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1,
346+
2.18938950816e+2, 2.77067027185e+2
347+
)
348+
)/(eta*eta)
282349
elseif eta < 100.0
283-
x=1.0/eta
284-
coeff3 = (@horner(x, -1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1, 1.99722374056e+2, -1.14311378756e+1) / @horner(x, 1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3, 7.11524019009e+3, 4.55746081453e+4))/(eta*eta)
350+
x = 1.0/eta
351+
coeff3 = (
352+
@horner(
353+
x,
354+
-1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1,
355+
1.99722374056e+2, -1.14311378756e+1
356+
) / @horner(
357+
x,
358+
1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3,
359+
7.11524019009e+3, 4.55746081453e+4
360+
)
361+
)/(eta*eta)
285362
else
286363
eta3 = eta*eta*eta
287364
coeff3 = -log(eta)/(12.0*eta3)

‎test/beta_inc.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -230,6 +230,11 @@
230230
@test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
231231
@test SpecialFunctions.loggammadiv(big(13.89), big(21.0001)) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
232232
@test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1)
233+
234+
@testset "Issue 334. Underflow without erroring" begin
235+
@test beta_inc(0.1, 4000, 0.2) == (1.0, 0.0)
236+
@test beta_inc(4000, 0.1, 0.2) == (0.0, 1.0)
237+
end
233238
end
234239

235240
@testset "inverse of incomplete beta" begin
@@ -289,4 +294,9 @@ end
289294
@test beta_inc_inv(a, b, p, 1 - p) === beta_inc_inv(a, b, p)
290295
end
291296
end
297+
298+
@testset "Avoid NaN when p=q=1" begin
299+
@test first(beta_inc_inv(1.0, 1.0, 1e-21)) ≈ 1e-21
300+
@test beta_inc_inv(1.0e30, 1.0, 0.49) == (1.0, 0.0)
301+
end
292302
end

‎test/gamma_inc.jl

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -201,6 +201,12 @@ end
201201
end
202202
end
203203
end
204+
205+
@testset "Issue 390 part 1" begin
206+
a = 1.0309015068677239
207+
q = 0.020202020202020204
208+
@test last(gamma_inc(a, gamma_inc_inv(a, 1 - q, q))) ≈ q
209+
end
204210
end
205211

206212
double(x::Real) = Float64(x)

0 commit comments

Comments
 (0)