From eb71c65a0d21d37f3d8ebacaefe5f5c9f8b448df Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Thu, 4 May 2023 01:07:20 -0400 Subject: [PATCH 01/15] Rewrite more concisely --- src/gamma_inc.jl | 142 +++++++++++++++++++++++++++-------------------- 1 file changed, 83 insertions(+), 59 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 3595ecaa..9a5ad041 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -421,39 +421,32 @@ 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(30) - flag = false - apn = a + 1.0 - t = x/apn - wk[1] = t - loop = 2 - for indx = 2:20 - apn += 1.0 - t *= x/apn - if t <= 1.0e-3 - loop = indx - flag = true - break - end - wk[indx] = t - end - if !flag - loop = 20 + 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 - tol = 0.5*acc #tolerance - while true - apn += 1.0 - t *= x/apn - sm += t - if t <= tol - break - end - end - for j = loop-1:-1:1 - sm += wk[j] + tolerance = 0.5acc + while t > tolerance + apn += 1.0 + t *= x / apn + sm += t end - p = (rgammax(a, x)/a)*(1.0 + sm) + + # Sum terms from small to large + sm += sum(reverse!(@view wk[1:i])) + p = (rgammax(a, x) / a) * (1.0 + sm) return (p, 1.0 - p) end """ @@ -470,40 +463,71 @@ External links: [DLMF](https://dlmf.nist.gov/8.11.2) See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc) """ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) - wk = zeros(30) - flag = false acc = acc0[ind + 1] - amn = a - 1.0 - t = amn/x - wk[1] = t - loop = 2 - for indx = 2:20 - amn -= 1.0 - t *= amn/x - if abs(t) <= 1.0e-3 - loop = indx - flag = true - break - end - wk[indx] = t - end - if !flag - loop = 20 + 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 true - if abs(t) < acc - break - end - amn -= 1.0 - t *= amn/x - sm += t - end - for j = loop-1:-1:1 - sm += wk[j] + while abs(t) > acc + amn -= 1.0 + t *= amn / x + sm += t end - q = (rgammax(a, x)/x)*(1.0 + sm) + + # Sum terms from small to large + sm += sum(reverse!(@view wk[1:i])) + q = (rgammax(a, x) / x) * (1.0 + sm) return (1.0 - q, q) + + + + + # wk = zeros(30) + # flag = false + # acc = acc0[ind + 1] + # amn = a - 1.0 + # t = amn/x + # wk[1] = t + # loop = 2 + # for indx = 2:20 + # amn -= 1.0 + # t *= amn/x + # if abs(t) <= 1.0e-3 + # loop = indx + # flag = true + # break + # end + # wk[indx] = t + # end + # if !flag + # loop = 20 + # end + # sm = t + # while true + # if abs(t) < acc + # break + # end + # amn -= 1.0 + # t *= amn/x + # sm += t + # end + # for j = loop-1:-1:1 + # sm += wk[j] + # end + # q = (rgammax(a, x)/x)*(1.0 + sm) + # return (1.0 - q, q) end """ gamma_inc_taylor_x(a,x,ind) From 2363f8f5dbab579e20281117d43b73f8516be40f Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Thu, 4 May 2023 01:10:32 -0400 Subject: [PATCH 02/15] Remove old code --- src/gamma_inc.jl | 38 -------------------------------------- 1 file changed, 38 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 9a5ad041..d0fb8eb2 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -490,44 +490,6 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) sm += sum(reverse!(@view wk[1:i])) q = (rgammax(a, x) / x) * (1.0 + sm) return (1.0 - q, q) - - - - - # wk = zeros(30) - # flag = false - # acc = acc0[ind + 1] - # amn = a - 1.0 - # t = amn/x - # wk[1] = t - # loop = 2 - # for indx = 2:20 - # amn -= 1.0 - # t *= amn/x - # if abs(t) <= 1.0e-3 - # loop = indx - # flag = true - # break - # end - # wk[indx] = t - # end - # if !flag - # loop = 20 - # end - # sm = t - # while true - # if abs(t) < acc - # break - # end - # amn -= 1.0 - # t *= amn/x - # sm += t - # end - # for j = loop-1:-1:1 - # sm += wk[j] - # end - # q = (rgammax(a, x)/x)*(1.0 + sm) - # return (1.0 - q, q) end """ gamma_inc_taylor_x(a,x,ind) From 0b0a61a446869132ebf0885ea5cf39492cecfcd5 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Thu, 4 May 2023 07:30:08 -0400 Subject: [PATCH 03/15] Remove sum (does pairwise addition?) --- src/gamma_inc.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index d0fb8eb2..5a8b8f5c 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -444,8 +444,11 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) sm += t end - # Sum terms from small to large - sm += sum(reverse!(@view wk[1:i])) + # sum terms from small to large + for v ∈ @view wk[i:-1:1] + sm += v + end + p = (rgammax(a, x) / a) * (1.0 + sm) return (p, 1.0 - p) end @@ -486,8 +489,11 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) sm += t end - # Sum terms from small to large - sm += sum(reverse!(@view wk[1:i])) + # sum terms from small to large + for v ∈ @view wk[i:-1:1] + sm += v + end + q = (rgammax(a, x) / x) * (1.0 + sm) return (1.0 - q, q) end From 46a7b39c8e333cc7836421efc5ce5bdd758c7818 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Fri, 5 May 2023 08:29:41 -0400 Subject: [PATCH 04/15] Remove view --- src/gamma_inc.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 5a8b8f5c..14994207 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -425,7 +425,7 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) apn = a # compute and store larger terms in wk, to add from small to large - t = 1 + t = 1.0 i = 0 while i < 20 i += 1 @@ -445,8 +445,8 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) end # sum terms from small to large - for v ∈ @view wk[i:-1:1] - sm += v + for j ∈ i:(-1):1 + sm += wk[j] end p = (rgammax(a, x) / a) * (1.0 + sm) From 4e2edc6ce3acd5be95d156957c148497a9ccfe11 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Fri, 5 May 2023 08:44:01 -0400 Subject: [PATCH 05/15] Missed a couple Co-authored-by: David Widmann --- src/gamma_inc.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 14994207..338b7567 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -471,7 +471,7 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) amn = a # compute and store larger terms in wk, to add from small to large - t = 1 + t = 1.0 i = 0 while i < 20 i += 1 @@ -490,8 +490,8 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) end # sum terms from small to large - for v ∈ @view wk[i:-1:1] - sm += v + for j in i:(-1):1 + sm += wk[j] end q = (rgammax(a, x) / x) * (1.0 + sm) From 7ac105ea6b9e029add028950f687e496bf575b52 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 9 May 2023 18:35:56 -0400 Subject: [PATCH 06/15] Move to the stack! --- src/gamma_inc.jl | 52 +++++++++++++++++++----------------------------- 1 file changed, 20 insertions(+), 32 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 338b7567..adb19cb6 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -421,22 +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.0 - i = 0 - while i < 20 - i += 1 - apn += 1.0 - t *= x / apn - t > 1.0e-3 || break - wk[i] = t - end + + # compute first 21 terms + ts = cumprod(ntuple(i -> x / (a + i), Val(21))) + last_t = findfirst(<(1.0e-3), ts) + last_t = last_t === nothing ? 20 : last_t - 1 + next_t = last_t + 1 # sum the smaller terms directly - sm = t + sm = t = ts[next_t] + apn = a + next_t tolerance = 0.5acc while t > tolerance apn += 1.0 @@ -445,8 +439,8 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) end # sum terms from small to large - for j ∈ i:(-1):1 - sm += wk[j] + for j ∈ last_t:(-1):1 + sm += ts[j] end p = (rgammax(a, x) / a) * (1.0 + sm) @@ -467,22 +461,16 @@ 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.0 - i = 0 - while i < 20 - i += 1 - amn -= 1.0 - t *= amn / x - abs(t) > 1.0e-3 || break - wk[i] = t - end + + # compute first 21 terms + ts = cumprod(ntuple(i -> (a - i) / x, Val(21))) + last_t = findfirst(x -> abs(x) < 1.0e-3, ts) + last_t = last_t === nothing ? 20 : last_t - 1 + next_t = last_t + 1 # sum the smaller terms directly - sm = t + sm = t = ts[next_t] + amn = a - next_t while abs(t) > acc amn -= 1.0 t *= amn / x @@ -490,8 +478,8 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) end # sum terms from small to large - for j in i:(-1):1 - sm += wk[j] + for j in last_t:(-1):1 + sm += ts[j] end q = (rgammax(a, x) / x) * (1.0 + sm) From c0953297621273a0db662c7238196463f77b0632 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 10:02:41 -0400 Subject: [PATCH 07/15] Try to add cumprod for compat with 1.3 --- src/gamma_inc.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index adb19cb6..c130bd01 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -485,6 +485,9 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) q = (rgammax(a, x) / x) * (1.0 + sm) return (1.0 - q, q) end +if VERSION < v"1.5.0" + cumprod(itr) = accumulate(mul_prod, itr) +end """ gamma_inc_taylor_x(a,x,ind) From c058a2ac0e77a359b5f47e4c1ea2e9889435e7bd Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 10:03:18 -0400 Subject: [PATCH 08/15] Add Base ref --- src/gamma_inc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index c130bd01..a0c8ec1c 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -486,7 +486,7 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) return (1.0 - q, q) end if VERSION < v"1.5.0" - cumprod(itr) = accumulate(mul_prod, itr) + Base.cumprod(itr) = accumulate(mul_prod, itr) end """ gamma_inc_taylor_x(a,x,ind) From ee41ab714f4b791221a420eef3e4ac9210907679 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 10:18:16 -0400 Subject: [PATCH 09/15] Add Base ref --- src/gamma_inc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index a0c8ec1c..6f1a2830 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -486,7 +486,7 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) return (1.0 - q, q) end if VERSION < v"1.5.0" - Base.cumprod(itr) = accumulate(mul_prod, itr) + Base.cumprod(itr) = accumulate(Base.mul_prod, itr) end """ gamma_inc_taylor_x(a,x,ind) From 83afebbb37819bb2574f6dcd970a69a0bd4c028c Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 11:34:43 -0400 Subject: [PATCH 10/15] Revert compat attempts --- src/gamma_inc.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 6f1a2830..adb19cb6 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -485,9 +485,6 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) q = (rgammax(a, x) / x) * (1.0 + sm) return (1.0 - q, q) end -if VERSION < v"1.5.0" - Base.cumprod(itr) = accumulate(Base.mul_prod, itr) -end """ gamma_inc_taylor_x(a,x,ind) From 34229ef341d73461b3995a734b41e8b883710f98 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 20:43:32 -0400 Subject: [PATCH 11/15] Bump CI --- src/gamma_inc.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index adb19cb6..ed6c2470 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -444,6 +444,7 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) end p = (rgammax(a, x) / a) * (1.0 + sm) + return (p, 1.0 - p) end """ @@ -481,8 +482,9 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) for j in last_t:(-1):1 sm += ts[j] end - + q = (rgammax(a, x) / x) * (1.0 + sm) + return (1.0 - q, q) end """ From 482516ce2c43b37e1916c78572061ea6b767460e Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 20:45:42 -0400 Subject: [PATCH 12/15] Change minimum supported Julia from 1.3 to 1.5 --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c9493aaa..5928260b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -14,7 +14,7 @@ jobs: fail-fast: false matrix: version: - - "1.3" + - "1.5" - "1" - "nightly" os: From 8d37a34c1df9635f57250eba175bb292641b9c42 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 10 May 2023 20:46:06 -0400 Subject: [PATCH 13/15] Change minimum supported Julia to 1.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cf84beb2..682973ee 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,7 @@ IrrationalConstants = "0.1, 0.2" LogExpFunctions = "0.3.2" OpenLibm_jll = "0.7, 0.8" OpenSpecFun_jll = "0.5" -julia = "1.3" +julia = "1.5" [extras] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 79bacae307ffc83cc7e158e24370798cfc03a62b Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Mon, 28 Oct 2024 15:30:13 -0400 Subject: [PATCH 14/15] Proper inversion of test Co-authored-by: CyHan --- src/gamma_inc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index 9af28f77..f2fd5c2b 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -478,7 +478,7 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) # sum the smaller terms directly sm = t = ts[next_t] amn = a - next_t - while abs(t) > acc + while abs(t) >= acc amn -= 1.0 t *= amn / x sm += t From 1c41a3db13b0306919e6ac0fef9414ea2e8c50f2 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Mon, 28 Oct 2024 21:33:31 -0400 Subject: [PATCH 15/15] Simplify/clarify --- src/gamma_inc.jl | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/src/gamma_inc.jl b/src/gamma_inc.jl index f2fd5c2b..7fbc9fb0 100644 --- a/src/gamma_inc.jl +++ b/src/gamma_inc.jl @@ -426,17 +426,15 @@ 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] + tolerance = 0.5acc # compute first 21 terms ts = cumprod(ntuple(i -> x / (a + i), Val(21))) - last_t = findfirst(<(1.0e-3), ts) - last_t = last_t === nothing ? 20 : last_t - 1 - next_t = last_t + 1 - + # sum the smaller terms directly - sm = t = ts[next_t] - apn = a + next_t - tolerance = 0.5acc + first_small_t = something(findfirst(<(1.0e-3), ts), 21) + sm = t = ts[first_small_t] + apn = a + first_small_t while t > tolerance apn += 1.0 t *= x / apn @@ -444,7 +442,8 @@ function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer) end # sum terms from small to large - for j ∈ last_t:(-1):1 + last_large_t = first_small_t - 1 + for j ∈ last_large_t:(-1):1 sm += ts[j] end @@ -471,21 +470,20 @@ function gamma_inc_asym(a::Float64, x::Float64, ind::Integer) # compute first 21 terms ts = cumprod(ntuple(i -> (a - i) / x, Val(21))) - last_t = findfirst(x -> abs(x) < 1.0e-3, ts) - last_t = last_t === nothing ? 20 : last_t - 1 - next_t = last_t + 1 # sum the smaller terms directly - sm = t = ts[next_t] - amn = a - next_t - while abs(t) >= acc + first_small_t = something(findfirst(x -> abs(x) < 1.0e-3, ts), 21) + sm = t = ts[first_small_t] + amn = a - first_small_t + while abs(t) ≥ acc amn -= 1.0 t *= amn / x sm += t end # sum terms from small to large - for j in last_t:(-1):1 + last_large_t = first_small_t - 1 + for j in last_large_t:(-1):1 sm += ts[j] end