Merge pull request #26 from heltonmc/bessely

bessely(nu, x) for positive and negative arguments

Author: heltonmc

.github/workflows/CI.yml

@@ -20,7 +20,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1'
           - '1.6'
           - 'nightly'
         os:

src/Bessels.jl

@@ -4,9 +4,11 @@ using SpecialFunctions: loggamma
 
 export besselj0
 export besselj1
+export besselj
 
 export bessely0
 export bessely1
+export bessely
 
 export besseli
 export besselix
@@ -40,6 +42,8 @@ include("recurrence.jl")
 include("misc.jl")
 include("Polynomials/besselj_polys.jl")
 include("asymptotics.jl")
+include("gamma.jl")
+
 #include("hankel.jl")
 
 end

src/U_polynomials.jl

@@ -1,3 +1,12 @@
+#besseljy_large_argument_min(::Type{Float32}) = 15.0f0
+besseljy_large_argument_min(::Type{Float64}) = 20.0
+besseljy_large_argument_min(::Type{T}) where T <: AbstractFloat = 40.0
+
+#besseljy_large_argument_cutoff(v, x::Float32) = (x > 1.2f0*v && x > besseljy_large_argument_min(Float32))
+besseljy_large_argument_cutoff(v, x::Float64) = (x > 1.65*v && x > besseljy_large_argument_min(Float64))
+besseljy_large_argument_cutoff(v, x::T) where T = (x > 4*v && x > besseljy_large_argument_min(T))
+
+
 function besseljy_large_argument(v, x::T) where T
     # gives both (besselj, bessely) for x > 1.6*v
     α, αp = _α_αp_asymptotic(v, x)
@@ -30,6 +39,7 @@ function _α_αp_asymptotic(v, x::Float64)
         return _α_αp_poly_30(v, x)
     end
 end
+#=
 function _α_αp_asymptotic(v, x::Float32)
     if x > 4*v
         return _α_αp_poly_5(v, x)
@@ -39,7 +49,7 @@ function _α_αp_asymptotic(v, x::Float32)
         return _α_αp_poly_30(v, x)
     end
 end
-
+=#
 # Float64
 # can only use for x > 20.0
 # 30 terms gives ~5e-16 relative error when x > 1.6*nu
@@ -74,6 +84,7 @@ end
 # a = 27.7479; b = 22.3588; c = 3.74567
 # this method requires significantly more terms when x gets closer to nu
 # so it becomes more efficient to use recurrence (or another algorithm) in this region
+#=
 function _α_αp_poly_5(v, x::T) where T
     xinv = inv(x)^2
     μ = 4 * T(v)^2
@@ -89,6 +100,7 @@ function _α_αp_poly_5(v, x::T) where T
     return α, αp
     return α, αp
 end
+=#
 function _α_αp_poly_10(v, x::T) where T
     xinv = inv(x)^2
     μ = 4 * v^2
@@ -109,6 +121,7 @@ function _α_αp_poly_10(v, x::T) where T
     return α, αp
     return α, αp
 end
+#=
 function _α_αp_poly_15(v, x::T) where T
     xinv = inv(x)^2
     μ = 4 * v^2
@@ -133,6 +146,7 @@ function _α_αp_poly_15(v, x::T) where T
     α = x * (evalpoly(xinv, (s0, -s1, -s2/3, -s3/5, -s4/7, -s5/9, -s6/11, -s7/13, -s8/15, -s9/17, -s10/19, -s11/21, -s12/23, -s13/25, -s14/27, -s15/29)))
     return α, αp
 end
+=#
 function _α_αp_poly_20(v, x::T) where T
     xinv = inv(x)^2
     μ = 4 * v^2

src/asymptotics.jl

@@ -198,15 +198,15 @@ end
 function _besseli_continued_fractions(nu, x::T) where T
     S = promote_type(T, Float64)
     xx = S(x)
-    knu, knum1 = up_recurrence(xx, besselk0(xx), besselk1(xx), nu)
+    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)
-    knu, knum1 = up_recurrence(xx, besselk0x(xx), besselk1x(xx), nu)
+    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)))

src/besseli.jl

@@ -150,13 +150,43 @@ function besselj1(x::Float32)
         return p * s
     end
 end
-"""
-    besselj(nu, x::T) where T <: Union{Float32, Float64}
 
-Bessel function of the first kind of order nu, ``J_{nu}(x)``.
-Nu must be real.
-"""
-function _besselj(nu, x)
+function besselj(nu::Real, x::T) where T
+    isinteger(nu) && return besselj(Int(nu), x)
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+
+    Jnu = besselj_positive_args(abs_nu, abs_x)
+    if nu >= zero(T)
+        return x >= zero(T) ? Jnu : Jnu * cispi(abs_nu)
+    else
+        Ynu = bessely_positive_args(abs_nu, abs_x)
+        spi, cpi = sincospi(abs_nu)
+        out = Jnu * cpi - Ynu * spi
+        return x >= zero(T) ? out : out * cispi(nu)
+    end
+end
+
+function besselj(nu::Integer, x::T) where T
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+    sg = iseven(abs_nu) ? 1 : -1
+
+    Jnu = besselj_positive_args(abs_nu, abs_x)
+    if nu >= zero(T)
+        return x >= zero(T) ? Jnu : Jnu * sg
+    else
+        if x >= zero(T)
+            return Jnu * sg
+        else
+            Ynu = bessely_positive_args(abs_nu, abs_x)
+            spi, cpi = sincospi(abs_nu)
+            return (cpi*Jnu - spi*Ynu) * sg
+        end
+    end
+end
+
+function besselj_positive_args(nu::Real, x::T) where T
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
 
@@ -172,10 +202,9 @@ function _besselj(nu, x)
     if nu >= x
         nu_shift = ceil(Int, debye_cutoff - nu)
         v = nu + nu_shift
-        arr = range(v, stop = nu, length = nu_shift + 1)
         jnu = besseljy_debye(v, x)[1]
         jnup1 = besseljy_debye(v+1, x)[1]
-        return besselj_down_recurrence(x, jnu, jnup1, arr)[2]
+        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[1]
     end
 
     # at this point x > nu and  x < nu * 1.65
@@ -192,14 +221,13 @@ function _besselj(nu, x)
         v2 = nu - nu_shift
         jnu = besseljy_large_argument(v2, x)[1]
         jnum1 = besseljy_large_argument(v2 - 1, x)[1]
-        return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
+        return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[1]
     else
         nu_shift = ceil(Int, debye_diff)
         v = nu + nu_shift
-        arr = range(v, stop = nu, length = nu_shift + 1)
         jnu = besseljy_debye(v, x)[1]
         jnup1 = besseljy_debye(v+1, x)[1]
-        return besselj_down_recurrence(x, jnu, jnup1, arr)[2]
+        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[1]
     end
 end
 
@@ -238,67 +266,3 @@ function log_besselj_small_arguments_orders(v, x::T) where T
     logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
     return exp(logout)
 end
-
-# For 0.0 <= x < 171.5
-# Mean ULP = 0.55
-# Max ULP = 2.4
-# Adapted from Cephes Mathematical Library (MIT license https://en.smath.com/view/CephesMathLibrary/license) by Stephen L. Moshier
-function gamma(x)
-    if x > 11.5
-        return large_gamma(x)
-    elseif x < 0.0
-        #p = floor(x)
-        #isequal(p, abs(x)) && return throw(DomainError(x, "NaN result for non-NaN input."))
-        # need reflection formula
-        return throw(DomainError(x, "Negative numbers are currently not implemented"))
-    elseif x <= 11.5
-        return small_gamma(x)
-    elseif isnan(x)
-        return x
-    end
-end
-function large_gamma(x)
-    isinf(x) && return x
-    T = Float64
-    w = inv(x)
-    s = (
-        8.333333333333331800504e-2, 3.472222222230075327854e-3, -2.681327161876304418288e-3, -2.294719747873185405699e-4,
-        7.840334842744753003862e-4, 6.989332260623193171870e-5, -5.950237554056330156018e-4, -2.363848809501759061727e-5,
-        7.147391378143610789273e-4
-    )
-    w = w * evalpoly(w, s) + one(T)
-    # lose precision on following block
-    y = exp((x)) 
-    # avoid overflow
-    v = x^(0.5 * x - 0.25)
-    y = v * (v / y)
-
-    return SQ2PI(T) * y * w
-end
-function small_gamma(x)
-    T = Float64
-    P = (
-        1.000000000000000000009e0, 8.378004301573126728826e-1, 3.629515436640239168939e-1, 1.113062816019361559013e-1,
-        2.385363243461108252554e-2, 4.092666828394035500949e-3, 4.542931960608009155600e-4, 4.212760487471622013093e-5
-    )
-    Q = (
-        9.999999999999999999908e-1, 4.150160950588455434583e-1, -2.243510905670329164562e-1, -4.633887671244534213831e-2,
-        2.773706565840072979165e-2, -7.955933682494738320586e-4, -1.237799246653152231188e-3, 2.346584059160635244282e-4,
-        -1.397148517476170440917e-5
-    )
-
-    z = one(T)
-    while x >= 3.0
-        x -= one(T)
-        z *= x
-    end
-    while x < 2.0
-        z /= x
-        x += one(T)
-    end
-
-    x -= T(2)
-    p = evalpoly(x, P)
-    q = evalpoly(x, Q)
-    return z * p / q
-end

src/besselj.jl

@@ -161,7 +161,7 @@ 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 up_recurrence(x, besselk0(x), besselk1(x), nu)[1]
+        return besselk_up_recurrence(x, besselk1(x), besselk0(x), 1, nu)[1]
     else
         return besselk_large_orders(nu, x)
     end
@@ -175,7 +175,7 @@ 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 up_recurrence(x, besselk0x(x), besselk1x(x), nu)[1]
+        return besselk_up_recurrence(x, besselk1x(x), besselk0x(x), 1, nu)[1]
     else
         return besselk_large_orders_scaled(nu, x)
     end