update docs and bessely0,1

Author: heltonmc

src/besselj.jl

@@ -10,13 +10,13 @@
 #    Polynomial coefficients are from [1] which is based on [2]
 #    For tiny arugments the power series expansion is used.
 #
-#    Branch 2: 5.0 < x < 75.0
+#    Branch 2: 5.0 < x < 25.0
 #              besselj0 = sqrt(2/(pi*x))*(cos(x - pi/4)*R7(x) - sin(x - pi/4)*R8(x))
 #    Hankel's asymptotic expansion is used
 #    where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
 #    See section 4 of [3] for more details and [1] for coefficients of polynomials
 # 
-#   Branch 3: x >= 75.0
+#   Branch 3: x >= 25.0
 #              besselj0 = sqrt(2/(pi*x))*beta(x)*(cos(x - pi/4 - alpha(x))
 #   See modified expansions given in [3]. Exact coefficients are used
 #
@@ -28,6 +28,12 @@
 # [3] Harrison, John. "Fast and accurate Bessel function computation." 
 #     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
 #
+
+"""
+    besselj0(x::T) where T <: Union{Float32, Float64}
+
+Bessel function of the first kind of order zero, ``J_0(x)``.
+"""
 function besselj0(x::Float64)
     T = Float64
     x = abs(x)
@@ -54,12 +60,12 @@ function besselj0(x::Float64)
         p = p * sc[2] - w * q * sc[1]
         return p * SQ2OPI(T) / sqrt(x)
     else
-        if x < 75.0
+        if x < 120.0
             p = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
             q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
         else
-            p = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288)
-            q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296)
+            p = (one(T), -1/16, 53/512, -4447/8192)
+            q = (-1/8, 25/384, -1073/5120, 375733/229376)
         end
         xinv = inv(x)
         x2 = xinv*xinv
@@ -74,8 +80,6 @@ function besselj0(x::Float64)
         return a * b
     end
 end
-
-
 function besselj0(x::Float32)
     T = Float32
     x = abs(x)
@@ -100,6 +104,11 @@ function besselj0(x::Float32)
     end
 end
 
+"""
+    besselj1(x::T) where T <: Union{Float32, Float64}
+
+Bessel function of the first kind of order one, ``J_1(x)``.
+"""
 function besselj1(x::Float64)
     T = Float64
     x = abs(x)

src/bessely.jl

@@ -10,13 +10,13 @@
 #    Polynomial coefficients are from [1] which is based on [2].
 #    For tiny arugments the power series expansion is used.
 #
-#    Branch 2: 5.0 < x < 75.0
+#    Branch 2: 5.0 < x < 25.0
 #              bessely0 = sqrt(2/(pi*x))*(sin(x - pi/4)*R7(x) - cos(x - pi/4)*R8(x))
 #    Hankel's asymptotic expansion is used
 #    where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
 #    See section 4 of [3] for more details and [1] for coefficients of polynomials
 # 
-#   Branch 3: x >= 75.0
+#   Branch 3: x >= 25.0
 #              bessely0 = sqrt(2/(pi*x))*beta(x)*(sin(x - pi/4 - alpha(x))
 #   See modified expansions given in [3]. Exact coefficients are used.
 #
@@ -28,6 +28,12 @@
 # [3] Harrison, John. "Fast and accurate Bessel function computation." 
 #     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
 #
+
+"""
+    bessely0(x::T) where T <: Union{Float32, Float64}
+
+Bessel function of the second kind of order zero, ``Y_0(x)``.
+"""
 function bessely0(x::T) where T <: Union{Float32, Float64}
     if x <= zero(x)
         if iszero(x)
@@ -42,29 +48,33 @@ function bessely0(x::T) where T <: Union{Float32, Float64}
 end
 function _bessely0_compute(x::Float64)
     T = Float64
-    if x <= 5
+    if x <= 5.0
         z = x * x
         w = evalpoly(z, YP_y0(T)) / evalpoly(z, YQ_y0(T))
         w += TWOOPI(T) * log(x) * besselj0(x)
         return w
-    elseif x < 75.0
-        w = T(5) / x
-        z = w*w
+    elseif x < 25.0
+        w = 5.0 / x
+        z = w * w
         p = evalpoly(z, PP_y0(T)) / evalpoly(z, PQ_y0(T))
         q = evalpoly(z, QP_y0(T)) / evalpoly(z, QQ_y0(T))
         xn = x - PIO4(T)
         sc = sincos(xn)
         p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
+        if x < 120.0
+            p = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
+            q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296, 24713030909/46137344, -7780757249041/436207616)
+        else
+            p = (one(T), -1/16, 53/512, -4447/8192)
+            q = (-1/8, 25/384, -1073/5120, 375733/229376)
+        end
         xinv = inv(x)
         x2 = xinv*xinv
 
-        p = (one(T), -1/16, 53/512, -4447/8192, 5066403/524288)
         p = evalpoly(x2, p)
         a = SQ2OPI(T) * sqrt(xinv) * p
-
-        q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296)
         xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))
 
         # the following computes b = sin(x + xn) more accurately
@@ -91,6 +101,12 @@ function _bessely0_compute(x::Float32)
         return p
     end
 end
+
+"""
+    bessely1(x::T) where T <: Union{Float32, Float64}
+
+Bessel function of the second kind of order one, ``Y_1(x)``.
+"""
 function bessely1(x::T) where T <: Union{Float32, Float64}
     if x <= zero(x)
         if iszero(x)
@@ -103,16 +119,15 @@ function bessely1(x::T) where T <: Union{Float32, Float64}
     end
     return _bessely1_compute(x)
 end
-
 function _bessely1_compute(x::Float64)
     T = Float64
     if x <= 5
         z = x * x
         w = x * (evalpoly(z, YP_y1(T)) / evalpoly(z, YQ_y1(T)))
         w += TWOOPI(T) * (besselj1(x) * log(x) - inv(x))
         return w
-    elseif x < 75.0
-        w = T(5) / x
+    elseif x < 25.0
+        w = 5.0 / x
         z = w * w
         p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
         q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
@@ -121,14 +136,18 @@ function _bessely1_compute(x::Float64)
         p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
+        if x < 120.0
+            p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296) 
+            q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144, -30413055339/46137344, 9228545313147/436207616)
+        else
+            p = (one(T), 3/16, -99/512, 6597/8192)
+            q = (3/8, -21/128, 1899/5120, -543483/229376)
+        end
         xinv = inv(x)
         x2 = xinv*xinv
 
-        p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288)
         p = evalpoly(x2, p)
         a = SQ2OPI(T) * sqrt(xinv) * p
-
-        q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144)
         xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
 
         # the following computes b = sin(x + xn) more accurately
@@ -137,7 +156,6 @@ function _bessely1_compute(x::Float64)
         return a * b
     end
 end
-
 function _bessely1_compute(x::Float32)
     T = Float32
     if x <= 2.0f0