Skip to content

Commit 9af6925

Browse files
andreasnoackandreasnoack
committed
Fix gamma_inc for larga a and small x (and very large x)
1 parent 2e90d1b commit 9af6925

File tree

2 files changed

+76
-65
lines changed

2 files changed

+76
-65
lines changed

‎src/gamma_inc.jl

Lines changed: 36 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -782,7 +782,7 @@ See also [`gamma(z)`](@ref SpecialFunctions.gamma), [`gamma_inc_inv(a,p,q)`](@re
782782
"""
783783
gamma_inc(a::Real,x::Real,ind::Integer=0) = _gamma_inc(promote(float(a),float(x))...,ind)
784784

785-
function _gamma_inc(a::Float64,x::Float64,ind::Integer)
785+
function _gamma_inc(a::Float64, x::Float64, ind::Integer)
786786
iop = ind + 1
787787
acc = acc0[iop]
788788
if a<0.0 || x<0.0
@@ -801,12 +801,19 @@ function _gamma_inc(a::Float64,x::Float64,ind::Integer)
801801
if a >= big1[iop]
802802
l = x/a
803803
if l == 0.0
804-
return (0.0,1.0)
804+
return (0.0, 1.0)
805805
end
806806
s = 1.0 - l
807807
z = -LogExpFunctions.logmxp1(l)
808808
if z >= 700.0/a
809-
return (0.0,1.0)
809+
if abs(s) <= 2*eps(Float64)
810+
throw(DomainError((a, x, ind, "P(a,x) or Q(a,x) is computationally indeterminant in this case.")))
811+
end
812+
if x <= a
813+
return (0.0, 1.0)
814+
else
815+
return (1.0, 0.0)
816+
end
810817
end
811818
y = a*z
812819
rta = sqrt(a)
@@ -820,7 +827,7 @@ function _gamma_inc(a::Float64,x::Float64,ind::Integer)
820827

821828
if abs(s) <= 0.4
822829
if abs(s) <= 2.0*eps() && a*eps()*eps() > 3.28e-3
823-
throw(DomainError((a,x,ind,"P(a,x) or Q(a,x) is computationally indeterminant in this case.")))
830+
throw(DomainError((a, x, ind, "P(a,x) or Q(a,x) is computationally indeterminant in this case.")))
824831
end
825832
c = exp(-y)
826833
w = 0.5*erfcx(sqrt(y))
@@ -830,33 +837,20 @@ function _gamma_inc(a::Float64,x::Float64,ind::Integer)
830837
z=-z
831838
end
832839
if iop == 1
833-
return gamma_inc_minimax(a,x,z)
840+
return gamma_inc_minimax(a, x, z)
834841
elseif iop == 2
835-
return gamma_inc_temme(a,x,z)
842+
return gamma_inc_temme(a, x, z)
836843
else
837-
t = @horner(z , d00 , d0[1] , d0[2] , d0[3])
844+
t = @horner(z, d00, d0[1], d0[2], d0[3])
838845
return gamma_inc_temme_1(a, x, z, ind)
839846
end
840-
end
841-
elseif a > x || x >= x0[iop] || !isinteger(2*a)
842-
r = rgammax(a,x)
843-
if r == 0.0
844-
if x <= a
845-
return (0.0,1.0)
846-
else
847-
return (1.0,0.0)
848-
end
849-
end
850-
if x <= max(a,alog10)
851-
return gamma_inc_taylor(a, x, ind)
852-
elseif x < x0[iop]
853-
return gamma_inc_cf(a, x, ind)
854847
else
855-
return gamma_inc_asym(a, x, ind)
848+
return _gamma_inc_choose_algorithm(a, x, ind)
856849
end
850+
elseif a > x || x >= x0[iop] || !isinteger(2*a)
851+
return _gamma_inc_choose_algorithm(a, x, ind)
857852
else
858853
return gamma_inc_fsum(a,x)
859-
860854
end
861855
elseif a == 0.5
862856
if x >= 0.25
@@ -868,15 +862,12 @@ function _gamma_inc(a::Float64,x::Float64,ind::Integer)
868862
elseif x < 1.1
869863
return gamma_inc_taylor_x(a, x, ind)
870864
end
871-
r = rgammax(a,x)
865+
r = rgammax(a, x)
872866
if r == 0.0
873867
return (1.0, 0.0)
874868
else
875869
return gamma_inc_cf(a, x, ind)
876870
end
877-
878-
879-
880871
end
881872

882873
function _gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall
@@ -888,6 +879,24 @@ end
888879
_gamma_inc(a::Float32,x::Float32,ind::Integer) = Float32.(_gamma_inc(Float64(a),Float64(x),ind))
889880
_gamma_inc(a::Float16,x::Float16,ind::Integer) = Float16.(_gamma_inc(Float64(a),Float64(x),ind))
890881

882+
function _gamma_inc_choose_algorithm(a::Float64, x::Float64, ind::Int)
883+
r = rgammax(a, x)
884+
if r == 0.0
885+
if x <= a
886+
return (0.0, 1.0)
887+
else
888+
return (1.0, 0.0)
889+
end
890+
end
891+
if x <= max(a, alog10)
892+
return gamma_inc_taylor(a, x, ind)
893+
elseif x < x0[ind + 1]
894+
return gamma_inc_cf(a, x, ind)
895+
else
896+
return gamma_inc_asym(a, x, ind)
897+
end
898+
end
899+
891900
#EFFICIENT AND ACCURATE ALGORITHMS FOR THECOMPUTATION AND INVERSION OF THE INCOMPLETEGAMMA FUNCTION RATIOS by Amparo Gil, Javier Segura, Nico M. Temme
892901
#SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981
893902
# arXiv:1306.1754

‎test/gamma_inc.jl

Lines changed: 40 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1,43 +1,45 @@
11
@testset "incomplete gamma ratios" begin
22
#Computed using Wolframalpha gamma(a,x)/gamma(a) ~ gamma_q(a,x,0) function.
3-
@test gamma_inc(10,10)[2] ≈ 0.45792971447185221
4-
@test gamma_inc(1,1)[2] ≈ 0.3678794411714423216
5-
@test gamma_inc(0.5,0.5)[2] ≈ 0.31731050786291410
6-
@test gamma_inc(BigFloat(30.5),BigFloat(30.5))[2] ≈ parse(BigFloat,"0.47591691193354987004") rtol=eps()
7-
@test gamma_inc(5.5,0.5)[2] ≈ 0.9999496100513121669
8-
@test gamma_inc(0.5,7.4)[2] ≈ 0.0001195355018130302
9-
@test gamma_inc(0.5,0.22)[2] ≈ 0.507122455359825146
10-
@test gamma_inc(0.5,0.8)[2] ≈ 0.20590321073206830887
11-
@test gamma_inc(11.5,0.5)[2] ≈ 0.999999999998406112
12-
@test gamma_inc(0.19,0.99)[2] ≈ 0.050147247342905857
13-
@test gamma_inc(0.9999,0.9999)[2] ≈ 0.3678730556923103
14-
@test gamma_inc(24,23.9999999999)[2] ≈ 0.472849720555859138
15-
@test gamma_inc(0.5,0.55)[2] ≈ 0.29426610430496289
16-
@test gamma_inc(Float32(0.5),Float32(0.55))[2] ≈ Float32(gamma_inc(0.5,0.55)[2])
17-
@test gamma_inc(Float16(0.5),Float16(0.55))[2] ≈ Float16(gamma_inc(0.5,0.55)[2])
18-
@test gamma_inc(30,29.99999)[2] ≈ 0.475717712451705704
19-
@test gamma_inc(30,29.9)[2] ≈ 0.482992166284958565
20-
@test gamma_inc(10,0.0001)[2] ≈ 1.0000
21-
@test gamma_inc(0.0001,0.0001)[2] ≈ 0.000862958131006599
22-
@test gamma_inc(0.0001,10.5)[1] ≈ 0.999999999758896146
23-
@test gamma_inc(1,1)[1] ≈ 0.63212055882855768
24-
@test gamma_inc(13,15.1)[2] ≈ 0.25940814264863701
25-
@test gamma_inc(0.6,1.3)[2] ≈ 0.136458554006505355
26-
@test gamma_inc((100),(80))[2] ≈ 0.9828916869648668
27-
@test gamma_inc((100),(80),1)[2] ≈ 0.9828916869
28-
@test Float16(gamma_inc((100),(80),2)[2]) ≈ Float16(.983)
29-
@test gamma_inc(13.5,15.1)[2] ≈ 0.305242642543419087
30-
@test gamma_inc(11,9)[1] ≈ 0.2940116796594881834
31-
@test gamma_inc(8,32)[1] ≈ 0.99999989060651042057
32-
@test gamma_inc(15,16)[2] ≈ 0.3675273597655649298
33-
@test gamma_inc(15.5,16)[2] ≈ 0.4167440299455427811
34-
@test gamma_inc(0.9,0.8)[1] ≈ 0.59832030278768172
35-
@test gamma_inc(1.7,2.5)[1] ≈ 0.78446115627678957
36-
@test gamma_inc(11.1,0.001)[2] ≈ 1.0000
37-
@test gamma_inc(1e7, (1e7)+1)[1] ≈ 0.5001682088254367
38-
@test gamma_inc(1e7, (1e7)+1)[2] ≈ 0.4998317911745633
39-
@test_throws DomainError gamma_inc(-1,2,2)
40-
@test_throws DomainError gamma_inc(0,0,1)
3+
@test gamma_inc(10, 10)[2] ≈ 0.45792971447185221
4+
@test gamma_inc(1, 1)[2] ≈ 0.3678794411714423216
5+
@test gamma_inc(0.5, 0.5)[2] ≈ 0.31731050786291410
6+
@test gamma_inc(BigFloat(30.5), BigFloat(30.5))[2] ≈ parse(BigFloat,"0.47591691193354987004") rtol=eps()
7+
@test gamma_inc(5.5, 0.5)[2] ≈ 0.9999496100513121669
8+
@test gamma_inc(0.5, 7.4)[2] ≈ 0.0001195355018130302
9+
@test gamma_inc(0.5, 0.22)[2] ≈ 0.507122455359825146
10+
@test gamma_inc(0.5, 0.8)[2] ≈ 0.20590321073206830887
11+
@test gamma_inc(11.5, 0.5)[2] ≈ 0.999999999998406112
12+
@test gamma_inc(0.19, 0.99)[2] ≈ 0.050147247342905857
13+
@test gamma_inc(0.9999, 0.9999)[2] ≈ 0.3678730556923103
14+
@test gamma_inc(24, 23.9999999999)[2] ≈ 0.472849720555859138
15+
@test gamma_inc(0.5, 0.55)[2] ≈ 0.29426610430496289
16+
@test gamma_inc(Float32(0.5), Float32(0.55))[2] ≈ Float32(gamma_inc(0.5,0.55)[2])
17+
@test gamma_inc(Float16(0.5), Float16(0.55))[2] ≈ Float16(gamma_inc(0.5,0.55)[2])
18+
@test gamma_inc(30, 29.99999)[2] ≈ 0.475717712451705704
19+
@test gamma_inc(30, 29.9)[2] ≈ 0.482992166284958565
20+
@test gamma_inc(10, 0.0001)[2] ≈ 1.0000
21+
@test gamma_inc(0.0001, 0.0001)[2] ≈ 0.000862958131006599
22+
@test gamma_inc(0.0001, 10.5)[1] ≈ 0.999999999758896146
23+
@test gamma_inc(1, 1)[1] ≈ 0.63212055882855768
24+
@test gamma_inc(13, 15.1)[2] ≈ 0.25940814264863701
25+
@test gamma_inc(0.6, 1.3)[2] ≈ 0.136458554006505355
26+
@test gamma_inc(100, 80)[2] ≈ 0.9828916869648668
27+
@test gamma_inc(100, 80, 1)[2] ≈ 0.9828916869
28+
@test Float16(gamma_inc(100, 80, 2)[2]) ≈ Float16(.983)
29+
@test gamma_inc(13.5, 15.1)[2] ≈ 0.305242642543419087
30+
@test gamma_inc(11, 9)[1] ≈ 0.2940116796594881834
31+
@test gamma_inc(8, 32)[1] ≈ 0.99999989060651042057
32+
@test gamma_inc(15, 16)[2] ≈ 0.3675273597655649298
33+
@test gamma_inc(15.5, 16)[2] ≈ 0.4167440299455427811
34+
@test gamma_inc(0.9, 0.8)[1] ≈ 0.59832030278768172
35+
@test gamma_inc(1.7, 2.5)[1] ≈ 0.78446115627678957
36+
@test gamma_inc(11.1, 0.001)[2] ≈ 1.0000
37+
@test gamma_inc(1e7, 1e7 + 1)[1] ≈ 0.5001682088254367
38+
@test gamma_inc(1e7, 1e7 + 1)[2] ≈ 0.4998317911745633
39+
@test gamma_inc(29.0, 0.3)[1] ≈ 5.80834761514062e-47
40+
@test gamma_inc(29.0, 1000.0)[2] == 0.0
41+
@test_throws DomainError gamma_inc(-1, 2, 2)
42+
@test_throws DomainError gamma_inc(0, 0, 1)
4143
end
4244
@testset "inverse of incomplete gamma ratios" begin
4345
#Compared with Scipy.special.gammaincinv

0 commit comments

Comments
 (0)