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

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

The iop=2 branch expanded c0/c1 in d00/d10 instead of z (NSWC gratio label
270), dropping the z-dependence of the expansion, and the ind=2 wide band
discarded its computed 4-term expansion (label 230) in favour of the
narrow-band 2-term form (label 280). Add reference-point regression tests
for both branches.

* refactor: share the Temme (P,Q) combine step; tighten wide/narrow-band tolerances

Extract the c/w/t -> (P,Q) combine into _gamma_inc_temme_combine, used by
both gamma_inc_temme_1 and the ind=2 wide-band branch in _gamma_inc that
previously duplicated it inline.

Measured actual error against the mpmath references: wide-band tops out
at 1.3e-4, narrow-band at 4.6e-5, so tighten rtol from 5e-4 to 2e-4/1e-4
respectively (~1.5-2x margin) to catch regressions sooner.

* Extend the Temme (P,Q) combine helper to all 5 call sites

Fixes the docstring being pushed off gamma_inc_temme_1 by the
helper's insertion point; also folds in gamma_inc_minimax and
gamma_inc_temme, which had the same duplicated block.

* Bump version from 2.8.1 to 2.8.2

---------

Co-authored-by: David Müller-Widmann <devmotion@users.noreply.github.com>

Author: gaoflow

Project.toml

@@ -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"

src/gamma_inc.jl

@@ -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)
 
@@ -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))
@@ -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"""
@@ -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"""
@@ -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
 
 """
@@ -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)

test/gamma_inc.jl

@@ -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
+    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