add neg orders and args

Author: heltonmc

src/Bessels.jl

@@ -4,6 +4,7 @@ using SpecialFunctions: loggamma
 
 export besselj0
 export besselj1
+export besselj
 
 export bessely0
 export bessely1

src/besselj.jl

@@ -150,13 +150,40 @@ 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
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+
+    Jnu = besselj_positive_args(abs_nu, abs_x)
+
+    if nu >= zero(T)
+        if x >= zero(T)
+            return Jnu
+        else
+            if isinteger(abs_nu)
+                return (-1)^Int(abs_nu) * Jnu
+            else
+                return cispi(abs_nu)*Jnu
+            end
+        end
+    else
+        if x >= zero(T)
+            isinteger(nu) && return (-1)^Int(abs_nu) * Jnu
+            Ynu = bessely_positive_args(abs_nu, abs_x)
+            spi, cpi = sincospi(abs_nu)
+            return cpi*Jnu - spi*Ynu
+        else
+            Ynu = bessely_positive_args(abs_nu, abs_x)
+            spi, cpi = sincospi(abs_nu)
+            out = cpi*Jnu - spi*Ynu
+            isinteger(nu) && return (-1)^Int(abs_nu) * out
+            return cispi(nu)*out
+        end
+    end
+end
+
+function besselj_positive_args(nu::Real, x::T) where T
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
 

src/bessely.jl

@@ -187,7 +187,40 @@ end
 Bessel function of the first kind of order nu, ``J_{nu}(x)``.
 Nu must be real.
 """
-function bessely(nu, x::T) where T
+
+
+function bessely(nu::Real, x::T) where T
+    abs_nu = abs(nu)
+    abs_x = abs(x)
+
+    Ynu = bessely_positive_args(abs_nu, abs_x)
+
+    if nu > zero(T)
+        if x > zero(T)
+            return Ynu
+        else
+            Jnu = besselj_positive_args(abs_nu, abs_x)
+            return cispi(-nu)*Ynu + 2im*cospi(nu)*Jnu
+        end
+    elseif nu < zero(T)
+        if x > zero(T)
+            isinteger(nu) && return (-1)^Int(abs_nu) * Ynu
+            Jnu = besselj_positive_args(abs_nu, abs_x)
+            spi, cpi = sincospi(abs_nu)
+            return cpi*Ynu + spi*Jnu
+        else
+            Jnu = besselj_positive_args(abs_nu, abs_x)
+            isinteger(nu) && return Ynu + 2im * Jnu
+            spi, cpi = sincospi(abs_nu)
+            out = cpi * (cispi(nu)*Ynu + 2im * cpi * Jnu)
+            return out + spi * cispi(abs_nu) * Jnu
+        end
+    else
+        return -T(Inf)
+    end
+end
+
+function bessely_positive_args(nu, x::T) where T
 
     # use forward recurrence if nu is an integer up until it becomes inefficient
     (isinteger(nu) && nu < 250) && return besselj_up_recurrence(x, bessely1(x), bessely0(x), 1, nu)[1]

test/besselj_test.jl

@@ -84,7 +84,7 @@ nu = [2, 4, 6, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100]
 for v in nu, xx in x
     xx *= v
     sf = SpecialFunctions.besselj(BigFloat(v), BigFloat(xx))
-    @test isapprox(Bessels._besselj(v, xx), Float64(sf), rtol=5e-14)
+    @test isapprox(besselj(v, xx), Float64(sf), rtol=5e-14)
 end
 
 # test half orders (SpecialFunctions does not give big float precision)
@@ -94,20 +94,38 @@ x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97,
 nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
 for v in nu, xx in x
     xx *= v
-    @test isapprox(Bessels._besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=1e-12)
+    @test isapprox(besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=1e-12)
 end
 
 ## test large orders
 nu = [150, 165.2, 200.0, 300.0, 500.0, 1000.0, 5000.2, 10000.0, 50000.0]
 x = [0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92,0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,0.995, 0.999, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2]
 for v in nu, xx in x
     xx *= v
-    @test isapprox(Bessels._besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=5e-11)
+    @test isapprox(besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=5e-11)
 end
 
 ## test large arguments
-@test isapprox(Bessels._besselj(10.0, 150.0), SpecialFunctions.besselj(10.0, 150.0), rtol=1e-12)
+@test isapprox(besselj(10.0, 150.0), SpecialFunctions.besselj(10.0, 150.0), rtol=1e-12)
 
 # test BigFloat for single point
-@test isapprox(Bessels._besselj(big"2000", big"1500.0"), SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
-@test isapprox(Bessels._besselj(big"20", big"1500.0"), SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
+@test isapprox(besselj(big"2000", big"1500.0"), SpecialFunctions.besselj(big"2000", big"1500"), rtol=5e-20)
+@test isapprox(besselj(big"20", big"1500.0"), SpecialFunctions.besselj(big"20", big"1500"), rtol=5e-20)
+
+
+# need to test accuracy of negative orders and negative arguments and all combinations within
+# SpecialFunctions.jl doesn't provide these so will hand check against hard values
+# values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica
+# need to also account for different branches when nu isa integer
+nu = -9.102; x = -12.48
+@test isapprox(besselj(nu, x), 0.09842356047575545808128 -0.03266486217437818486161im, rtol=8e-15)
+nu = -5.0; x = -5.1
+@test isapprox(besselj(nu, x), 0.2740038554704588327387, rtol=8e-15)
+nu = -7.3; x = 19.1
+@test isapprox(besselj(nu, x), 0.1848055978553359009813, rtol=8e-15)
+nu = -14.0; x = 21.3
+@test isapprox(besselj(nu, x), -0.1962844898264965120021, rtol=8e-15)
+nu = 13.0; x = -8.5
+@test isapprox(besselj(nu, x), -0.006128034621313167000171, rtol=8e-15)
+nu = 17.45; x = -16.23
+@test isapprox(besselj(nu, x), -0.01607335977752705869797 -0.1014831996412783806255im, rtol=8e-15)

test/bessely_test.jl

@@ -86,3 +86,20 @@ for v in nu, xx in x
     xx *= v
     @test isapprox(bessely(v, xx), SpecialFunctions.bessely(v, xx), rtol=5e-12)
 end
+
+# need to test accuracy of negative orders and negative arguments and all combinations within
+# SpecialFunctions.jl doesn't provide these so will hand check against hard values
+# values taken from https://keisan.casio.com/exec/system/1180573474 which match mathematica
+# need to also account for different branches when nu isa integer
+nu = -2.3; x = -2.4
+@test isapprox(bessely(nu, x), -0.0179769671833457636186 + 0.7394120337538928700168im, rtol=8e-15)
+nu = -4.0; x = -12.6
+@test isapprox(bessely(nu, x), -0.02845106816742465563357 +0.4577922229605476882792im, rtol=8e-15)
+nu = -6.2; x = 18.6
+@test isapprox(bessely(nu, x), -0.05880321550673669650027, rtol=8e-15)
+nu = -8.0; x = 23.2
+@test isapprox(bessely(nu, x),-0.166071277370329242677, rtol=8e-15)
+nu = 11.0; x = -8.2
+@test isapprox(bessely(nu, x),1.438494049708244558901 -0.06222860017222637350092im, rtol=8e-15)
+nu = 13.678; x = -12.98
+@test isapprox(bessely(nu, x),-0.2227392320508850571009 -0.2085585256158188848322im, rtol=8e-15)