add large args j1, y1

Author: heltonmc

src/besselj.jl

@@ -81,14 +81,26 @@ function besselj1(x::Float64)
         w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
         w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
         return w
-    else
+    elseif x < 100.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))
         xn = x - THPIO4(T)
         p = p * cos(xn) - w * q * sin(xn)
         return p * SQ2OPI(T) / sqrt(x)
+    else
+        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))
+        b = cos(x + xn)
+        return a * b
     end
 end
 

src/bessely.jl

@@ -23,14 +23,26 @@ function _bessely0_compute(x::Float64)
         w = evalpoly(z, YP_y0(T)) / evalpoly(z, YQ_y0(T))
         w += TWOOPI(T) * log(x) * besselj0(x)
         return w
-    else
+    elseif x < 100.0
         w = T(5) / 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)
         p = p * sin(xn) + w * q * cos(xn);
         return p * SQ2OPI(T) / sqrt(x)
+    else
+        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))
+        b = sin(x + xn)
+        return a * b
     end
 end
 function _bessely0_compute(x::Float32)
@@ -71,14 +83,26 @@ 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
-    else
+    elseif x < 100.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)
         return p * SQ2OPI(T) / sqrt(x)
+    else
+        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))
+        b = sin(x + xn)
+        return a * b
     end
 end