merge F32 & F64 besseli1

Author: heltonmc

src/besseli.jl

@@ -34,8 +34,7 @@ function besseli0x(x::T) where T <: Union{Float32, Float64}
         return evalpoly(inv(x), besseli0_large_coefs(T)) / sqrt(x)
     end
 end
-function besseli1(x::Float32)
-    T = Float32
+function besseli1(x::T) where T <: Union{Float32, Float64}
     z = abs(x)
     if z < 7.75
         a = z * z / 4
@@ -51,46 +50,14 @@ function besseli1(x::Float32)
     end
     return z
 end
-function besseli1(x::Float64)
-    T = Float64
-    z = abs(x)
-    if z < 7.75
-        a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
-        z = z * evalpoly(a, inner) / 2
-    else
-        a = exp(z / 2)
-        s = a * evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
-        z =  a * s
-    end
-    if x < zero(x)
-        z = -z
-    end
-    return z
-end
-function besseli1x(x::Float32)
-    T = Float32
-    z = abs(x)
-    if z < 7.75
-        a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
-        z = z * evalpoly(a, inner) / 2 * exp(-z)
-    else
-        z =  evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
-    end
-    if x < zero(x)
-        z = -z
-    end
-    return z
-end
-function besseli1x(x::Float64)
-    T = Float64
+function besseli1x(x::T) where T <: Union{Float32, Float64}
+    T == Float32 ? branch = 50 : branch = 500
     z = abs(x)
     if z < 7.75
         a = z * z / 4
         inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
         z = z * evalpoly(a, inner) / 2 * exp(-z)
-    elseif z < 500
+    elseif z < branch
         z = evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
     else
         z = evalpoly(inv(z), besseli1_large_coefs(T)) / sqrt(z)

src/constants.jl

@@ -166,6 +166,10 @@ const besseli1_large_coefs(::Type{Float64}) = (
     3.989422804014314820e-01, -1.496033551467584157e-01, -4.675105322571775911e-02,
     -4.090421597376992892e-02, -5.843630344778927582e-02
 )
+const besseli1_large_coefs(::Type{Float32}) = (
+    3.989422804014314820f-01, -1.496033551467584157f-01, -4.675105322571775911f-02,
+    -4.090421597376992892f-02, -5.843630344778927582f-02
+)
 const A_k0(::Type{Float64}) = (
     1.37446543561352307156E-16, 4.25981614279661018399E-14, 1.03496952576338420167E-11,
     1.90451637722020886025E-9, 2.53479107902614945675E-7, 2.28621210311945178607E-5,