Incorporate further improvement from @helto

Author: BioTurboNick

src/gamma_inc.jl

@@ -421,35 +421,16 @@ See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
 """
 function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer)
     acc = acc0[ind + 1]
-    wk = zeros(20)
-    apn = a
-
-    # compute and store larger terms in wk, to add from small to large
-    t = 1
-    i = 0
-    while i < 20
-        i += 1
-        apn += 1.0
-        t *= x / apn
-        t > 1.0e-3 || break
-        wk[i] = t
-    end
-
-    # sum the smaller terms directly
-    sm = t
+    MaxIter = 5000
+    t = 1.0
+    s = 0.0
     tolerance = 0.5acc
-    while t > tolerance
-        apn += 1.0
-        t *= x / apn
-        sm += t
-    end
-
-    # sum terms from small to large
-    for v ∈ @view wk[i:-1:1]
-        sm += v
+    for i in 1:MaxIter
+        s += t
+        t < tolerance && break
+        t *= x / (a + i)
     end
-    
-    p = (rgammax(a, x) / a) * (1.0 + sm)
+    p = (rgammax(a, x) / a) * s
     return (p, 1.0 - p)
 end
 """
@@ -467,34 +448,15 @@ See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
 """
 function gamma_inc_asym(a::Float64, x::Float64, ind::Integer)
     acc = acc0[ind + 1]
-    wk = zeros(20)
-    amn = a
-
-    # compute and store larger terms in wk, to add from small to large
-    t = 1
-    i = 0
-    while i < 20
-        i += 1
-        amn -= 1.0
-        t *= amn / x
-        abs(t) > 1.0e-3 || break
-        wk[i] = t
-    end
-
-    # sum the smaller terms directly
-    sm = t
-    while abs(t) > acc
-        amn -= 1.0
-        t *= amn / x
-        sm += t
+    MaxIter = 5000
+    t = 1.0
+    s = 0.0
+    for i in 1:MaxIter
+        s += t
+        abs(t) < acc && break
+        t *= (a - i) / x
     end
-
-    # sum terms from small to large
-    for v ∈ @view wk[i:-1:1]
-        sm += v
-    end
-
-    q = (rgammax(a, x) / x) * (1.0 + sm)
+    q = (rgammax(a, x) / x) * s
     return (1.0 - q, q)
 end
 """