diff --git a/Project.toml b/Project.toml index adc14fe6..5d3f0419 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SpecialFunctions" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.8.0" +version = "2.8.1" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index be9759d0..2be2a788 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -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) diff --git a/test/gamma_inc.jl b/test/gamma_inc.jl index e2bc545b..6d3679d1 100644 --- a/test/gamma_inc.jl +++ b/test/gamma_inc.jl @@ -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 + 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