add domain errors bessely

Author: heltonmc

src/j0.jl

@@ -20,6 +20,13 @@
 #
 #
 function besselj0(x::Float64)
+    x = abs(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
+    end
+
     if x <= 5.0
         z = x * x
         if x < 1.0e-5
@@ -78,8 +85,11 @@ end
 
 
 function besselj0(x::Float32)
-    if x < 0.0f0
-        x *= -1
+    x = abs(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
     end
 
     if x <= 2.0f0
@@ -97,7 +107,7 @@ function besselj0(x::Float32)
         p = (z - DR1) * evalpoly(z, JP)
         return p
     else
-        q = 1.0f0 / x
+        q = inv(x)
         w = sqrt(q)
 
         MO = (
@@ -142,17 +152,22 @@ end
 # Q0(x) = 1/x (-.125 + 1/x^2 R(1/x^2))
 #
 function besselj0(x::BigFloat)
-    xx = abs(x)
-    if iszero(xx)
-        return one(BigFloat)
-    elseif xx <= 2.0
-        z = xx * xx
+    x = abs(x)
+    T = eltype(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 2.0
+        z = x * x
         p = z * z * evalpoly(z, J0_2N) / evalpoly(z, J0_2D)
         p -= 0.25 * z
-        p += 1.0
+        p += one(T)
         return p
     else
-        xinv = 1.0 / xx
+        xinv = inv(x)
         z  = xinv * xinv
         if xinv <= 0.25
             if xinv <= 0.125
@@ -191,17 +206,17 @@ function besselj0(x::BigFloat)
         p = 1.0 + z * p
         q = z * xinv * q
         q = q - 0.125 * xinv
-        c = cos(xx)
-        s = sin(xx)
+        c = cos(x)
+        s = sin(x)
         ss = s - c
         cc = s + c
-        z = -cos(xx + xx)
+        z = -cos(x + x)
         if (s * c) < 0
             cc = z / ss;
         else
             ss = z / cc;
         end
-        z = ONEOSQPI(BigFloat) * (p * cc - q * ss) / sqrt(xx)
+        z = ONEOSQPI(T) * (p * cc - q * ss) / sqrt(x)
         return z
     end
 end

src/j1.jl

@@ -1,7 +1,11 @@
 function besselj1(x::Float32)
-    if x < zero(x)
-        x = -x
+    x = abs(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
     end
+
     if x <= 2.0f0
         z = x * x
         Z1 = 1.46819706421238932572f1
@@ -19,15 +23,17 @@ function besselj1(x::Float32)
 end
 
 function besselj1(x::Float64)
-    w = x
-    if x < 0
-        w = -x
+    x = abs(x)
+    if iszero(x)
+        return one(x)
+    elseif isinf(x)
+        return zero(x)
     end
 
-    if w <= 5.0
+    if x <= 5.0
         z = x * x
         w = evalpoly(z, RP) / evalpoly(z, RQ)
-        w = w * x * (z - 1.46819706421238932572E1) * (z - 4.92184563216946036703E1)
+        w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
         return w
     else
         w = 5.0 / x
@@ -41,10 +47,17 @@ function besselj1(x::Float64)
 end
 
 function bessely1(x::Float32)
-    if x <= 2.0f0
-        if x <= zero(x)
-            return error("Domain error: x <= 0")
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf32
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
         end
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 2.0f0
         z = x * x
         YO1 =  4.66539330185668857532f0
         w = (z - YO1) * x * evalpoly(z, YP32)
@@ -62,10 +75,17 @@ function bessely1(x::Float32)
 end
 
 function bessely1(x::Float64)
-    if x <= 5.0
-        if x< 0.0
-            return NaN
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf64
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
         end
+    elseif isinf(x)
+        return zero(x)
+    end
+
+    if x <= 5.0
         z = x * x
         w = x * (evalpoly(z, YP) / evalpoly(z, YQ))
         w += TWOOPI(Float64) * (besselj1(x) * log(x) - inv(x))

src/y0.jl

@@ -21,10 +21,16 @@
 # Ported to Julia in December 2021 by Michael Helton
 #
 function bessely0(x::Float64)
-    if x <= 5.0
-        if x <= 0.0
-            return NaN64
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf64
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
         end
+    elseif isinf(x)
+        return zero(x)
+    end
+    if x <= 5.0
         z = x * x
 
         YP = (
@@ -78,14 +84,16 @@ end
 
 
 function bessely0(x::Float32)
-    if x <= 2.0f0
-        if x <= 0.0f0
-            if x == 0.0f0
-                return -Inf32
-            end
-            return NaN32
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf32
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
         end
-
+    elseif isinf(x)
+        return zero(x)
+    end
+    if x <= 2.0f0
         z = x * x
         YZ1 =  0.43221455686510834878f0
         YP = (
@@ -122,16 +130,16 @@ end
 
 
 function bessely0(x::BigFloat)
-    if isnan(x)
-        return NaN
-    elseif x < zero(x)
-        return NaN
-    elseif iszero(x)
-        return -Inf
+    if x <= zero(x)
+        if iszero(x)
+            return -Inf
+        else
+            return throw(DomainError(x, "NaN result for non-NaN input."))
+        end
     elseif isinf(x)
         return zero(x)
     end
-
+  
     xx = abs(x)
 
     if xx <= 2.0
@@ -141,7 +149,7 @@ function bessely0(x::BigFloat)
         return p
     end
 
-    xinv = 1.0 / xx
+    xinv = inv(xx)
     z = xinv * xinv
     if xinv <= 0.25
         if xinv <= 0.125

test/bessely_test.jl

@@ -16,6 +16,14 @@ y0_big = bessely0.(big.(x))
 
 @test isapprox(y0_big, y0_SpecialFunctions, atol=1.5e-34)
 
+@test_throws DomainError bessely0(-1.0)
+
+@test bessely0(zero(Float32)) == -Inf32
+@test bessely0(zero(Float64)) == -Inf64
+@test bessely0(zero(BigFloat)) == -Inf
+
+@test bessely0(Inf32) == zero(Float32)
+@test bessely0(Inf64) == zero(Float64)
 
 y1_SpecialFunctions = SpecialFunctions.bessely1.(big.(x))
 @assert y1_SpecialFunctions[1] isa BigFloat
@@ -31,4 +39,11 @@ y1_32 = bessely1.(Float32.(x))
 @test y1_64 ≈ y1_SpecialFunctions
 @test y1_32 ≈ y1_SpecialFunctions
 
+@test_throws DomainError bessely1(-1.0)
+
+@test bessely1(zero(Float32)) == -Inf32
+@test bessely1(zero(Float64)) == -Inf64
+
+@test bessely1(Inf32) == zero(Float32)
+@test bessely1(Inf64) == zero(Float64)