Merge pull request #7 from heltonmc/besselk

Faster implementations of Besselk0 and Besselk1

Author: heltonmc

src/Bessels.jl

@@ -12,6 +12,8 @@ export besseli0x
 export besseli1
 export besseli1x
 
+export besselk
+export besselkx
 export besselk0
 export besselk0x
 export besselk1
@@ -29,7 +31,6 @@ include("Float128/besselk.jl")
 include("Float128/bessely.jl")
 include("Float128/constants.jl")
 
-include("chebyshev.jl")
 include("math_constants.jl")
 #include("parse.jl")
 #include("hankel.jl")

src/besselk.jl

@@ -1,60 +1,175 @@
-#=
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 2000 by Stephen L. Moshier
-https://github.com/jeremybarnes/cephes/blob/master/bessel/k0.c
-https://github.com/jeremybarnes/cephes/blob/master/bessel/k1.c
-=#
-function besselk0(x::T) where T <: Union{Float32, Float64}
-    if x <= zero(x)
-        return throw(DomainError(x, "NaN result for non-NaN input."))
-    end
+#    Modified Bessel functions of the second kind of order zero and one
+#                       besselk0, besselk1
+#
+#    Scaled modified Bessel functions of the second kind of order zero and one
+#                       besselk0x, besselk1x
+#
+#    (Scaled) Modified Bessel functions of the second kind of order nu
+#                       besselk, besselkx
+#
+#
+#    Calculation of besselk0 is done in two branches using polynomial approximations [1]
+#
+#    Branch 1: x < 1.0 
+#              besselk0(x) + log(x)*besseli0(x) = P7(x^2)
+#                            besseli0(x) = [x/2]^2 * P6([x/2]^2) + 1
+#    Branch 2: x >= 1.0
+#              sqrt(x) * exp(x) * besselk0(x) = P22(1/x) / Q2(1/x)
+#    where P7, P6, P22, and Q2 are 7, 6, 22, and 2 degree polynomials respectively.
+#
+#
+#
+#    Calculation of besselk1 is done in two branches using polynomial approximations [2]
+#
+#    Branch 1: x < 1.0 
+#              besselk1(x) - log(x)*besseli1(x) - 1/x = x*P8(x^2)
+#                            besseli1(x) = [x/2]^2 * (1 + 0.5 * (*x/2)^2 + (x/2)^4 * P5([x/2]^2))
+#    Branch 2: x >= 1.0
+#              sqrt(x) * exp(x) * besselk1(x) = P22(1/x) / Q2(1/x)
+#    where P8, P5, P22, and Q2 are 8, 5, 22, and 2 degree polynomials respectively.
+#
+#
+#    The polynomial coefficients are taken from boost math library [3].
+#    Evaluation of the coefficients using Remez.jl is prohibitive due to the increase
+#    precision required in ArbNumerics.jl. 
+#
+#    Horner's scheme is then used to evaluate all polynomials.
+#
+#    Calculation of besselk and besselkx can be done with recursion starting with
+#    besselk0 and besselk1 and using upward recursion for any value of nu (order).
+#
+#                    K_{nu+1} = (2*nu/x)*K_{nu} + K_{nu-1}
+#
+#    When nu is large, a large amount of recurrence is necesary.
+#    At this point as nu -> Inf it is more efficient to use a uniform expansion.
+#    The boundary of the expansion depends on the accuracy required.
+#    For more information see [4]. This approach is not yet implemented, so recurrence
+#    is used for all values of nu.
+#
+# 
+# [1] "Rational Approximations for the Modified Bessel Function of the Second Kind 
+#     - K0(x) for Computations with Double Precision" by Pavel Holoborodko     
+# [2] "Rational Approximations for the Modified Bessel Function of the Second Kind 
+#     - K1(x) for Computations with Double Precision" by Pavel Holoborodko
+# [3] https://github.com/boostorg/math/tree/develop/include/boost/math/special_functions/detail
+# [4] "Computation of Bessel Functions of Complex Argument and Large ORder" by Donald E. Amos
+#      Sandia National Laboratories
+#
+"""
+    besselk0(x::T) where T <: Union{Float32, Float64}
 
-    if x <= 2
-        y = muladd(x, x, T(-2))
-        y = chbevl(y, A_k0(T)) - log(T(.5) * x) * besseli0(x)
-        return y
+Modified Bessel function of the second kind of order zero, ``K_0(x)``.
+"""
+function besselk0(x::T) where T <: Union{Float32, Float64}
+    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
+    if x <= one(T)
+        a = x * x / 4
+        s = muladd(evalpoly(a, P1_k0(T)), inv(evalpoly(a, Q1_k0(T))), T(Y_k0))
+        a = muladd(s, a, 1)
+        return muladd(-a, log(x), evalpoly(x * x, P2_k0(T)))
     else
-        z = T(8) / x - T(2)
-        return exp(-x) * chbevl(z, B_k0(T)) / sqrt(x)
+        s = exp(-x / 2)
+        a = muladd(evalpoly(inv(x), P3_k0(T)), inv(evalpoly(inv(x), Q3_k0(T))), one(T)) * s / sqrt(x)
+        return a * s
     end
 end
+"""
+    besselk0x(x::T) where T <: Union{Float32, Float64}
 
+Scaled modified Bessel function of the second kind of order zero, ``K_0(x)*e^{x}``.
+"""
 function besselk0x(x::T) where T <: Union{Float32, Float64}
-    if x <= zero(x)
-        return throw(DomainError(x, "NaN result for non-NaN input."))
-    end
-    if x <= 2
-        y =  muladd(x, x, T(-2))
-        y = chbevl(y, A_k0(T)) - log(T(.5) * x) * besseli0(x)
-        return y * exp(x)
+    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
+    if x <= one(T)
+        a = x * x / 4
+        s = muladd(evalpoly(a, P1_k0(T)), inv(evalpoly(a, Q1_k0(T))), T(Y_k0))
+        a = muladd(s, a, 1)
+        return muladd(-a, log(x), evalpoly(x * x, P2_k0(T))) * exp(x)
     else
-        z = T(8) / x - T(2)
-        return chbevl(z, B_k0(T)) / sqrt(x)
+        return muladd(evalpoly(inv(x), P3_k0(T)), inv(evalpoly(inv(x), Q3_k0(T))), one(T)) / sqrt(x)
     end
 end
+"""
+    besselk1(x::T) where T <: Union{Float32, Float64}
+
+Modified Bessel function of the second kind of order one, ``K_1(x)``.
+"""
 function besselk1(x::T) where T <: Union{Float32, Float64}
-    z = T(.5) * x
-    if x <= zero(x)
-        return throw(DomainError(x, "NaN result for non-NaN input."))
-    end
-    if x <= 2
-        y = muladd(x, x, T(-2))
-        y = log(z) * besseli1(x) + chbevl(y, A_k1(T)) / x
-        return y
+    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
+    if x <= one(T)
+        z = x * x
+        a = z / 4
+        pq = muladd(evalpoly(a, P1_k1(T)), inv(evalpoly(a, Q1_k1(T))), Y_k1(T))
+        pq = muladd(pq * a, a, (a / 2 + 1))
+        a = pq * x / 2
+        pq = muladd(evalpoly(z, P2_k1(T)) / evalpoly(z, Q2_k1(T)), x, inv(x))
+        return muladd(a, log(x), pq)
     else
-        return exp(-x) * chbevl(T(8) / x - T(2), B_k1(T)) / sqrt(x)
+        s = exp(-x / 2)
+        a = muladd(evalpoly(inv(x), P3_k1(T)), inv(evalpoly(inv(x), Q3_k1(T))), Y2_k1(T)) * s / sqrt(x)
+        return a * s
     end
 end
+"""
+    besselk1x(x::T) where T <: Union{Float32, Float64}
+
+Scaled modified Bessel function of the second kind of order one, ``K_1(x)*e^{x}``.
+"""
 function besselk1x(x::T) where T <: Union{Float32, Float64}
-    z = T(.5) * x
-    if x <= zero(x)
-        return throw(DomainError(x, "NaN result for non-NaN input."))
-    end
-    if x <= 2
-        y = muladd(x, x, T(-2))
-        y = log(z) * besseli1(x) + chbevl(y, A_k1(T)) / x
-        return y * exp(x)
+    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
+    if x <= one(T)
+        z = x * x
+        a = z / 4
+        pq = muladd(evalpoly(a, P1_k1(T)), inv(evalpoly(a, Q1_k1(T))), Y_k1(T))
+        pq = muladd(pq * a, a, (a / 2 + 1))
+        a = pq * x / 2
+        pq = muladd(evalpoly(z, P2_k1(T)) / evalpoly(z, Q2_k1(T)), x, inv(x))
+        return muladd(a, log(x), pq) * exp(x)
     else
-        return chbevl(T(8) / x - T(2), B_k1(T)) / sqrt(x)
+        return muladd(evalpoly(inv(x), P3_k1(T)), inv(evalpoly(inv(x), Q3_k1(T))), Y2_k1(T)) / sqrt(x)
+    end
+end
+#=
+Recurrence is used for all values of nu. 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.
+
+In the future, a more efficient algorithm for large nu should be incorporated.
+=#
+"""
+    besselk(x::T) where T <: Union{Float32, Float64}
+
+Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
+"""
+function besselk(nu::Int, x::T) where T <: Union{Float32, Float64}
+    return three_term_recurrence(x, besselk0(x), besselk1(x), nu)
+end
+"""
+    besselk(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}
+    return three_term_recurrence(x, besselk0x(x), besselk1x(x), nu)
+end
+
+@inline function three_term_recurrence(x, k0, k1, nu)
+    nu == 0 && return k0
+    nu == 1 && return k1
+
+    # this prevents us from looping through large values of nu when the loop will always return zero
+    (iszero(k0) || iszero(k1)) && return zero(x) 
+
+    k2 = k0
+    x2 = 2 / x
+    for n in 1:nu-1
+        a = x2 * n
+        k2 = muladd(a, k1, k0)
+        k0 = k1
+        k1 = k2
     end
+    return k2
 end