add support for negative args

Author: heltonmc

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))

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)
 

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

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
 
- 
+ ### 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

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