Skip to content

Fix the iop=2 Temme expansion variable and the ind=2 wide-band truncation #538

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
devmotion merged 5 commits into from
Aug 5, 2026
+80 −32
Merged

Fix the iop=2 Temme expansion variable and the ind=2 wide-band truncation #538

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e2faf49
Fix the iop=2 Temme expansion variable and the ind=2 wide-band trunca…
gaoflow Jul 31, 2026
2b880b8
refactor: share the Temme (P,Q) combine step; tighten wide/narrow-ban…
gaoflow Aug 3, 2026
1397469
Extend the Temme (P,Q) combine helper to all 5 call sites
gaoflow Aug 4, 2026
0c44e79
Bump version from 2.8.1 to 2.8.2
devmotion Aug 5, 2026
e148d24
Merge branch 'master' into fix-temme-iop2-and-iop3-wide
devmotion Aug 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "SpecialFunctions"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.8.1"
version = "2.8.2"

[deps]
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
Expand Down
51 changes: 20 additions & 31 deletions src/gamma_inc.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -542,6 +542,19 @@ function gamma_inc_taylor_x(a::Float64, x::Float64, ind::Integer)
end
end

# Combine the Temme series value `t` with the Gaussian-tail term `c*w` into (P, Q);
# shared by gamma_inc_minimax, gamma_inc_temme, gamma_inc_temme_1's three branches,
# and the ind=2 wide-band case in _gamma_inc.
function _gamma_inc_temme_combine(c::Float64, w::Float64, t::Float64, rta::Float64, l::Float64)
if l < 1.0
p = c*(w - rt2pin*t/rta)
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/rta)
return (1.0 - q, q)
end
end

@doc raw"""
gamma_inc_minimax(a,x,z)

Expand Down Expand Up @@ -580,13 +593,7 @@ function gamma_inc_minimax(a::Float64, x::Float64, z::Float64)
c6 = @evalpoly(z, d60, d6)

t = @evalpoly(1.0/a, c0, c1, c2, c3, c4, c5, c6, d70, d80)
if l < 1.0
p = c*(w - rt2pin*t/sqrt(a))
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/sqrt(a))
return (1.0 - q, q)
end
return _gamma_inc_temme_combine(c, w, t, sqrt(a), l)
end
#---USING THE MINIMAX APPROXIMATIONS---
c0 = @evalpoly(z, -.333333333333333E+00, -.159840143443990E+00, -.335378520024220E-01, -.231272501940775E-02)/(@evalpoly(z, 1.0, .729520430331981E+00, .238549219145773E+00, .376245718289389E-01, .239521354917408E-02, -.939001940478355E-05, .633763414209504E-06))
Expand All @@ -600,13 +607,7 @@ function gamma_inc_minimax(a::Float64, x::Float64, z::Float64)
c8 = @evalpoly(z, -.686013280418038E-03, .878371203603888E-03)

t = @evalpoly(1.0/a, c0, c1, c2, c3, c4, c5, c6, c7, c8)
if l < 1.0
p = c*(w - rt2pin*t/sqrt(a))
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/sqrt(a))
return (1.0 - q, q)
end
return _gamma_inc_temme_combine(c, w, t, sqrt(a), l)
end

@doc raw"""
Expand Down Expand Up @@ -635,13 +636,7 @@ function gamma_inc_temme(a::Float64, x::Float64, z::Float64)
c1 = @evalpoly(z, d10, d1[1], d1[2], d1[3], d1[4])
c2 = @evalpoly(z, d20, d2[1])
t = @evalpoly(1.0/a, c0, c1, c2)
if l < 1.0
p = c*(w - rt2pin*t/sqrt(a))
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/sqrt(a))
return (1.0 - q, q)
end
return _gamma_inc_temme_combine(c, w, t, sqrt(a), l)
end

@doc raw"""
Expand Down Expand Up @@ -681,21 +676,15 @@ function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer)
t = @evalpoly(u, c0, c1, c2, c3, c4, c5, c6, d70, d80)

elseif iop == 2
c0 = @evalpoly(d00, d0[1], d0[2])
c1 = @evalpoly(d10, d1[1])
c0 = @evalpoly(z, d00, d0[1], d0[2])
c1 = @evalpoly(z, d10, d1[1])
t = @evalpoly(u, c0, c1, d20)

else
t = @evalpoly(z, d00, d0[1])

end
if l < 1.0
p = c*(w - rt2pin*t/sqrt(a))
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/sqrt(a))
return (1.0 - q, q)
end
return _gamma_inc_temme_combine(c, w, t, sqrt(a), l)
end

"""
Expand Down Expand Up @@ -919,7 +908,7 @@ function _gamma_inc(a::Float64, x::Float64, ind::Integer)
return gamma_inc_temme(a, x, z)
else
t = @evalpoly(z, d00, d0[1], d0[2], d0[3])
return gamma_inc_temme_1(a, x, z, ind)
return _gamma_inc_temme_combine(c, w, t, rta, l)
end
else
return _gamma_inc_choose_algorithm(a, x, ind)
Expand Down
59 changes: 59 additions & 0 deletions test/gamma_inc.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,65 @@ end
end
end

# Temme expansions for ind=1 (iop==2) and the ind=2 wide band (iop==3): reference
# values from mpmath at 80+ digits (DLMF 8.7.1 series and 8.9.2 continued fraction).
@testset "incomplete gamma ratios: Temme iop 2/3 expansions" begin
@testset "ind=1 iop=2: a=$a, x=$x" for (a, x, p, q) in (
(14.0, 13.911135637064119, 0.5261046500921642, 0.4738953499078358),
(14.0, 14.088864362935881, 0.5449404484556186, 0.4550595515443813),
(100.0, 99.7625, 0.5038214442269839, 0.4961785557730161),
(100.0, 100.2375, 0.522753673458071, 0.477246326541929),
(1000.0, 999.24895905571, 0.4947284864885412, 0.5052715135114588),
(1000.0, 1000.7510409442899, 0.5136748881127599, 0.4863251118872402),
(10000.0, 9997.625, 0.4918547737457385, 0.5081452262542615),
(10000.0, 10002.374999999998, 0.5108025931978544, 0.48919740680214563),
(10000000.0, 9999924.895905571, 0.4905680282181972, 0.5094319717818028),
(10000000.0, 10000075.104094429, 0.5095160050557511, 0.4904839949442488),
(14.0, 13.906552106765321, 0.525615807496999, 0.474384192503001),
)
@test gamma_inc(a, x, 1)[1] ≈ p rtol=5e-7
@test gamma_inc(a, x, 1)[2] ≈ q rtol=5e-7
Comment on lines +103 to +104

@devmotion devmotion Aug 3, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these tolerances (and the ones below) as strict as possible, to increase the chances of detecting regressions?

end
@testset "ind=2 wide band: a=$a, x=$x" for (a, x, p, q) in (
(10.0, 6.5, 0.12261595071018043, 0.8773840492898196),
(10.0, 13.5, 0.8647360050511872, 0.1352639949488128),
(25.0, 16.25, 0.02613839150563149, 0.9738616084943685),
(25.0, 33.75, 0.9499593482838966, 0.05004065171610339),
(100.0, 65.0, 3.372775722939482e-05, 0.9999662722427706),
(100.0, 135.0, 0.999292148107214, 0.0007078518927860074),
(1000.0, 650.0, 2.957759012559303e-37, 1.0),
(1000.0, 1350.0, 1.0, 7.635408516228592e-24),
)
@test gamma_inc(a, x, 2)[1] ≈ p rtol=2e-4
@test gamma_inc(a, x, 2)[2] ≈ q rtol=2e-4
end
@testset "ind=2 narrow band: a=$a, x=$x" for (a, x, p, q) in (
(10.0, 9.778640563788214, 0.5140881193058611, 0.485911880694139),
(10.0, 10.221359436211785, 0.5694405131522309, 0.43055948684776907),
(100000.0, 99977.86405637882, 0.47251429019620883, 0.5274857098037912),
(100000.0, 100022.13594362116, 0.5283205849613148, 0.4716794150386851),
)
@test gamma_inc(a, x, 2)[1] ≈ p rtol=1e-4
@test gamma_inc(a, x, 2)[2] ≈ q rtol=1e-4
end
@testset "ind=0 unchanged: a=$a, x=$x" for (a, x, p, q) in (
(14.0, 13.911135637064119, 0.5261046500921642, 0.4738953499078358),
(14.0, 14.088864362935881, 0.5449404484556186, 0.4550595515443813),
(100.0, 99.7625, 0.5038214442269839, 0.4961785557730161),
(100.0, 100.2375, 0.522753673458071, 0.477246326541929),
(1000.0, 999.24895905571, 0.4947284864885412, 0.5052715135114588),
(1000.0, 1000.7510409442899, 0.5136748881127599, 0.4863251118872402),
(10000.0, 9997.625, 0.4918547737457385, 0.5081452262542615),
(10000.0, 10002.374999999998, 0.5108025931978544, 0.48919740680214563),
(10000000.0, 9999924.895905571, 0.4905680282181972, 0.5094319717818028),
(10000000.0, 10000075.104094429, 0.5095160050557511, 0.4904839949442488),
(14.0, 13.906552106765321, 0.525615807496999, 0.474384192503001),
)
@test gamma_inc(a, x)[1] ≈ p rtol=1e-12
@test gamma_inc(a, x)[2] ≈ q rtol=1e-12
end
end

@testset "inverse of incomplete gamma ratios" begin
#Compared with Scipy.special.gammaincinv
@test gamma_inc_inv(1.0,0.5,0.5) ≈ 0.69314718055994529
Expand Down
Loading