use sincos for medium args

Author: heltonmc

src/besselj.jl

@@ -114,7 +114,8 @@ function besselj1(x::Float64)
         p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
         q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
         xn = x - THPIO4(T)
-        p = p * cos(xn) - w * q * sin(xn)
+        sc = sincos(xn)
+        p = p * sc[2] - w * q * sc[1]
         return p * SQ2OPI(T) / sqrt(x)
     else
         xinv = inv(x)

src/bessely.jl

@@ -29,7 +29,8 @@ function _bessely0_compute(x::Float64)
         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)
-        p = p * sin(xn) + w * q * cos(xn);
+        sc = sincos(xn)
+        p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
         xinv = inv(x)
@@ -89,13 +90,14 @@ function _bessely1_compute(x::Float64)
         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 < 100.0
+    elseif x < 75.0
         w = T(5) / 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))
         xn = x - THPIO4(T)
-        p = p * sin(xn) + w * q * cos(xn)
+        sc = sincos(xn)
+        p = p * sc[1] + w * q * sc[2]
         return p * SQ2OPI(T) / sqrt(x)
     else
         xinv = inv(x)