diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 5a8b8f5c..3e10fc67 100644 --- a/src/gamma_inc.jl +++ b/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 """