From 52587f51aeedce7be918cf118e0f54921a457b18 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 7 Nov 2022 09:34:18 -0500 Subject: [PATCH 1/2] improve accuracy of `besselk_power_series` `sinpi(1-x)==sinpi(x)` --- src/besselk.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index d070bfe..9d2224f 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -453,8 +453,7 @@ function besselk_power_series(v, x::ComplexOrReal{T}) where T # use the reflection identify to calculate gamma(-v) # use relation gamma(v)*v = gamma(v+1) to avoid two gamma calls gam_v = gamma(v) - xp1 = abs(v) + one(S) - gam_nv = π / (sinpi(xp1) * gam_v * v) + gam_nv = π / (sinpi(v) * gam_v * v) gam_1mv = -gam_nv * v gam_1mnv = gam_v * v From 9b06a72112709dc3b80e382835747bbdf84051fd Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 7 Nov 2022 10:42:44 -0500 Subject: [PATCH 2/2] fix for positive `v` --- src/besselk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/besselk.jl b/src/besselk.jl index 9d2224f..fa8e750 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -453,7 +453,7 @@ function besselk_power_series(v, x::ComplexOrReal{T}) where T # use the reflection identify to calculate gamma(-v) # use relation gamma(v)*v = gamma(v+1) to avoid two gamma calls gam_v = gamma(v) - gam_nv = π / (sinpi(v) * gam_v * v) + gam_nv = π / (sinpi(-abs(v)) * gam_v * v) gam_1mv = -gam_nv * v gam_1mnv = gam_v * v