add mysin, mycos

Author: heltonmc

src/Bessels.jl

@@ -36,6 +36,7 @@ include("Float128/constants.jl")
 include("math_constants.jl")
 include("U_polynomials.jl")
 include("recurrence.jl")
+include("misc.jl")
 #include("hankel.jl")
 
 end

src/besselj.jl

@@ -11,13 +11,13 @@
 #    See [1] for more details and [2] for coefficients of polynomials.
 #    For tiny arugments the power series expansion is used.
 #
-#    Branch 2: 5.0 < x < 80.0
+#    Branch 2: 5.0 < x < 75.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 >= 80.0
+#   Branch 3: x >= 75.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,7 +28,7 @@
 # [2] Cephes Math Library Release 2.8:  June, 2000 by Stephen L. Moshier
 # [3] Harrison, John. "Fast and accurate Bessel function computation." 
 #     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
-
+#
 function besselj0(x::T) where T
     x = abs(x)
     isinf(x) && return zero(x)
@@ -64,12 +64,9 @@ function besselj0(x::T) where T
         q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296)
         xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))
 
-        # the following computes b = cos(x + xn)
-        # but uses cos(x)*cos(xn) - sin(x)*sin(xn)
-        # to improve accuracy when x >> xn
-        c1 = sincos(x)
-        c2 = sincos(xn)
-        b = c1[2] * c2[2] - c1[1] * c2[1]
+        # the following computes b = cos(x + xn) more accurately
+        # see src/misc.jl
+        b = mycos(x, xn)
         return a * b
     end
 end
@@ -127,12 +124,9 @@ function besselj1(x::Float64)
         q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144)
         xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
 
-        # the following computes b = cos(x + xn)
-        # but uses cos(x)*cos(xn) - sin(x)*sin(xn)
-        # to improve accuracy when x >> xn
-        c1 = sincos(x)
-        c2 = sincos(xn)
-        b = c1[2] * c2[2] - c1[1] * c2[1]
+        # the following computes b = cos(x + xn) more accurately
+        # see src/misc.jl
+        b = mycos(x, xn)
         return a * b
     end
 end

src/bessely.jl

@@ -1,9 +1,34 @@
-#=
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 2000 by Stephen L. Moshier
-https://github.com/jeremybarnes/cephes/blob/master/bessel/j0.c
-https://github.com/jeremybarnes/cephes/blob/master/bessel/j1.c
-=#
+#    Bessel functions of the second kind of order zero and one
+#                       bessely0, bessely1
+#
+#    Calculation of bessely0 is done in three branches using polynomial approximations
+#
+#    Branch 1: x <= 5.0
+#              bessely0 = R(x^2) + 2*log(x)*besselj0(x) / pi
+#    where r1 and r2 are zeros of J0
+#    and P3 and Q8 are a 3 and 8 degree polynomial respectively
+#    Polynomial coefficients are from [1] which is based on [2]
+#    See [1] for more details and [2] for coefficients of polynomials.
+#    For tiny arugments the power series expansion is used.
+#
+#    Branch 2: 5.0 < x < 75.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
+#              bessely0 = sqrt(2/(pi*x))*beta(x)*(sin(x - pi/4 - alpha(x))
+#   See modified expansions given in [3]. Exact coefficients are used.
+#
+#   Calculation of bessely1 is done in a similar way as bessely0.
+#   See [3] for details on similarities.
+# 
+# [1] https://github.com/deepmind/torch-cephes
+# [2] Cephes Math Library Release 2.8:  June, 2000 by Stephen L. Moshier
+# [3] Harrison, John. "Fast and accurate Bessel function computation." 
+#     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
+#
 function bessely0(x::T) where T <: Union{Float32, Float64}
     if x <= zero(x)
         if iszero(x)
@@ -43,12 +68,9 @@ function _bessely0_compute(x::Float64)
         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)
-        # but uses cos(x)*sin(xn) + sin(x)*cos(xn)
-        # to improve accuracy when x >> xn
-        c1 = sincos(x)
-        c2 = sincos(xn)
-        b = c1[2] * c2[1] + c1[1] * c2[2]
+        # the following computes b = sin(x + xn) more accurately
+        # see src/misc.jl
+        b = mysin(x, xn)
         return a * b
     end
 end
@@ -110,12 +132,9 @@ function _bessely1_compute(x::Float64)
         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)
-        # but uses cos(x)*sin(xn) + sin(x)*cos(xn)
-        # to improve accuracy when x >> xn
-        c1 = sincos(x)
-        c2 = sincos(xn)
-        b = c1[2] * c2[1] + c1[1] * c2[2]
+        # the following computes b = sin(x + xn) more accurately
+        # see src/misc.jl
+        b = mysin(x, xn)
         return a * b
     end
 end
@@ -124,7 +143,7 @@ function _bessely1_compute(x::Float32)
     T = Float32
     if x <= 2.0f0
         z = x * x
-        YO1 =  4.66539330185668857532f0
+        YO1 = 4.66539330185668857532f0
         w = (z - YO1) * x * evalpoly(z, YP32)
         w += TWOOPI(Float32) * (besselj1(x) * log(x) - inv(x))
         return w

src/misc.jl

@@ -0,0 +1,39 @@
+# function to more accurately compute cos(x + xn)
+# when x >> xn
+# see https://github.com/heltonmc/Bessels.jl/pull/13
+# written by @oscardssmith
+function mycos(x, xn)
+    s = x + xn
+    n, r = Base.Math.rem_pio2_kernel(s)
+    lo = r.lo - ((s - x) - xn)
+    hi = r.hi + lo
+    y = Base.Math.DoubleFloat64(hi, r.hi-hi+lo)
+    n = n&3
+    if n == 0
+        return Base.Math.cos_kernel(y)
+    elseif n == 1
+        return -Base.Math.sin_kernel(y)
+    elseif n == 2
+        return -Base.Math.cos_kernel(y)
+    else
+        return Base.Math.sin_kernel(y)
+    end
+end
+# function to more accurately compute sin(x + xn)
+function mysin(x, xn)
+    s = x + xn
+    n, r = Base.Math.rem_pio2_kernel(s)
+    lo = r.lo - ((s - x) - xn)
+    hi = r.hi + lo
+    y = Base.Math.DoubleFloat64(hi, r.hi-hi+lo)
+    n = n&3
+    if n == 0
+        return Base.Math.sin_kernel(y)
+    elseif n == 1
+        return Base.Math.cos_kernel(y)
+    elseif n == 2
+        return -Base.Math.sin_kernel(y)
+    else
+        return -Base.Math.cos_kernel(y)
+    end
+end