remove branch, update coeff

Author: heltonmc

src/besseli.jl

@@ -11,18 +11,14 @@ http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bess
 Actual coefficients used are from the boost math library.
 https://github.com/boostorg/math/tree/develop/include/boost/math/special_functions/detail
 =#
-
 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)
-    elseif x < branch
-        return exp(x) * evalpoly(inv(x), P2_i0(T)) / sqrt(x)
+        return muladd(a, evalpoly(a, besseli0_small_coefs(T)), 1)
     else
         a = exp(x / 2)
-        s = a * evalpoly(inv(x), P3_i0(T)) / sqrt(x)
+        s = a * evalpoly(inv(x), besseli0_med_coefs(T)) / sqrt(x)
         return a * s
     end
 end
@@ -31,23 +27,23 @@ function besseli0x(x::T) where T <: Union{Float32, Float64}
     x = abs(x)
     if x < 7.75
         a = x * x / 4
-        return muladd(a, evalpoly(a, P1_i0(T)), 1) * exp(-x)
+        return muladd(a, evalpoly(a, besseli0_small_coefs(T)), 1) * exp(-x)
     elseif x < branch
-        return evalpoly(inv(x), P2_i0(T)) / sqrt(x)
+        return evalpoly(inv(x), besseli0_med_coefs(T)) / sqrt(x)
     else
-        return evalpoly(inv(x), P3_i0(T)) / sqrt(x)
+        return evalpoly(inv(x), besseli0_large_coefs(T)) / sqrt(x)
     end
 end
 function besseli1(x::Float32)
     T = Float32
     z = abs(x)
     if z < 7.75
         a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
+        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), P2_i1(T)) / sqrt(z)
+        s = a * evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
         z =  a * s
     end
     if x < zero(x)
@@ -60,13 +56,11 @@ function besseli1(x::Float64)
     z = abs(x)
     if z < 7.75
         a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
+        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
         z = z * evalpoly(a, inner) / 2
-    elseif z < 500
-        return exp(z) * evalpoly(inv(z), P2_i1(T)) / sqrt(z)
     else
         a = exp(z / 2)
-        s = a * evalpoly(inv(z), P3_i1(T)) / sqrt(z)
+        s = a * evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
         z =  a * s
     end
     if x < zero(x)
@@ -79,10 +73,10 @@ function besseli1x(x::Float32)
     z = abs(x)
     if z < 7.75
         a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
+        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), P2_i1(T)) / sqrt(z)
+        z =  evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
     end
     if x < zero(x)
         z = -z
@@ -94,12 +88,12 @@ function besseli1x(x::Float64)
     z = abs(x)
     if z < 7.75
         a = z * z / 4
-        inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
+        inner = (one(T), T(0.5), evalpoly(a, besseli1_small_coefs(T)))
         z = z * evalpoly(a, inner) / 2 * exp(-z)
     elseif z < 500
-        z = evalpoly(inv(z), P2_i1(T)) / sqrt(z)
+        z = evalpoly(inv(z), besseli1_med_coefs(T)) / sqrt(z)
     else
-        z = evalpoly(inv(z), P3_i1(T)) / sqrt(z)
+        z = evalpoly(inv(z), besseli1_large_coefs(T)) / sqrt(z)
     end
     if x < zero(x)
         z = -z

src/constants.jl

@@ -104,23 +104,23 @@ const QQ_j0(::Type{Float64}) = (
     7.24046774195652478189E3, 3.88240183605401609683E3, 8.56430025976980587198E2,
     6.43178256118178023184E1, 1.00000000000000000000E0
 )
-const P1_i0(::Type{Float32}) = (
+const besseli0_small_coefs(::Type{Float32}) = (
     1.00000003928615375f0, 2.49999576572179639f-1, 2.77785268558399407f-2,
     1.73560257755821695f-3, 6.96166518788906424f-5, 1.89645733877137904f-6,  
     4.29455004657565361f-8, 3.90565476357034480f-10, 1.48095934745267240f-11
 )
-const P1_i0(::Type{Float64}) = (
+const besseli0_small_coefs(::Type{Float64}) = (
    1.00000000000000000e0, 2.49999999999999909e-1, 2.77777777777782257e-2,
    1.73611111111023792e-3, 6.94444444453352521e-5, 1.92901234513219920e-6,
    3.93675991102510739e-8, 6.15118672704439289e-10, 7.59407002058973446e-12,
    7.59389793369836367e-14, 6.27767773636292611e-16, 4.34709704153272287e-18,
    2.63417742690109154e-20, 1.13943037744822825e-22, 9.07926920085624812e-25
 )
-const P2_i0(::Type{Float32}) = (
+const besseli0_med_coefs(::Type{Float32}) = (
    3.98942651588301770f-1, 4.98327234176892844f-2, 2.91866904423115499f-2,
    1.35614940793742178f-2, 1.31409251787866793f-1
 )
-const P2_i0(::Type{Float64}) = (
+const besseli0_med_coefs(::Type{Float64}) = (
    3.98942280401425088e-01, 4.98677850604961985e-02, 2.80506233928312623e-02,
    2.92211225166047873e-02, 4.44207299493659561e-02, 1.30970574605856719e-01,
    -3.35052280231727022e+00, 2.33025711583514727e+02, -1.13366350697172355e+04,
@@ -130,29 +130,29 @@ const P2_i0(::Type{Float64}) = (
    -8.13426467865659318e+14, 2.02391097391687777e+15, -3.08675715295370878e+15,
    2.17587543863819074e+15
 )
-const P3_i0(::Type{Float32}) =  (
+const besseli0_large_coefs(::Type{Float32}) =  (
    3.98942391532752700f-1, 4.98455950638200020f-2, 2.94835666900682535f-2
 )
-const P3_i0(::Type{Float64}) = (
+const besseli0_large_coefs(::Type{Float64}) = (
    3.98942280401432905e-01, 4.98677850491434560e-02, 2.80506308916506102e-02,
    2.92179096853915176e-02, 4.53371208762579442e-02
 )
-const P1_i1(::Type{Float32}) = (
+const besseli1_small_coefs(::Type{Float32}) = (
     8.333333221f-2, 6.944453712f-3, 3.472097211f-4, 1.158047174f-5,
     2.739745142f-7, 5.135884609f-9, 5.262251502f-11, 1.331933703f-12
 )
-const P2_i1(::Type{Float32}) = (
+const besseli1_med_coefs(::Type{Float32}) = (
     3.98942115977513013f-1, -1.49581264836620262f-1, -4.76475741878486795f-2,
     -2.65157315524784407f-2, -1.47148600683672014f-1
 )
-const P1_i1(::Type{Float64}) = (
+const besseli1_small_coefs(::Type{Float64}) = (
     8.333333333333333803e-02, 6.944444444444341983e-03, 3.472222222225921045e-04,
     1.157407407354987232e-05, 2.755731926254790268e-07, 4.920949692800671435e-09,
     6.834657311305621830e-11, 7.593969849687574339e-13, 6.904822652741917551e-15,
     5.220157095351373194e-17, 3.410720494727771276e-19, 1.625212890947171108e-21,
     1.332898928162290861e-23
 )
-const P2_i1(::Type{Float64}) = (
+const besseli1_med_coefs(::Type{Float64}) = (
     3.989422804014406054e-01, -1.496033551613111533e-01, -4.675104253598537322e-02,
     -4.090895951581637791e-02, -5.719036414430205390e-02, -1.528189554374492735e-01,
     3.458284470977172076e+00, -2.426181371595021021e+02, 1.178785865993440669e+04,
@@ -162,7 +162,7 @@ const P2_i1(::Type{Float64}) = (
     8.325554073334618015e+14, -2.067285045778906105e+15, 3.146401654361325073e+15,
     -2.213318202179221945e+15
 )
-const P3_i1(::Type{Float64}) = (
+const besseli1_large_coefs(::Type{Float64}) = (
     3.989422804014314820e-01, -1.496033551467584157e-01, -4.675105322571775911e-02,
     -4.090421597376992892e-02, -5.843630344778927582e-02
 )