Skip to content

Expand the minimax Temme series in 1/a rather than z #534

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 2 commits into from
Jul 31, 2026
Merged

Expand the minimax Temme series in 1/a rather than z #534

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b63e681
Expand the minimax Temme series in 1/a, not z
gaoflow Jul 27, 2026
8a25843
Bump version to 2.8.1
gaoflow Jul 31, 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 src/gamma_inc.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -579,7 +579,7 @@ function gamma_inc_minimax(a::Float64, x::Float64, z::Float64)
c5 = @evalpoly(z, d50, d5)
c6 = @evalpoly(z, d60, d6)

t = @evalpoly(z, c0, c1, c2, c3, c4, c5, c6, d70, d80)
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)
Expand Down
27 changes: 27 additions & 0 deletions test/gamma_inc.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,33 @@
@test_throws DomainError gamma_inc(1.0, -Inf)
end

# Reference values from mpmath at 50+ digits, cross-checked against the DLMF 8.7.1
# series and the 8.9.2 continued fraction.
@testset "incomplete gamma ratios: Temme minimax series" begin
# a >= 25 and e0/sqrt(a) < |1 - x/a| <= 1e-3, i.e. the gamma_inc_minimax series branch.
@testset "a=$a, x=$x" for (a, x, p, q) in (
(25.0, 24.97525, 0.5246323715433058506793, 0.4753676284566941493207),
(25.0, 25.02475, 0.5285687428518152079350, 0.4714312571481847920650),
(25.0, 24.9925, 0.5260050200501651608654, 0.4739949799498348391346),
(50.0, 50.0495, 0.5215950013484762091395, 0.4784049986515237908605),
(50.0, 49.975, 0.5173998411006592979910, 0.4826001588993407020090),
(134.0, 134.05, 0.5132100266748297555471, 0.4867899733251702444529),
(134.0, 133.86734, 0.5069170257485132637813, 0.4930829742514867362187),
(400.0, 400.396, 0.5145421195364036819570, 0.4854578804635963180430),
(400.0, 399.92, 0.5050535378251085332799, 0.4949464621748914667201),
(2000.0, 2001.98, 0.5206211429087569237891, 0.4793788570912430762109),
(2000.0, 1998.5, 0.4895906656999191847208, 0.5104093343000808152792),
(10000.0, 9990.1, 0.4618797888026423050784, 0.5381202111973576949216),
(10000.0, 10005.0, 0.5212634681357695888112, 0.4787365318642304111888),
(100000.0, 100099.0, 0.6232456486261727122096, 0.3767543513738272877904),
(100000.0, 99901.0, 0.3774766852491783530590, 0.6225233147508216469410),
(100000.0, 99985.0, 0.4815027198019921823376, 0.5184972801980078176624),
)
@test gamma_inc(a, x)[1] ≈ p rtol=1e-13
@test gamma_inc(a, x)[2] ≈ q rtol=1e-13
Comment on lines +64 to +83

@devmotion devmotion Jul 27, 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.

AFAIU these tests are failing on master?

@gaoflow gaoflow Jul 27, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Yes — all 32 assertions fail on master (02a51af), worst relative error 1.16e-5 at a=25, x=24.97525, shrinking to ~7.8e-10 by a=1e5. With the one-line change they pass at rtol=1e-13.

The red CI isn't them, though: that's the pre-existing gamma_inc allocations testset on Julia pre. It fails 12 times on master too — https://github.com/JuliaMath/SpecialFunctions.jl/actions/runs/27935023667 on 02a51af (2026-06-22) has the same 12 failures at test/gamma_inc.jl:301/303/305/310, all three OSes, everything else green. Here they show as :328/:330/:332/:337 only because this PR inserts 27 lines above them.

@gaoflow gaoflow Jul 31, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Done — bumped to 2.8.1 in 8a25843, ready to tag after merge.

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