add neg numbers support

Author: heltonmc

src/besseli.jl

@@ -14,6 +14,7 @@ https://github.com/boostorg/math/tree/develop/include/boost/math/special_functio
 
 function besseli0(x::T) where T <: Union{Float32, Float64}
     T == Float32 ? branch = 50 : branch = 500
+    x = abs(x)
     if x < 7.75
         a = x * x / 4
         return muladd(a, evalpoly(a, P1_i0(T)), 1)
@@ -27,6 +28,7 @@ function besseli0(x::T) where T <: Union{Float32, Float64}
 end
 function besseli0x(x::T) where T <: Union{Float32, Float64}
     T == Float32 ? branch = 50 : branch = 500
+    x = abs(x)
     if x < 7.75
         a = x * x / 4
         return muladd(a, evalpoly(a, P1_i0(T)), 1) * exp(-x)
@@ -38,49 +40,69 @@ function besseli0x(x::T) where T <: Union{Float32, Float64}
 end
 function besseli1(x::Float32)
     T = Float32
-    if x < 7.75
-        a = x * x / 4
+    z = abs(x)
+    if z < 7.75
+        a = z * z / 4
         inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
-        return x * evalpoly(a, inner) / 2
+        z = z * evalpoly(a, inner) / 2
     else
-        a = exp(x / 2)
-        s = a * evalpoly(inv(x), P2_i1(T)) / sqrt(x)
-        return a * s
+        a = exp(z / 2)
+        s = a * evalpoly(inv(z), P2_i1(T)) / sqrt(z)
+        z =  a * s
     end
+    if x < zero(x)
+        z = -z
+    end
+    return z
 end
 function besseli1(x::Float64)
     T = Float64
-    if x < 7.75
-        a = x * x / 4
+    z = abs(x)
+    if z < 7.75
+        a = z * z / 4
         inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
-        return x * evalpoly(a, inner) / 2
-    elseif x < 500
-        return exp(x) * evalpoly(inv(x), P2_i1(T)) / sqrt(x)
+        z = z * evalpoly(a, inner) / 2
+    elseif z < 500
+        return exp(z) * evalpoly(inv(z), P2_i1(T)) / sqrt(z)
     else
-        a = exp(x / 2)
-        s = a * evalpoly(inv(x), P3_i1(T)) / sqrt(x)
-        return a * s
+        a = exp(z / 2)
+        s = a * evalpoly(inv(z), P3_i1(T)) / sqrt(z)
+        z =  a * s
+    end
+    if x < zero(x)
+        z = -z
     end
+    return z
 end
 function besseli1x(x::Float32)
     T = Float32
-    if x < 7.75
-        a = x * x / 4
+    z = abs(x)
+    if z < 7.75
+        a = z * z / 4
         inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
-        return x * evalpoly(a, inner) / 2 * exp(-x)
+        z = z * evalpoly(a, inner) / 2 * exp(-z)
     else
-        return evalpoly(inv(x), P2_i1(T)) / sqrt(x)
+        z =  evalpoly(inv(z), P2_i1(T)) / sqrt(z)
     end
+    if x < zero(x)
+        z = -z
+    end
+    return z
 end
 function besseli1x(x::Float64)
     T = Float64
-    if x < 7.75
-        a = x * x / 4
+    z = abs(x)
+    if z < 7.75
+        a = z * z / 4
         inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
-        return x * evalpoly(a, inner) / 2 * exp(-x)
-    elseif x < 500
-        return evalpoly(inv(x), P2_i1(T)) / sqrt(x)
+        z = z * evalpoly(a, inner) / 2 * exp(-z)
+    elseif z < 500
+        z = evalpoly(inv(z), P2_i1(T)) / sqrt(z)
     else
-        return evalpoly(inv(x), P3_i1(T)) / sqrt(x)
+        z = evalpoly(inv(z), P3_i1(T)) / sqrt(z)
+    end
+    if x < zero(x)
+        z = -z
     end
+    return z
 end

test/besseli_test.jl

@@ -1,6 +1,7 @@
 # general array for testing input to SpecialFunctions.jl
-x32 = 1e-6:0.01:60.0
-x64 = 1e-6:0.01:600.0
+x32 = [-1.0; 0.0; 1e-6; 0.01:0.01:60.0]
+x64 = [-1.0; 0.0; 1e-6; 0.01:0.01:600.0]
+
 ### Tests for besseli0
 i032_SpecialFunctions = SpecialFunctions.besseli.(0, x32)
 i064_SpecialFunctions = SpecialFunctions.besseli.(0, x64)