add tests

Author: heltonmc

src/besselj.jl

@@ -159,7 +159,7 @@ Nu must be real.
 function _besselj(nu, x)
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
-    if x < 30.0
+    if x < 20.0
         if nu > 60.0
             return log_besselj_small_arguments_orders(nu, x)
         else
@@ -215,9 +215,10 @@ end
 
 # this needs a better way to sum these as it produces large errors
 # only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))
-# power series has premature underflow for large orders 
+# power series has premature underflow for large orders
+# gives wrong answers for x > 20.0 (might want to fix this)
 function besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 1_000
+    MaxIter = 1000
     out = zero(T)
     a = (x/2)^v / gamma(v + one(T))
     t2 = (x/2)^2
@@ -230,10 +231,10 @@ function besselj_small_arguments_orders(v, x::T) where T
 end
 
 # this needs a better way to sum these as it produces large errors
-# perhaps use when v is small i believe v also has to be positive for this to work
+# use when v is large and x is small
 # need for bessely 
 function log_besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 500
+    MaxIter = 1000
     out = zero(T)
     a   = one(T)
     x2 = (x/2)^2
@@ -242,7 +243,6 @@ function log_besselj_small_arguments_orders(v, x::T) where T
         a = -a * x2 * inv((i + one(T)) * (v + i + one(T)))
         (abs(a) < eps(T) * abs(out)) && break
     end
-    @show out
     logout = -loggamma(v + 1) + fma(v, log(x/2), log(out))
     return exp(logout)
 end
@@ -284,6 +284,7 @@ function gamma(x)
     end
 end
 function large_gamma(x)
+    isinf(x) && return x
     T = Float64
     w = inv(x)
     s = (
@@ -294,15 +295,10 @@ function large_gamma(x)
     w = w * evalpoly(w, s) + one(T)
     # lose precision on following block
     y = exp((x)) 
-    if x > 143.01608
-        isinf(x) && return x
-        # avoid overflow
-        v = x^(0.5 * x - 0.25)
-        y = v * (v / y)
-    else
-        y = x^(x - 0.5) / y
-        # (x - 0.5) * log(x) - x
-    end
+    # avoid overflow
+    v = x^(0.5 * x - 0.25)
+    y = v * (v / y)
+    
     return SQ2PI(T) * y * w
 end
 function small_gamma(x)

test/besselj_test.jl

@@ -78,11 +78,25 @@ j1_32 = besselj1.(Float32.(x))
 #@test besselj(-5, 6.1) ≈ SpecialFunctions.besselj(-5, 6.1)
 
 ## test the small value approximation using power series 
-
 nu = [0.5, 1.5, 3.0, 10.0, 22.2, 35.0, 52.1, 100.0, 250.2, 500.0, 1000.0]
-x = [0.001, 1.0, 5.0, 15.0, 29.9]
+x = [0.001, 1.0, 5.0, 15.0, 19.9]
 
 for v in nu, x in x
+    @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
+end
+
+## test the large argument asymptotic expansion
+vnu = ((0.5, 21.0), (0.5, 45.0), (3.5, 21.0), (10.0, 21.0), (21.2, 45.0), (43.2, 90.0), (100.1, 205.2))
+
+for (v, x) in vnu
+    @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
+end
+
+## test the debye uniform asymptotic expansion for x < nu
+vnu = ((60.5, 21.0), (100.5, 45.0), (150.5, 61.0))
+
+for (v, x) in vnu
     @show v, x
     @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
 end
+