From aeddfe77733fbb6bf0f81699fa5b6c499c10abac Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 12:14:40 -0400 Subject: [PATCH 01/14] continued fraction --- src/U_polynomials.jl | 8 ++--- src/besseli.jl | 42 ++++++++++---------------- src/besselk.jl | 71 ++++++++++++++++++++++++++++++++++++++++++-- src/recurrence.jl | 4 +-- 4 files changed, 90 insertions(+), 35 deletions(-) diff --git a/src/U_polynomials.jl b/src/U_polynomials.jl index 0096d2f..b9cf29c 100644 --- a/src/U_polynomials.jl +++ b/src/U_polynomials.jl @@ -70,10 +70,10 @@ function Uk_poly_Hankel(p, v, p2, x::T) where T <: Float64 end Uk_poly_Hankel(p, v, p2, x) = Uk_poly_Jn(p, v, p2, x::BigFloat) -Uk_poly_In(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[1] -Uk_poly_In(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[1] -Uk_poly_Kn(p, v, p2, ::Type{Float32}) = Uk_poly3(p, v, p2)[2] -Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly5(p, v, p2)[2] +Uk_poly_In(p, v, p2, ::Type{Float32}) = Uk_poly5(p, v, p2)[1] +Uk_poly_In(p, v, p2, ::Type{Float64}) = Uk_poly10(p, v, p2)[1] +Uk_poly_Kn(p, v, p2, ::Type{Float32}) = Uk_poly5(p, v, p2)[2] +Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly10(p, v, p2)[2] @inline function split_evalpoly(x, P) # polynomial P must have an even number of terms diff --git a/src/besseli.jl b/src/besseli.jl index 1841aeb..fd99c49 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -139,18 +139,15 @@ Nu must be real. function besseli(nu, x::T) where T <: Union{Float32, Float64} nu == 0 && return besseli0(x) nu == 1 && return besseli1(x) - - if x > maximum((T(30), nu^2 / 4)) + + if x > maximum((T(30), nu^2 / 6)) return T(besseli_large_argument(nu, x)) - elseif x <= 2 * sqrt(nu + 1) - return T(besseli_small_arguments(nu, x)) - elseif nu < 100 - return T(_besseli_continued_fractions(nu, x)) - else + elseif nu > 25.0 || x > 35.0 return T(besseli_large_orders(nu, x)) + else + return T(besseli_power_series(nu, x)) end end - """ besselix(nu, x::T) where T <: Union{Float32, Float64} @@ -264,25 +261,16 @@ function besseli_large_argument_scaled(v, z::T) where T return res * coef end -function besseli_small_arguments(v, z::T) where T - S = promote_type(T, Float64) - x = S(z) - if v < 20 - coef = (x / 2)^v / factorial(v) - else - vinv = inv(v) - coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1)) - coef *= evalpoly(vinv, (1, -1/12, 1/288, 139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600)) - end - - MaxIter = 1000 - out = one(S) - zz = x^2 / 4 - a = one(S) - for k in 1:MaxIter - a *= zz / (k * (k + v)) +function besseli_power_series(v, x::T) where T + MaxIter = 3000 + out = zero(T) + xs = (x/2)^v + a = xs / gamma(v + one(T)) + t2 = (x/2)^2 + for i in 0:MaxIter out += a - a <= eps(T) && break + abs(a) < eps(T) * abs(out) && break + a *= inv((v + i + one(T)) * (i + one(T))) * t2 end - return coef * out + return out end diff --git a/src/besselk.jl b/src/besselk.jl index 7ca1dab..5baa4b2 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -153,6 +153,7 @@ In other words, for large values of x and small to moderate values of nu, this routine will underflow to zero. For small to medium values of x and large values of nu this will overflow and return Inf. =# +#= """ besselk(x::T) where T <: Union{Float32, Float64} @@ -166,7 +167,7 @@ function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} return besselk_large_orders(nu, x) end end - +=# """ besselk(x::T) where T <: Union{Float32, Float64} @@ -203,4 +204,70 @@ function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, p2 = v^2/fma(max(v,x), max(v,x), min(v,x)^2) return T(coef*Uk_poly_Kn(p, v, p2, T)) -end \ No newline at end of file +end + +function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} + (isinteger(nu) && nu < 250) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1] + + if nu > 25.0 || x > 35.0 + return besselk_large_orders(nu, x) + else + return besselk_continued_fraction(nu, x) + end +end + +# could also use the continued fraction for inu/inmu1 +# but it seems like adapting the besseli_power series +# to give both nu and nu-1 is faster + +#inum1 = besseli_power_series(nu-1, x) +#H_inu = steed(nu, x) +#inu = besseli_power_series(nu, x)#inum1 * H_inu +function besselk_continued_fraction(nu, x) + inu, inum1 = besseli_power_series_inu_inum1(nu, x) + H_knu = besselk_ratio_knu_knup1(nu-1, x) + return 1 / (x * (inum1 + inu / H_knu)) +end + +function besseli_power_series_inu_inum1(v, x::T) where T + MaxIter = 3000 + out = zero(T) + out2 = zero(T) + x2 = x / 2 + xs = (x2)^v + gmx = xs / gamma(v) + a = gmx / v + b = gmx / x2 + t2 = x2*x2 + for i in 0:MaxIter + out += a + out2 += b + abs(a) < eps(T) * abs(out) && break + a *= inv((v + i + one(T)) * (i + one(T))) * t2 + b *= inv((v + i) * (i + one(T))) * t2 + end + return out, out2 +end + + +# slightly modified version of https://github.com/heltonmc/Bessels.jl/issues/17#issuecomment-1195726642 from @cgeoga +function besselk_ratio_knu_knup1(v, x::T) where T + MaxIter = 1000 + # do the modified Lentz method: + (hn, Dn, Cn) = (1e-50, zero(v), 1e-50) + + jf = one(T) + vv = v*v + for _ in 1:MaxIter + an = (vv - ((2*jf - 1)^2) * T(0.25)) + bn = 2 * (x + jf) + Cn = an / Cn + bn + Dn = inv(muladd(an, Dn, bn)) + del = Dn * Cn + hn *= del + abs(del - 1) < eps(T) && break + jf += one(T) + end + xinv = inv(x) + return xinv * (v + x + 1/2) + xinv * hn +end diff --git a/src/recurrence.jl b/src/recurrence.jl index fd808ca..f32c611 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -36,7 +36,7 @@ end return jnup1, jnu end -#= + # currently not used # backward recurrence relation for besselk and besseli # outputs both (bessel(x, nu_end), bessel(x, nu_end-1) @@ -50,4 +50,4 @@ end end return jnup1, jnu end -=# + From 9f6cb5cbc72219dda76176355e6998015c454735 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 14:23:57 -0400 Subject: [PATCH 02/14] add power series Co-authored-by: Chris Geoga --- src/besselk.jl | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/src/besselk.jl b/src/besselk.jl index 5baa4b2..e3b5813 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -271,3 +271,52 @@ function besselk_ratio_knu_knup1(v, x::T) where T xinv = inv(x) return xinv * (v + x + 1/2) + xinv * hn end + +# originally contribued by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl +function _besselk_ser(v, x, maxit, tol, modify) + T = promote_type(typeof(x), typeof(v)) + out = zero(T) + oneT = one(T) + twoT = oneT+oneT + # precompute a handful of things: + xd2 = x/twoT + xd22 = xd2*xd2 + half = oneT/twoT + if modify + e2v = exp2(v) + xd2_v = (x^(2*v))/e2v + xd2_nv = e2v + else + lxd2 = log(xd2) + xd2_v = exp(v*lxd2) + xd2_nv = exp(-v*lxd2) + end + gam_v = gamma(v) + gam_nv = gamma(-v) + gam_1mv = -gam_nv*v # == gamma(one(T)-v) + gam_1mnv = gam_v*v # == gamma(one(T)+v) + xd2_pow = oneT + fact_k = oneT + floatk = convert(T, 0) + (gpv, gmv) = (gam_1mnv, gam_1mv) + # One final re-compression of a few things: + _t1 = gam_v*xd2_nv*gam_1mv + _t2 = gam_nv*xd2_v*gam_1mnv + # now the loop using Oana's series expansion, with term function manually + # inlined for max speed: + for _j in 0:maxit + t1 = half*xd2_pow + tmp = _t1/(gmv*fact_k) + tmp += _t2/(gpv*fact_k) + term = t1*tmp + out += term + ((abs(term) < tol) && _j>5) && return out + # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls: + (gpv, gmv) = (gpv*(oneT+v+floatk), gmv*(oneT-v+floatk)) + xd2_pow *= xd22 + fact_k *= (floatk+1) + floatk += 1.0 + end + throw(error("$maxit iterations reached without achieving atol $tol.")) +end + From 6863aeb164bcfdc03bd2fe7877b5ac93bdfffb73 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 14:47:38 -0400 Subject: [PATCH 03/14] update power series --- src/besselk.jl | 89 ++++++++++++++++++++++++-------------------------- 1 file changed, 42 insertions(+), 47 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index e3b5813..c90d425 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -211,6 +211,8 @@ function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} if nu > 25.0 || x > 35.0 return besselk_large_orders(nu, x) + elseif x < 2.0 || nu > 1.6x - 1.0 + return besselk_power_series(nu, x) else return besselk_continued_fraction(nu, x) end @@ -272,51 +274,44 @@ function besselk_ratio_knu_knup1(v, x::T) where T return xinv * (v + x + 1/2) + xinv * hn end -# originally contribued by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl -function _besselk_ser(v, x, maxit, tol, modify) - T = promote_type(typeof(x), typeof(v)) - out = zero(T) - oneT = one(T) - twoT = oneT+oneT - # precompute a handful of things: - xd2 = x/twoT - xd22 = xd2*xd2 - half = oneT/twoT - if modify - e2v = exp2(v) - xd2_v = (x^(2*v))/e2v - xd2_nv = e2v - else - lxd2 = log(xd2) - xd2_v = exp(v*lxd2) - xd2_nv = exp(-v*lxd2) - end - gam_v = gamma(v) - gam_nv = gamma(-v) - gam_1mv = -gam_nv*v # == gamma(one(T)-v) - gam_1mnv = gam_v*v # == gamma(one(T)+v) - xd2_pow = oneT - fact_k = oneT - floatk = convert(T, 0) - (gpv, gmv) = (gam_1mnv, gam_1mv) - # One final re-compression of a few things: - _t1 = gam_v*xd2_nv*gam_1mv - _t2 = gam_nv*xd2_v*gam_1mnv - # now the loop using Oana's series expansion, with term function manually - # inlined for max speed: - for _j in 0:maxit - t1 = half*xd2_pow - tmp = _t1/(gmv*fact_k) - tmp += _t2/(gpv*fact_k) - term = t1*tmp - out += term - ((abs(term) < tol) && _j>5) && return out - # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls: - (gpv, gmv) = (gpv*(oneT+v+floatk), gmv*(oneT-v+floatk)) - xd2_pow *= xd22 - fact_k *= (floatk+1) - floatk += 1.0 +# originally contributed by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl +# Equation 3.2 from Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." +# arXiv preprint arXiv:2201.00090 (2022). +function besselk_power_series(v, x::T) where T + MaxIter = 1000 + # precompute a handful of things: + xd2 = x / 2 + xd22 = xd2 * xd2 + half = one(T) / 2 + # (x/2)^(±v). Writing the literal power doesn't seem to change anything here, + # and I think this is faster: + lxd2 = log(xd2) + xd2_v = exp(v*lxd2) + xd2_nv = exp(-v*lxd2) + # use the gamma function a couple times to start: + gam_v = gamma(v) + xp1 = abs(v) + one(T) + gam_nv = π / sinpi(xp1) / _gamma(xp1) + gam_1mv = -gam_nv*v # == gamma(one(T)-v) + gam_1mnv = gam_v*v # == gamma(one(T)+v) + (gpv, gmv) = (gam_1mnv, gam_1mv) + # One final re-compression of a few things: + _t1 = gam_v*xd2_nv*gam_1mv + _t2 = gam_nv*xd2_v*gam_1mnv + # A couple series-specific accumulators: + (xd2_pow, fact_k, floatk, out) = (one(T), one(T), zero(T), zero(T)) + for _ in 0:MaxIter + t1 = half*xd2_pow + tmp = _t1/(gmv*fact_k) + tmp += _t2/(gpv*fact_k) + term = t1*tmp + out += term + abs(term/out) < eps(T) && return out + # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls: + (gpv, gmv) = (gpv*(one(T)+v+floatk), gmv*(one(T)-v+floatk)) + xd2_pow *= xd22 + fact_k *= (floatk+1) + floatk += one(T) + end + return out end - throw(error("$maxit iterations reached without achieving atol $tol.")) -end - From 71de923b0143ead2b77db80ea4f548139bfa51d8 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 14:52:02 -0400 Subject: [PATCH 04/14] rm extra assignment --- src/besselk.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index c90d425..c02cd8b 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -294,7 +294,7 @@ function besselk_power_series(v, x::T) where T gam_nv = π / sinpi(xp1) / _gamma(xp1) gam_1mv = -gam_nv*v # == gamma(one(T)-v) gam_1mnv = gam_v*v # == gamma(one(T)+v) - (gpv, gmv) = (gam_1mnv, gam_1mv) + # One final re-compression of a few things: _t1 = gam_v*xd2_nv*gam_1mv _t2 = gam_nv*xd2_v*gam_1mnv @@ -302,13 +302,13 @@ function besselk_power_series(v, x::T) where T (xd2_pow, fact_k, floatk, out) = (one(T), one(T), zero(T), zero(T)) for _ in 0:MaxIter t1 = half*xd2_pow - tmp = _t1/(gmv*fact_k) - tmp += _t2/(gpv*fact_k) + tmp = _t1/(gam_1mv*fact_k) + tmp += _t2/(gam_1mnv*fact_k) term = t1*tmp out += term abs(term/out) < eps(T) && return out # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls: - (gpv, gmv) = (gpv*(one(T)+v+floatk), gmv*(one(T)-v+floatk)) + (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T)+v+floatk), gam_1mv*(one(T)-v+floatk)) xd2_pow *= xd22 fact_k *= (floatk+1) floatk += one(T) From 3252399a6c84012c7c9d3d3a3cfa88079cb62cd7 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 14:57:59 -0400 Subject: [PATCH 05/14] rm extra division --- src/besselk.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index c02cd8b..71ec9e4 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -279,7 +279,6 @@ end # arXiv preprint arXiv:2201.00090 (2022). function besselk_power_series(v, x::T) where T MaxIter = 1000 - # precompute a handful of things: xd2 = x / 2 xd22 = xd2 * xd2 half = one(T) / 2 @@ -291,10 +290,10 @@ function besselk_power_series(v, x::T) where T # use the gamma function a couple times to start: gam_v = gamma(v) xp1 = abs(v) + one(T) - gam_nv = π / sinpi(xp1) / _gamma(xp1) + gam_nv = π / (sinpi(xp1) * _gamma(xp1)) gam_1mv = -gam_nv*v # == gamma(one(T)-v) gam_1mnv = gam_v*v # == gamma(one(T)+v) - + # One final re-compression of a few things: _t1 = gam_v*xd2_nv*gam_1mv _t2 = gam_nv*xd2_v*gam_1mnv From 277bd35ca742e975a0ad5ccb26e049cb193f7e60 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 15:16:01 -0400 Subject: [PATCH 06/14] actually optimize gamma call away --- src/besselk.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 71ec9e4..fc2791d 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -211,7 +211,7 @@ function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} if nu > 25.0 || x > 35.0 return besselk_large_orders(nu, x) - elseif x < 2.0 || nu > 1.6x - 1.0 + elseif x < 2.0 return besselk_power_series(nu, x) else return besselk_continued_fraction(nu, x) @@ -287,12 +287,15 @@ function besselk_power_series(v, x::T) where T lxd2 = log(xd2) xd2_v = exp(v*lxd2) xd2_nv = exp(-v*lxd2) - # use the gamma function a couple times to start: + gam_v = gamma(v) + # use the reflection identify to calculate gamma(-v) + # use relation gamma(v)*v = gamma(v+1) xp1 = abs(v) + one(T) - gam_nv = π / (sinpi(xp1) * _gamma(xp1)) - gam_1mv = -gam_nv*v # == gamma(one(T)-v) - gam_1mnv = gam_v*v # == gamma(one(T)+v) + gam_nv = π / (sinpi(xp1) * gam_v*v) + + gam_1mv = -gam_nv*v + gam_1mnv = gam_v*v # One final re-compression of a few things: _t1 = gam_v*xd2_nv*gam_1mv From b4a8832825258e56373c7248e5a6c4dce65f1834 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 16:21:17 -0400 Subject: [PATCH 07/14] update docs for power series --- src/besselk.jl | 78 ++++++++++++++++++++++++++++---------------------- 1 file changed, 44 insertions(+), 34 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index fc2791d..8469d69 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -274,46 +274,56 @@ function besselk_ratio_knu_knup1(v, x::T) where T return xinv * (v + x + 1/2) + xinv * hn end -# originally contributed by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl -# Equation 3.2 from Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." -# arXiv preprint arXiv:2201.00090 (2022). + +##### +##### Power series for K_{nu}(x) +##### + +# Use power series form of K_v(x) which is accurate for small x (x<2) or when nu > X +# We use the form as described by Equation 3.2 from reference [7]. +# This method was originally contributed by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl +# A modified form appears below. See more discussion at https://github.com/heltonmc/Bessels.jl/pull/29 +# This is only accurate for noninteger orders (nu) and no checks are performed. +# +# [7] Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." +# arXiv preprint arXiv:2201.00090 (2022). +""" + besselk_power_series(nu, x::T) where T <: Float64 + +Computes ``K_{nu}(x)`` using the power series when nu is not an integer. +In general, this is most accurate for small arguments and when nu > x. +No checks are performed on nu so this is not accurate when nu is an integer. +""" function besselk_power_series(v, x::T) where T MaxIter = 1000 - xd2 = x / 2 - xd22 = xd2 * xd2 - half = one(T) / 2 - # (x/2)^(±v). Writing the literal power doesn't seem to change anything here, - # and I think this is faster: - lxd2 = log(xd2) - xd2_v = exp(v*lxd2) - xd2_nv = exp(-v*lxd2) + z = x / 2 + zz = z * z + + logz = log(z) + xd2_v = exp(v*logz) + xd2_nv = inv(xd2_v) - gam_v = gamma(v) # use the reflection identify to calculate gamma(-v) - # use relation gamma(v)*v = gamma(v+1) + # use relation gamma(v)*v = gamma(v+1) to avoid two gamma calls + gam_v = gamma(v) xp1 = abs(v) + one(T) - gam_nv = π / (sinpi(xp1) * gam_v*v) - - gam_1mv = -gam_nv*v - gam_1mnv = gam_v*v + gam_nv = π / (sinpi(xp1) * gam_v * v) + gam_1mv = -gam_nv * v + gam_1mnv = gam_v * v - # One final re-compression of a few things: - _t1 = gam_v*xd2_nv*gam_1mv - _t2 = gam_nv*xd2_v*gam_1mnv - # A couple series-specific accumulators: - (xd2_pow, fact_k, floatk, out) = (one(T), one(T), zero(T), zero(T)) - for _ in 0:MaxIter - t1 = half*xd2_pow - tmp = _t1/(gam_1mv*fact_k) - tmp += _t2/(gam_1mnv*fact_k) - term = t1*tmp + _t1 = gam_v * xd2_nv * gam_1mv + _t2 = gam_nv * xd2_v * gam_1mnv + (xd2_pow, fact_k, out) = (one(T), one(T), zero(T)) + for k in 0:MaxIter + t1 = xd2_pow * T(0.5) + tmp = muladd(_t1, gam_1mnv, _t2 * gam_1mv) + tmp *= inv(gam_1mv * gam_1mnv * fact_k) + term = t1 * tmp out += term - abs(term/out) < eps(T) && return out - # Use the trick that gamma(1+k+1+v) == gamma(1+k+v)*(1+k+v) to skip gamma calls: - (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T)+v+floatk), gam_1mv*(one(T)-v+floatk)) - xd2_pow *= xd22 - fact_k *= (floatk+1) - floatk += one(T) + abs(term / out) < eps(T) && return out + (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T) + v + k), gam_1mv*(one(T) - v + k)) + xd2_pow *= zz + fact_k *= k + one(T) end return out - end +end From bc4ec0e05855748dd3c7ca65c0fe23bb8141b2b4 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 17:22:59 -0400 Subject: [PATCH 08/14] add docs and besselk tests --- src/besseli.jl | 2 +- src/besselk.jl | 151 ++++++++++++++++++++++++------------------- test/besselk_test.jl | 21 ++++-- 3 files changed, 104 insertions(+), 70 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index fd99c49..f617764 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -161,7 +161,7 @@ function besselix(nu, x::T) where T <: Union{Float32, Float64} if x > maximum((T(30), nu^2 / 4)) return T(besseli_large_argument_scaled(nu, x)) elseif x <= 2 * sqrt(nu + 1) - return T(besseli_small_arguments(nu, x)) * exp(-x) + return T(besseli_power_series(nu, x)) * exp(-x) elseif nu < 100 return T(_besseli_continued_fractions_scaled(nu, x)) else diff --git a/src/besselk.jl b/src/besselk.jl index 8469d69..11e15cd 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -145,43 +145,63 @@ function besselk1x(x::T) where T <: Union{Float32, Float64} return muladd(evalpoly(inv(x), P3_k1(T)), inv(evalpoly(inv(x), Q3_k1(T))), Y2_k1(T)) / sqrt(x) end end -#= -If besselk0(x) or besselk1(0) is equal to zero -this will underflow and always return zero even if besselk(x, nu) -is larger than the smallest representing floating point value. -In other words, for large values of x and small to moderate values of nu, -this routine will underflow to zero. -For small to medium values of x and large values of nu this will overflow and return Inf. -=# -#= + """ besselk(x::T) where T <: Union{Float32, Float64} Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ -function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} - T == Float32 ? branch = 20 : branch = 50 - if nu < branch - return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1] - else - return besselk_large_orders(nu, x) - end +function besselk(nu, x::T) where T <: Union{Float32, Float64} + # check to make sure nu isn't zero + iszero(nu) && return besselk0(x) + + # use uniform debye expansion if x or nu is large + besselik_debye_cutoff(nu, x) && return besselk_large_orders(nu, x) + + # for integer nu use forward recurrence starting with K_0 and K_1 + isinteger(nu) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1] + + # for small x and nu > x use power series + besselk_power_series_cutoff(nu, x) && return besselk_power_series(nu, x) + + # for rest of values use the continued fraction approach + return besselk_continued_fraction(nu, x) end -=# """ - besselk(x::T) where T <: Union{Float32, Float64} + besselkx(x::T) where T <: Union{Float32, Float64} Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x}``. """ -function besselkx(nu::Int, x::T) where T <: Union{Float32, Float64} - T == Float32 ? branch = 20 : branch = 50 - if nu < branch - return besselk_up_recurrence(x, besselk1x(x), besselk0x(x), 1, nu)[1] - else - return besselk_large_orders_scaled(nu, x) - end +function besselkx(nu, x::T) where T <: Union{Float32, Float64} + # check to make sure nu isn't zero + iszero(nu) && return besselk0x(x) + + # use uniform debye expansion if x or nu is large + besselik_debye_cutoff(nu, x) && return besselk_large_orders_scaled(nu, x) + + # for integer nu use forward recurrence starting with K_0 and K_1 + isinteger(nu) && return besselk_up_recurrence(x, besselk1x(x), besselk0x(x), 1, nu)[1] + + # for small x and nu > x use power series + besselk_power_series_cutoff(nu, x) && return besselk_power_series(nu, x) * exp(x) + + # for rest of values use the continued fraction approach + return besselk_continued_fraction(nu, x) * exp(x) end -function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat} + +##### +##### Debye's uniform asymptotic for K_{nu}(x) +##### + +# Implements the uniform asymptotic expansion https://dlmf.nist.gov/10.41 +# In general this is valid when either x or nu is gets large +# see the file src/U_polynomials.jl for more details +""" + besselk_large_orders(nu, x::T) + +Debey's uniform asymptotic expansion for large order valid when v-> ∞ or x -> ∞ +""" +function besselk_large_orders(v, x::T) where T S = promote_type(T, Float64) x = S(x) z = x / v @@ -193,7 +213,7 @@ function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFlo return T(coef*Uk_poly_Kn(p, v, p2, T)) end -function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat} +function besselk_large_orders_scaled(v, x::T) where T S = promote_type(T, Float64) x = S(x) z = x / v @@ -205,42 +225,44 @@ function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, return T(coef*Uk_poly_Kn(p, v, p2, T)) end +besselik_debye_cutoff(nu, x::Float64) = nu > 25.0 || x > 35.0 +besselik_debye_cutoff(nu, x::Float32) = nu > 15.0 || x > 20.0 -function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat} - (isinteger(nu) && nu < 250) && return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1] - - if nu > 25.0 || x > 35.0 - return besselk_large_orders(nu, x) - elseif x < 2.0 - return besselk_power_series(nu, x) - else - return besselk_continued_fraction(nu, x) - end -end +##### +##### Continued fraction with Wronskian for K_{nu}(x) +##### -# could also use the continued fraction for inu/inmu1 -# but it seems like adapting the besseli_power series -# to give both nu and nu-1 is faster +# Use the ratio K_{nu+1}/K_{nu} and I_{nu-1}, I_{nu} +# along with the Wronskian (NIST https://dlmf.nist.gov/10.28.E2) to calculate K_{nu} +# Inu and Inum1 are generated from the power series form where K_{nu_1}/K_{nu} +# is calculated with continued fractions. +# The continued fraction K_{nu_1}/K_{nu} method is a slightly modified form +# https://github.com/heltonmc/Bessels.jl/issues/17#issuecomment-1195726642 by @cgeoga +# +# It is also possible to use continued fraction to calculate inu/inmu1 such as +# inum1 = besseli_power_series(nu-1, x) +# H_inu = steed(nu, x) +# inu = besseli_power_series(nu, x)#inum1 * H_inu +# but it appears to be faster to modify the series to calculate both Inu and Inum1 -#inum1 = besseli_power_series(nu-1, x) -#H_inu = steed(nu, x) -#inu = besseli_power_series(nu, x)#inum1 * H_inu function besselk_continued_fraction(nu, x) inu, inum1 = besseli_power_series_inu_inum1(nu, x) H_knu = besselk_ratio_knu_knup1(nu-1, x) return 1 / (x * (inum1 + inu / H_knu)) end +# a modified version of the I_{nu} power series to compute both I_{nu} and I_{nu-1} +# use this along with the continued fractions for besselk function besseli_power_series_inu_inum1(v, x::T) where T MaxIter = 3000 out = zero(T) out2 = zero(T) x2 = x / 2 - xs = (x2)^v + xs = x2^v gmx = xs / gamma(v) a = gmx / v b = gmx / x2 - t2 = x2*x2 + t2 = x2 * x2 for i in 0:MaxIter out += a out2 += b @@ -251,20 +273,19 @@ function besseli_power_series_inu_inum1(v, x::T) where T return out, out2 end - -# slightly modified version of https://github.com/heltonmc/Bessels.jl/issues/17#issuecomment-1195726642 from @cgeoga +# computes K_{nu+1}/K_{nu} using continued fractions and the modified Lentz method +# generally slow to converge for small x function besselk_ratio_knu_knup1(v, x::T) where T MaxIter = 1000 - # do the modified Lentz method: (hn, Dn, Cn) = (1e-50, zero(v), 1e-50) - jf = one(T) - vv = v*v + jf = one(T) + vv = v * v for _ in 1:MaxIter - an = (vv - ((2*jf - 1)^2) * T(0.25)) - bn = 2 * (x + jf) - Cn = an / Cn + bn - Dn = inv(muladd(an, Dn, bn)) + an = (vv - ((2*jf - 1)^2) * T(0.25)) + bn = 2 * (x + jf) + Cn = an / Cn + bn + Dn = inv(muladd(an, Dn, bn)) del = Dn * Cn hn *= del abs(del - 1) < eps(T) && break @@ -274,7 +295,6 @@ function besselk_ratio_knu_knup1(v, x::T) where T return xinv * (v + x + 1/2) + xinv * hn end - ##### ##### Power series for K_{nu}(x) ##### @@ -315,15 +335,16 @@ function besselk_power_series(v, x::T) where T _t2 = gam_nv * xd2_v * gam_1mnv (xd2_pow, fact_k, out) = (one(T), one(T), zero(T)) for k in 0:MaxIter - t1 = xd2_pow * T(0.5) - tmp = muladd(_t1, gam_1mnv, _t2 * gam_1mv) - tmp *= inv(gam_1mv * gam_1mnv * fact_k) - term = t1 * tmp - out += term - abs(term / out) < eps(T) && return out - (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T) + v + k), gam_1mv*(one(T) - v + k)) - xd2_pow *= zz - fact_k *= k + one(T) + t1 = xd2_pow * T(0.5) + tmp = muladd(_t1, gam_1mnv, _t2 * gam_1mv) + tmp *= inv(gam_1mv * gam_1mnv * fact_k) + term = t1 * tmp + out += term + abs(term / out) < eps(T) && return out + (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T) + v + k), gam_1mv*(one(T) - v + k)) + xd2_pow *= zz + fact_k *= k + one(T) end return out end +besselk_power_series_cutoff(nu, x) = x < 2.0 || nu > 1.6x - 1.0 diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 1638cf8..f2af6d9 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -68,8 +68,10 @@ k1x_32 = besselk1x.(Float32.(x)) @test besselk(100, 234.0) ≈ SpecialFunctions.besselk(100, 234.0) # test small arguments and order -m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:700.0] -@test [besselk(m, x) for m in m, x in x] ≈ [SpecialFunctions.besselk(m, x) for m in m, x in x] +m = 0:40; x = [1e-6; 1e-4; 1e-3; 1e-2; 0.1; 1.0:2.0:500.0] +for m in m, x in x + @test besselk(m, x) ≈ SpecialFunctions.besselk(m, x) +end # test medium arguments and order m = 30:200; x = 5.0:5.0:100.0 @@ -97,8 +99,19 @@ t = [besselk(m, x) for m in m, x in x] #@test isinf(besselk(250, 5.0)) ### Tests for besselkx -@test besselkx(0, 12.0) == besselk0x(12.0) -@test besselkx(1, 89.0) == besselk1x(89.0) +@test besselkx(0, 12.0) ≈ besselk0x(12.0) +@test besselkx(1, 89.0) ≈ besselk1x(89.0) @test besselkx(15, 82.123) ≈ SpecialFunctions.besselk(15, 82.123)*exp(82.123) @test besselkx(105, 182.123) ≈ SpecialFunctions.besselk(105, 182.123)*exp(182.123) + +## Tests for besselk + +## test all numbers and orders for 0 Date: Wed, 3 Aug 2022 18:03:09 -0400 Subject: [PATCH 09/14] add tests, fix F32 --- src/besseli.jl | 156 ++++++++++++++++++++++++++++--------------- src/gamma.jl | 2 + test/besseli_test.jl | 12 ++++ 3 files changed, 116 insertions(+), 54 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index f617764..9bc174d 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -137,17 +137,19 @@ Modified Bessel function of the first kind of order nu, ``I_{nu}(x)``. Nu must be real. """ function besseli(nu, x::T) where T <: Union{Float32, Float64} - nu == 0 && return besseli0(x) - nu == 1 && return besseli1(x) + iszero(nu) && return besseli0(x) + isone(nu) && return besseli1(x) - if x > maximum((T(30), nu^2 / 6)) - return T(besseli_large_argument(nu, x)) - elseif nu > 25.0 || x > 35.0 - return T(besseli_large_orders(nu, x)) - else - return T(besseli_power_series(nu, x)) - end + # use large argument expansion if x >> nu + besseli_large_argument_cutoff(nu, x) && return besseli_large_argument(nu, x) + + # use uniform debye expansion if x or nu is large + besselik_debye_cutoff(nu, x) && return T(besseli_large_orders(nu, x)) + + # for rest of values use the power series + return besseli_power_series(nu, x) end + """ besselix(nu, x::T) where T <: Union{Float32, Float64} @@ -155,19 +157,31 @@ Scaled modified Bessel function of the first kind of order nu, ``I_{nu}(x)*e^{-x Nu must be real. """ function besselix(nu, x::T) where T <: Union{Float32, Float64} - nu == 0 && return besseli0x(x) - nu == 1 && return besseli1x(x) - - if x > maximum((T(30), nu^2 / 4)) - return T(besseli_large_argument_scaled(nu, x)) - elseif x <= 2 * sqrt(nu + 1) - return T(besseli_power_series(nu, x)) * exp(-x) - elseif nu < 100 - return T(_besseli_continued_fractions_scaled(nu, x)) - else - return besseli_large_orders_scaled(nu, x) - end + iszero(nu) && return besseli0x(x) + isone(nu) && return besseli1x(x) + + # use large argument expansion if x >> nu + besseli_large_argument_cutoff(nu, x) && return besseli_large_argument_scaled(nu, x) + + # use uniform debye expansion if x or nu is large + besselik_debye_cutoff(nu, x) && return T(besseli_large_orders_scaled(nu, x)) + + # for rest of values use the power series + return besseli_power_series(nu, x) * exp(-x) end + +##### +##### Debye's uniform asymptotic for I_{nu}(x) +##### + +# Implements the uniform asymptotic expansion https://dlmf.nist.gov/10.41 +# In general this is valid when either x or nu is gets large +# see the file src/U_polynomials.jl for more details +""" + besseli_large_orders(nu, x::T) + +Debey's uniform asymptotic expansion for large order valid when v-> ∞ or x -> ∞ +""" function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64} S = promote_type(T, Float64) x = S(x) @@ -180,6 +194,7 @@ function besseli_large_orders(v, x::T) where T <: Union{Float32, Float64} return coef*Uk_poly_In(p, v, p2, T) end + function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64} S = promote_type(T, Float64) x = S(x) @@ -192,39 +207,18 @@ function besseli_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64} return T(coef*Uk_poly_In(p, v, p2, T)) end -function _besseli_continued_fractions(nu, x::T) where T - S = promote_type(T, Float64) - xx = S(x) - knum1, knu = besselk_up_recurrence(xx, besselk1(xx), besselk0(xx), 1, nu-1) - # if knu or knum1 is zero then besseli will likely overflow - (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error")) - return 1 / (x * (knum1 + knu / steed(nu, x))) -end -function _besseli_continued_fractions_scaled(nu, x::T) where T - S = promote_type(T, Float64) - xx = S(x) - knum1, knu = besselk_up_recurrence(xx, besselk1x(xx), besselk0x(xx), 1, nu-1) - # if knu or knum1 is zero then besseli will likely overflow - (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error")) - return 1 / (x * (knum1 + knu / steed(nu, x))) -end -function steed(n, x::T) where T - MaxIter = 1000 - xinv = inv(x) - xinv2 = 2 * xinv - d = x / (n + n) - a = d - h = a - b = muladd(2, n, 2) * xinv - for _ in 1:MaxIter - d = inv(b + d) - a *= muladd(b, d, -1) - h = h + a - b = b + xinv2 - abs(a / h) <= eps(T) && break - end - return h -end + +##### +##### Large argument expansion (x>>nu) for I_{nu}(x) +##### + +# Implements the uniform asymptotic expansion https://dlmf.nist.gov/10.40.E1 +# In general this is valid when x > nu^2 +""" + besseli_large_orders(nu, x::T) + +Debey's uniform asymptotic expansion for large order valid when v-> ∞ or x -> ∞ +""" function besseli_large_argument(v, z::T) where T MaxIter = 1000 a = exp(z / 2) @@ -260,7 +254,19 @@ function besseli_large_argument_scaled(v, z::T) where T end return res * coef end +besseli_large_argument_cutoff(nu, x) = x > maximum((30.0, nu^2 / 6)) + +##### +##### Power series for I_{nu}(x) +##### + +# Use power series form of I_v(x) which is generally accurate across all values though slower for larger x +# https://dlmf.nist.gov/10.25.E2 +""" + besseli_power_series(nu, x::T) where T <: Float64 +Computes ``I_{nu}(x)`` using the power series for any value of nu. +""" function besseli_power_series(v, x::T) where T MaxIter = 3000 out = zero(T) @@ -274,3 +280,45 @@ function besseli_power_series(v, x::T) where T end return out end + +#= +# the following is a deprecated version of the continued fraction approach +# using K0 and K1 as starting values then forward recurrence up till nu +# then using the wronskian to getting I_{nu} +# in general this method is slow and depends on starting values of K0 and K1 +# which is not very flexible for arbitary orders + +function _besseli_continued_fractions(nu, x::T) where T + S = promote_type(T, Float64) + xx = S(x) + knum1, knu = besselk_up_recurrence(xx, besselk1(xx), besselk0(xx), 1, nu-1) + # if knu or knum1 is zero then besseli will likely overflow + (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error")) + return 1 / (x * (knum1 + knu / steed(nu, x))) +end +function _besseli_continued_fractions_scaled(nu, x::T) where T + S = promote_type(T, Float64) + xx = S(x) + knum1, knu = besselk_up_recurrence(xx, besselk1x(xx), besselk0x(xx), 1, nu-1) + # if knu or knum1 is zero then besseli will likely overflow + (iszero(knu) || iszero(knum1)) && return throw(DomainError(x, "Overflow error")) + return 1 / (x * (knum1 + knu / steed(nu, x))) +end +function steed(n, x::T) where T + MaxIter = 1000 + xinv = inv(x) + xinv2 = 2 * xinv + d = x / (n + n) + a = d + h = a + b = muladd(2, n, 2) * xinv + for _ in 1:MaxIter + d = inv(b + d) + a *= muladd(b, d, -1) + h = h + a + b = b + xinv2 + abs(a / h) <= eps(T) && break + end + return h +end +=# diff --git a/src/gamma.jl b/src/gamma.jl index 56283a1..488ad00 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -8,6 +8,8 @@ function gamma(x) return _gamma(x) end end +# only have a Float64 implementations +gamma(x::Float32) = Float32(gamma(Float64(x))) function _gamma(x) if x > 11.5 return large_gamma(x) diff --git a/test/besseli_test.jl b/test/besseli_test.jl index 6a5bcb6..e4398c8 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -90,4 +90,16 @@ t = [besseli(m, x) for m in m, x in x] t = [besselix(m, x) for m in m, x in x] @test t[10] isa Float64 @test t ≈ [SpecialFunctions.besselix(m, x) for m in m, x in x] + +## Tests for besselk + +## test all numbers and orders for 0 Date: Wed, 3 Aug 2022 18:27:11 -0400 Subject: [PATCH 10/14] add docs --- src/besseli.jl | 29 +++++++++++++++++++++++++++++ src/besselk.jl | 36 ++++++++++++++++++++++++++++++++++-- 2 files changed, 63 insertions(+), 2 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 9bc174d..e988378 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -130,6 +130,35 @@ function besseli1x(x::T) where T <: Union{Float32, Float64} return z end +# Modified Bessel functions of the first kind of order nu +# besseli(nu, x) +# +# A numerical routine to compute the modified Bessel function of the first kind I_{ν}(x) [1] +# for real orders and arguments of positive or negative value. The routine is based on several +# publications [2-6] that calculate I_{ν}(x) for positive arguments and orders where +# reflection identities are used to compute negative arguments and orders. +# +# In particular, the reflectance identities for negative noninteger orders I_{−ν}(x) = I_{ν}(x) + 2 / πsin(πν)*Kν(x) +# and for negative integer orders I_{−n}(x) = I_n(x) are used. +# For negative arguments of integer order, In(−x) = (−1)^n In(x) is used and for +# noninteger orders, Iν(−x) = exp(iπν) Iν(x) is used. For negative orders and arguments the previous identities are combined. +# +# The identities are computed by calling the `besseli_positive_args(nu, x)` function which computes I_{ν}(x) +# for positive arguments and orders. For large orders, Debye's uniform asymptotic expansions are used where large arguments (x>>nu) +# a large argument expansion is used. The rest of the values are computed using the power series. + +# [1] https://dlmf.nist.gov/10.40.E1 +# [2] Amos, Donald E. "Computation of modified Bessel functions and their ratios." Mathematics of computation 28.125 (1974): 239-251. +# [3] Gatto, M. A., and J. B. Seery. "Numerical evaluation of the modified Bessel functions I and K." +# Computers & Mathematics with Applications 7.3 (1981): 203-209. +# [4] Temme, Nico M. "On the numerical evaluation of the modified Bessel function of the third kind." +# Journal of Computational Physics 19.3 (1975): 324-337. +# [5] Amos, DEv. "Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order." +# ACM Transactions on Mathematical Software (TOMS) 12.3 (1986): 265-273. +# [6] Segura, Javier, P. Fernández de Córdoba, and Yu L. Ratis. "A code to evaluate modified bessel functions based on thecontinued fraction method." +# Computer physics communications 105.2-3 (1997): 263-272. +# + """ besseli(nu, x::T) where T <: Union{Float32, Float64} diff --git a/src/besselk.jl b/src/besselk.jl index 11e15cd..8c51adc 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -146,6 +146,40 @@ function besselk1x(x::T) where T <: Union{Float32, Float64} end end +# Modified Bessel functions of the second kind of order nu +# besselk(nu, x) +# +# A numerical routine to compute the modified Bessel function of the second kind K_{ν}(x) [1] +# for real orders and arguments of positive or negative value. The routine is based on several +# publications [2-8] that calculate K_{ν}(x) for positive arguments and orders where +# reflection identities are used to compute negative arguments and orders. +# +# In particular, the reflectance identities for negative orders I_{−ν}(x) = I_{ν}(x). +# For negative arguments of integer order, Kn(−x) = (−1)^n Kn(x) − im * π In(x) is used and for +# noninteger orders, Kν(−x) = exp(−iπν)*Kν(x) − im π Iν(x) is used. For negative orders and arguments the previous identities are combined. +# +# The identities are computed by calling the `besseli_positive_args(nu, x)` function which computes K_{ν}(x) +# for positive arguments and orders. For large orders, Debye's uniform asymptotic expansions are used. +# For small value and when nu > ~x the power series is used. The rest of the values are computed using slightly different methods. +# The power series for besseli is modified to give both I_{v} and I_{v-1} where the ratio K_{v+1} / K_{v} is computed using continued fractions [8]. +# The wronskian connection formula is then used to compute K_v. + +# [1] http://dlmf.nist.gov/10.27.E4 +# [2] Amos, Donald E. "Computation of modified Bessel functions and their ratios." Mathematics of computation 28.125 (1974): 239-251. +# [3] Gatto, M. A., and J. B. Seery. "Numerical evaluation of the modified Bessel functions I and K." +# Computers & Mathematics with Applications 7.3 (1981): 203-209. +# [4] Temme, Nico M. "On the numerical evaluation of the modified Bessel function of the third kind." +# Journal of Computational Physics 19.3 (1975): 324-337. +# [5] Amos, DEv. "Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order." +# ACM Transactions on Mathematical Software (TOMS) 12.3 (1986): 265-273. +# [6] Segura, Javier, P. Fernández de Córdoba, and Yu L. Ratis. "A code to evaluate modified bessel functions based on thecontinued fraction method." +# Computer physics communications 105.2-3 (1997): 263-272. +# [7] Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." +# arXiv preprint arXiv:2201.00090 (2022). +# [8] Cuyt, A. A., Petersen, V., Verdonk, B., Waadeland, H., & Jones, W. B. (2008). +# Handbook of continued fractions for special functions. Springer Science & Business Media. +# + """ besselk(x::T) where T <: Union{Float32, Float64} @@ -305,8 +339,6 @@ end # A modified form appears below. See more discussion at https://github.com/heltonmc/Bessels.jl/pull/29 # This is only accurate for noninteger orders (nu) and no checks are performed. # -# [7] Geoga, Christopher J., et al. "Fitting Mat\'ern Smoothness Parameters Using Automatic Differentiation." -# arXiv preprint arXiv:2201.00090 (2022). """ besselk_power_series(nu, x::T) where T <: Float64 From 7278cfef29c028869066f2722b354f952da6ab98 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 19:09:48 -0400 Subject: [PATCH 11/14] add support for negative args --- src/U_polynomials.jl | 10 ---------- src/besseli.jl | 45 ++++++++++++++++++++++++++++++++++++++++---- src/besselk.jl | 40 ++++++++++++++++++++++++++++++++++----- test/besseli_test.jl | 20 +++++++++++++++++++- test/besselk_test.jl | 20 ++++++++++++++++++++ 5 files changed, 115 insertions(+), 20 deletions(-) diff --git a/src/U_polynomials.jl b/src/U_polynomials.jl index b9cf29c..b6cd84d 100644 --- a/src/U_polynomials.jl +++ b/src/U_polynomials.jl @@ -97,16 +97,6 @@ Uk_poly_Kn(p, v, p2, ::Type{Float64}) = Uk_poly10(p, v, p2)[2] end end -function Uk_poly3(p, v, p2) - u0 = 1.0 - u1 = evalpoly(p2, (0.125, -0.20833333333333334)) - u2 = evalpoly(p2, (0.0703125, -0.4010416666666667, 0.3342013888888889)) - u3 = evalpoly(p2, (0.0732421875, -0.8912109375, 1.8464626736111112, -1.0258125964506173)) - - Poly = (u0, u1, u2, u3) - - return split_evalpoly(-p/v, Poly) -end function Uk_poly5(p, v, p2) u0 = 1.0 u1 = evalpoly(p2, (0.125, -0.20833333333333334)) diff --git a/src/besseli.jl b/src/besseli.jl index e988378..72ac052 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -160,12 +160,49 @@ end # """ - besseli(nu, x::T) where T <: Union{Float32, Float64} + besselk(x::T) where T <: Union{Float32, Float64} -Modified Bessel function of the first kind of order nu, ``I_{nu}(x)``. -Nu must be real. +Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. +""" +function besseli(nu::Real, x::T) where T + isinteger(nu) && return besseli(Int(nu), x) + abs_nu = abs(nu) + abs_x = abs(x) + + if nu >= 0 + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + else + return cispi(abs_nu) * besseli_positive_args(abs_nu, abs_x) + end + else + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + 2 / π * sinpi(abs_nu) * besselk_positive_args(abs_nu, abs_x) + else + Iv = besseli_positive_args(abs_nu, abs_x) + Kv = besselk_positive_args(abs_nu, abs_x) + return cispi(abs_nu) * Iv + 2 / π * sinpi(abs_nu) * (cispi(-abs_nu) * Kv - im * π * Iv) + end + end +end +function besseli(nu::Integer, x::T) where T + abs_nu = abs(nu) + abs_x = abs(x) + sg = iseven(abs_nu) ? 1 : -1 + + if x >= 0 + return besseli_positive_args(abs_nu, abs_x) + else + return sg * besseli_positive_args(abs_nu, abs_x) + end +end + +""" + besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64} + +Modified Bessel function of the first kind of order nu, ``I_{nu}(x)`` for positive arguments. """ -function besseli(nu, x::T) where T <: Union{Float32, Float64} +function besseli_positive_args(nu, x::T) where T <: Union{Float32, Float64} iszero(nu) && return besseli0(x) isone(nu) && return besseli1(x) diff --git a/src/besselk.jl b/src/besselk.jl index 8c51adc..0fed8ab 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -185,8 +185,38 @@ end Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. """ -function besselk(nu, x::T) where T <: Union{Float32, Float64} - # check to make sure nu isn't zero +function besselk(nu::Real, x::T) where T + isinteger(nu) && return besselk(Int(nu), x) + abs_nu = abs(nu) + abs_x = abs(x) + + if x >= 0 + return besselk_positive_args(abs_nu, abs_x) + else + return cispi(-abs_nu)*besselk_positive_args(abs_nu, abs_x) - besseli_positive_args(abs_nu, abs_x) * im * π + end +end +function besselk(nu::Integer, x::T) where T + abs_nu = abs(nu) + abs_x = abs(x) + sg = iseven(abs_nu) ? 1 : -1 + + if x >= 0 + return besselk_positive_args(abs_nu, abs_x) + else + return sg * besselk_positive_args(abs_nu, abs_x) - im * π * besseli_positive_args(abs_nu, abs_x) + end +end + +""" + besselk_positive_args(x::T) where T <: Union{Float32, Float64} + +Modified Bessel function of the second kind of order nu, ``K_{nu}(x)`` valid for postive arguments and orders. +""" +function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} + iszero(x) && return T(Inf) + + # dispatch to avoid uniform expansion when nu = 0 iszero(nu) && return besselk0(x) # use uniform debye expansion if x or nu is large @@ -207,7 +237,7 @@ end Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x}``. """ function besselkx(nu, x::T) where T <: Union{Float32, Float64} - # check to make sure nu isn't zero + # dispatch to avoid uniform expansion when nu = 0 iszero(nu) && return besselk0x(x) # use uniform debye expansion if x or nu is large @@ -333,11 +363,11 @@ end ##### Power series for K_{nu}(x) ##### -# Use power series form of K_v(x) which is accurate for small x (x<2) or when nu > X +# Use power series form of K_v(x) which is accurate for small x (x<2) or when nu > x # We use the form as described by Equation 3.2 from reference [7]. # This method was originally contributed by @cgeoga https://github.com/cgeoga/BesselK.jl/blob/main/src/besk_ser.jl # A modified form appears below. See more discussion at https://github.com/heltonmc/Bessels.jl/pull/29 -# This is only accurate for noninteger orders (nu) and no checks are performed. +# This is only valid for noninteger orders (nu) and no checks are performed. # """ besselk_power_series(nu, x::T) where T <: Float64 diff --git a/test/besseli_test.jl b/test/besseli_test.jl index e4398c8..ab8ddeb 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -102,4 +102,22 @@ for v in nu, xx in x @test isapprox(besseli(v, xx), Float64(sf), rtol=2e-13) end - \ No newline at end of file + ### tests for negative arguments + +(v, x) = 12.0, 3.2 +@test besseli(v,x) ≈ 7.1455266650203694069897133431E-7 + +(v, x) = -8.0, 4.2 +@test besseli(v,x) ≈ 0.0151395115677545706602449919627 + +(v, x) = 12.3, 8.2 +@test besseli(v,x) ≈ 0.113040018422133018059759059298 + +(v, x) = -12.3, 8.2 +@test besseli(v,x) ≈ 0.267079696793126091886043602895 + +(v, x) = -14.0, -9.9 +@test besseli(v,x) ≈ 0.2892290867115615816280234648 + +(v, x) = -14.6, -10.6 +@test besseli(v,x) ≈ -0.157582642056898598750175404443 - 0.484989503203097528858271270828*im diff --git a/test/besselk_test.jl b/test/besselk_test.jl index f2af6d9..ff3040f 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -115,3 +115,23 @@ for v in nu, xx in x sf = SpecialFunctions.besselk(v, xx) @test isapprox(besselk(v, xx), Float64(sf), rtol=2e-13) end + +### tests for negative arguments + +(v, x) = 12.0, 3.2 +@test besselk(v,x) ≈ 56331.504348755621996013084096 + +(v, x) = -8.0, 4.2 +@test besselk(v,x) ≈ 3.65165949039881135495282699061 + +(v, x) = 12.3, 8.2 +@test besselk(v,x) ≈ 0.299085139926840649406079315812 + +(v, x) = -12.3, 8.2 +@test besselk(v,x) ≈ 0.299085139926840649406079315812 + +(v, x) = -14.0, -9.9 +@test besselk(v,x) ≈ 0.100786833375605803570325345603 - 0.90863997401752715470886289641*im + +(v, x) = -14.6, -10.6 +@test besselk(v,x) ≈ -0.0180385087581148387140033906859 - 1.54653251445680014758965158559*im From aea768b96a0731c5c04f100a45ef95e39b50f094 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 19:25:21 -0400 Subject: [PATCH 12/14] fix coverage --- src/besselk.jl | 2 +- src/recurrence.jl | 4 ++-- test/besselk_test.jl | 4 ++++ 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/besselk.jl b/src/besselk.jl index 0fed8ab..52be2c8 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -215,7 +215,7 @@ Modified Bessel function of the second kind of order nu, ``K_{nu}(x)`` valid for """ function besselk_positive_args(nu, x::T) where T <: Union{Float32, Float64} iszero(x) && return T(Inf) - + # dispatch to avoid uniform expansion when nu = 0 iszero(nu) && return besselk0(x) diff --git a/src/recurrence.jl b/src/recurrence.jl index f32c611..fd808ca 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -36,7 +36,7 @@ end return jnup1, jnu end - +#= # currently not used # backward recurrence relation for besselk and besseli # outputs both (bessel(x, nu_end), bessel(x, nu_end-1) @@ -50,4 +50,4 @@ end end return jnup1, jnu end - +=# diff --git a/test/besselk_test.jl b/test/besselk_test.jl index ff3040f..5788947 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -104,6 +104,10 @@ t = [besselk(m, x) for m in m, x in x] @test besselkx(15, 82.123) ≈ SpecialFunctions.besselk(15, 82.123)*exp(82.123) @test besselkx(105, 182.123) ≈ SpecialFunctions.besselk(105, 182.123)*exp(182.123) +@test besselkx(4, 3.0) ≈ SpecialFunctions.besselk(4, 3.0)*exp(3.0) +@test besselkx(3.2, 1.1) ≈ SpecialFunctions.besselk(3.2, 1.1)*exp(1.1) +@test besselkx(8.2, 9.1) ≈ SpecialFunctions.besselk(8.2, 9.1)*exp(9.1) + ## Tests for besselk From ab6a5e5161a3ebe67b26ad6e34a15ba0b4cd22d5 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 19:38:17 -0400 Subject: [PATCH 13/14] add more negative arg tests --- test/besseli_test.jl | 8 +++++++- test/besselk_test.jl | 6 ++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/test/besseli_test.jl b/test/besseli_test.jl index ab8ddeb..d25f3dd 100644 --- a/test/besseli_test.jl +++ b/test/besseli_test.jl @@ -105,7 +105,13 @@ end ### tests for negative arguments (v, x) = 12.0, 3.2 -@test besseli(v,x) ≈ 7.1455266650203694069897133431E-7 +@test besseli(v,x) ≈ 7.1455266650203694069897133431e-7 + +(v,x) = 13.0, -1.0 +@test besseli(v,x) ≈ -1.995631678207200756444e-14 + +(v,x) = 12.6, -3.0 +@test besseli(v,x) ≈ -2.725684975265100582482e-8 + 8.388795775899337839603e-8 * im (v, x) = -8.0, 4.2 @test besseli(v,x) ≈ 0.0151395115677545706602449919627 diff --git a/test/besselk_test.jl b/test/besselk_test.jl index 5788947..a908bad 100644 --- a/test/besselk_test.jl +++ b/test/besselk_test.jl @@ -125,6 +125,12 @@ end (v, x) = 12.0, 3.2 @test besselk(v,x) ≈ 56331.504348755621996013084096 +(v,x) = 13.0, -1.0 +@test besselk(v,x) ≈ -1921576392792.994084565 - 6.26946181952681217841e-14*im + +(v,x) = 12.6, -3.0 +@test besselk(v,x) ≈ -135222.7926354826692727 - 416172.9627453652120223*im + (v, x) = -8.0, 4.2 @test besselk(v,x) ≈ 3.65165949039881135495282699061 From dfbd16db42a754230034d7dc263a3ca8a06f2bf4 Mon Sep 17 00:00:00 2001 From: Michael Helton Date: Wed, 3 Aug 2022 19:43:55 -0400 Subject: [PATCH 14/14] fix doc string error --- src/besseli.jl | 4 ++-- src/besselk.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/besseli.jl b/src/besseli.jl index 72ac052..714d0e1 100644 --- a/src/besseli.jl +++ b/src/besseli.jl @@ -160,9 +160,9 @@ end # """ - besselk(x::T) where T <: Union{Float32, Float64} + besseli(x::T) where T <: Union{Float32, Float64} -Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``. +Modified Bessel function of the second kind of order nu, ``I_{nu}(x)``. """ function besseli(nu::Real, x::T) where T isinteger(nu) && return besseli(Int(nu), x) diff --git a/src/besselk.jl b/src/besselk.jl index 52be2c8..7a8caa5 100644 --- a/src/besselk.jl +++ b/src/besselk.jl @@ -402,7 +402,7 @@ function besselk_power_series(v, x::T) where T tmp *= inv(gam_1mv * gam_1mnv * fact_k) term = t1 * tmp out += term - abs(term / out) < eps(T) && return out + abs(term / out) < eps(T) && break (gam_1mnv, gam_1mv) = (gam_1mnv*(one(T) + v + k), gam_1mv*(one(T) - v + k)) xd2_pow *= zz fact_k *= k + one(T)