refactor F64, BigFloat

Author: heltonmc

src/asymptotics.jl

@@ -180,13 +180,6 @@ function _α_αp_poly_30(v, x::T) where T
 end
 
 # the following are the coefficients in exact fractions which can be used for arbitrary precision calculations
-function besselj_large_argument(v, x::BigFloat)
-    T = BigFloat
-    α, αp = _α_αp_asymptotic(v, x)
-    xn = fma(v, PIO2(T), PIO4(T))
-    b = SQ2OPI(T) / sqrt(αp * x)
-    return cos(α - xn)*b
-end
 function _α_αp_asymptotic(v, x::BigFloat)
     T = BigFloat
     xinv = inv(x)^2

src/besselj.jl

@@ -203,17 +203,16 @@ end
 
 # for moderate size arguments of x and v this has relative errors ~9e-15
 # for large arguments relative errors ~1e-13
-function besselj_large_argument(v, x::T) where T <: Float64
+function besselj_large_argument(v, x::T) where T
     α, αp = _α_αp_asymptotic(v, x)
-    xn = fma(v, PIO2(T), PIO4(T))
     b = SQ2OPI(T) / sqrt(αp * x)
-    return cos_sum(α, -xn)*b
-end
-function besselj_large_argument(v, x::T) where T <: Float32
-    α, αp = _α_αp_asymptotic(v, x)
-    xn = fma(v, PIO2(T), PIO4(T))
-    b = SQ2OPI(T) / sqrt(αp * x)
-    return cos(Float64(α) - Float64(xn))*b
+
+    S, C = sincos(PIO2(T)*v)
+    Sα, Cα = sincos(α)
+    s1 = (C - S) * Cα
+    s2 = (C + S) * Sα
+
+    return SQ2O2(T) * (s1 + s2) * b
 end
 
 # generally can only use for x < 4.0
@@ -258,14 +257,12 @@ function besselj_debye(v, x)
     S = promote_type(T, Float64)
     x = S(x)
 
-    vmx = fma(v,v, -x^2)
-    vdx = v/x
-    b = sqrt(vmx)
-
-    n = v * log(vdx + sqrt(vdx^2 - 1)) - b
+    vmx = (v + x) * (v − x)
+    vs = sqrt(vmx)
+    n  = fma(v, -log(x / (v + vs)), -vs)
 
-    coef = SQ1O2PI(S) * exp(-n) / sqrt(b)
-    p = v / b
+    coef = SQ1O2PI(S) * exp(-n) / sqrt(vs)
+    p = v / vs
     p2  = v^2 / vmx
 
     return coef * Uk_poly_Jn(p, v, p2, T)

src/math_constants.jl

@@ -5,6 +5,7 @@ const SQPIO2(::Type{BigFloat}) = big"1.25331413731550025120788264240552262650349
 const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008"
 const SQ2OPI(::Type{BigFloat}) = big"0.7978845608028653558798921198687637369517172623298693153318516593413158517986017"
 const PIO4(::Type{BigFloat}) = big"0.7853981633974483096156608458198757210492923498437764552437361480769541015715495"
+const SQ2O2(::Type{BigFloat}) = big"0.707106781186547524400844362104849039284835937688474036588339868995366239231051"
 
 const PIO4(::Type{Float64}) = .78539816339744830962
 const PIO2(::Type{Float64}) = 1.5707963267948966
@@ -14,7 +15,9 @@ const TWOOPI(::Type{Float64}) = 0.6366197723675814
 const SQPIO2(::Type{Float64}) = 1.25331413731550025
 const SQ1O2PI(::Type{Float64}) = 0.3989422804014327
 const SQ2PI(::Type{Float64}) = 2.5066282746310007
+const SQ2O2(::Type{Float64}) = 0.7071067811865476
 
 const PIO4(::Type{Float32}) = 0.78539816339744830962f0
 const TWOOPI(::Type{Float32}) = 0.636619772367581343075535f0
 const THPIO4(::Type{Float32}) = 2.35619449019234492885f0
+const SQ2O2(::Type{Float32}) = 0.7071067811865476f0