update besselj1

Author: heltonmc

src/besselj.jl

@@ -58,15 +58,12 @@ function besselj0(x::Float64)
         x2 = xinv * xinv
         p = (one(T), -9/128, 3675/32768, - 2401245/4194304, 13043905875/2147483648, - 30241281245175/274877906944, 213786613951685775/70368744177664, -1070401384414690453125/9007199254740992, 57673297952355815927071875/9223372036854775808)
         q = (-1/8, 75/1024, - 59535/262144, 57972915/33554432, - 418854310875/17179869184, 1212400457192925/2199023255552, - 10278202593831046875/562949953421312, 60013837619516978071875/72057594037927936, - 3694483615889146090857721875/73786976294838206464)
-
         p = evalpoly(x2, p)
         q = evalpoly(x2, q) * xinv
+        alpha = atan(-q/p)
 
         pq = (one(T), -1/16, 53/512, -4447/8192, 3066403/524288, -896631415/8388608, 796754802993/268435456, -500528959023471/4294967296)
         beta = evalpoly(x2, pq)
-
-        alpha = atan(-q/p)
-        
         return SQ2OPI(T) * sqrt(xinv) * beta * cos_sum(x, - T(pi)/4 - alpha)
     else
         xinv = inv(x)
@@ -119,7 +116,7 @@ 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
-    elseif x < 75.0
+    elseif x < 25.0
         w = 5.0 / x
         z = w * w
         p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
@@ -128,6 +125,18 @@ function besselj1(x::Float64)
         sc = sincos(xn)
         p = p * sc[2] - w * q * sc[1]
         return p * SQ2OPI(T) / sqrt(x)
+    elseif x < 75.0
+        xinv = inv(x)
+        x2 = xinv * xinv
+        p = (one(T), 15/128, -4725/32768, 2837835/4194304, -14783093325/2147483648, 33424574007825/274877906944, -232376754295310625/70368744177664, 1149690375852815671875/9007199254740992)
+        q = (3/8, -105/1024, 72765/262144, -66891825/33554432, 468131288625/17179869184, -1327867167401775/2199023255552, 11100458801337530625/562949953421312, -64152722972587114490625/72057594037927936)
+        p = evalpoly(x2, p)
+        q = evalpoly(x2, q) * xinv
+        alpha = atan(-q/p)
+
+        pq = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288, 1113686901/8388608, -951148335159/268435456, 581513783771781/4294967296, -3198424285940846836612419/9223372036854775808)
+        beta = evalpoly(x2, pq)
+        return SQ2OPI(T) * sqrt(xinv) * beta * cos_sum(x, - 3*T(pi)/4 - alpha)
     else
         xinv = inv(x)
         x2 = xinv*xinv