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.

Author: gaoflow

src/gamma_inc.jl

@@ -681,8 +681,8 @@ 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
@@ -919,7 +919,13 @@ 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)
+                    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
             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=5e-4
+        @test gamma_inc(a, x, 2)[2] ≈ q rtol=5e-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=5e-4
+        @test gamma_inc(a, x, 2)[2] ≈ q rtol=5e-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