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

* Expand the minimax Temme series in 1/a, not z

gamma_inc_minimax evaluates the |1 - x/a| <= 1e-3 branch of

    T(a,lambda) = sum_k c_k(z) a^-k

with @evalpoly(z, c0, ..., d80).  The outer sum is in a^-k, as the
function's own docstring states and as the other three code paths
built from the same d tables already do: the rational branch ten lines
below uses @evalpoly(1.0/a, ...), gamma_inc_temme uses @evalpoly(1.0/a, ...)
and gamma_inc_temme_1 uses @evalpoly(u, ...) with u = 1.0/a.

In the affected band |z| is 2-5% of 1/a, so every correction term past
c0 is weighted 20-40x too small.  Over 300 points with a in [25, 1.2e6]
the worst relative error drops from 1.16e-5 to 8.2e-14; afterwards it
sits at or below the 1-ulp conditioning limit of x.  No point regresses,
and every other branch is bit-identical.

The band was covered by exactly one existing assertion, gamma_inc(1e7,
1e7 + 1), where |z| and 1/a happen to coincide to 1e-7 and the defect
cancels.  The new testset pins 16 points across the band against mpmath.

* Bump version to 2.8.1

Author: gaoflow

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"

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)

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