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.

Author: gaoflow

src/gamma_inc.jl

@@ -659,6 +659,19 @@ Used instead of it's previous function when ``\sigma \leq e_{0}/\sqrt{a}``.
 
 External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8)
 """
+# Combine the Temme series value `t` with the Gaussian-tail term `c*w` into (P, Q);
+# shared by gamma_inc_temme_1's three branches and the ind=2 wide-band case in
+# _gamma_inc, which recomputes c/w/t itself instead of calling gamma_inc_temme_1.
+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
+
 function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer)
     iop = ind + 1
     l = x/a
@@ -689,13 +702,7 @@ function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer)
         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,13 +926,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])
-                    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
+                    return _gamma_inc_temme_combine(c, w, t, rta, l)
                 end
             else
                 return _gamma_inc_choose_algorithm(a, x, ind)

test/gamma_inc.jl

@@ -113,17 +113,17 @@ end
         (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
+        @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=5e-4
-        @test gamma_inc(a, x, 2)[2] ≈ q rtol=5e-4
+        @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),