remove arctan function

Author: heltonmc

src/besselj.jl

@@ -53,27 +53,19 @@ function besselj0(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), -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)
-        return SQ2OPI(T) * sqrt(xinv) * beta * cos_sum(x, - T(pi)/4 - alpha)
     else
+        if x < 75.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)
+        end
         xinv = inv(x)
         x2 = xinv*xinv
 
-        p = (one(T), -1/16, 53/512, -4447/8192, 3066403/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 = cos(x + xn) more accurately
@@ -82,6 +74,8 @@ function besselj0(x::Float64)
         return a * b
     end
 end
+
+
 function besselj0(x::Float32)
     T = Float32
     x = abs(x)
@@ -125,27 +119,19 @@ 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
+        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 = cos(x + xn) more accurately
@@ -154,7 +140,6 @@ function besselj1(x::Float64)
         return a * b
     end
 end
-
 function besselj1(x::Float32)
     x = abs(x)
     isinf(x) && return zero(x)