add docs and constants

Author: heltonmc

src/Float128/besselj.jl

@@ -1,4 +1,4 @@
-function besselj0(x::BigFloat)
+#=function besselj0(x::BigFloat)
     x = abs(x)
     T = eltype(x)
     if iszero(x)
@@ -67,3 +67,4 @@ function besselj0(x::BigFloat)
         return z
     end
 end
+=#

src/besselj.jl

@@ -1,28 +1,48 @@
-#=
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1987, 2000 by Stephen L. Moshier
-https://github.com/jeremybarnes/cephes/blob/master/bessel/j0.c
-https://github.com/jeremybarnes/cephes/blob/master/bessel/j1.c
-=#
-function besselj0(x::Float64)
-    T = Float64
+#    Bessel functions of the first kind of order zero and one
+#                       besselj0, besselj1
+#
+#    Calculation of besselj0 is done in three branches using polynomial approximations
+#
+#    Branch 1: x <= 5.0
+#              besselj0 = (x^2 - r1^2)*(x^2 - r2^2)*P3(x^2) / Q8(x^2)
+#    where r1 and r2 are zeros of J0
+#    and P3 and Q8 are a 3 and 8 degree polynomial respectively
+#    Polynomial coefficients are from [1] which is based on [2]
+#    See [1] for more details and [2] for coefficients of polynomials
+#
+#    Branch 2: 5.0 < x < 80.0
+#              besselj0 = sqrt(2/(pi*x))*(cos(x - pi/4)*R7(x) - sin(x - pi/4)*R8(x))
+#    Hankel's asymptotic expansion is used
+#    where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
+#    See section 4 of [3] for more details and [1] for coefficients of polynomials
+# 
+#   Branch 3: x >= 80.0
+#              besselj0 = sqrt(2/(pi*x))*beta(x)*(cos(x - pi/4 - alpha(x))
+#   See modified expansions given in [3]. Exact coefficients are used
+#
+#   Calculation of besselj1 is done in a similar way as besselj0.
+#   See [3] for details on similarities.
+# 
+# [1] https://github.com/deepmind/torch-cephes
+# [2] Cephes Math Library Release 2.8:  June, 2000 by Stephen L. Moshier
+# [3] Harrison, John. "Fast and accurate Bessel function computation." 
+#     2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
+
+function besselj0(x::T) where T
     x = abs(x)
-    iszero(x) && return one(x)
     isinf(x) && return zero(x)
 
     if x <= 5
         z = x * x
         if x < 1.0e-5
             return 1.0 - z / 4.0
         end
-
-        DR1 = 5.78318596294678452118E0
-        DR2 = 3.04712623436620863991E1
-
+        DR1 = 5.78318596294678452118e0
+        DR2 = 3.04712623436620863991e1
         p = (z - DR1) * (z - DR2)
         p = p * evalpoly(z, RP_j0(T)) / evalpoly(z, RQ_j0(T))
         return p
-    elseif x < 100.0
+    elseif x < 80.0
         w = 5.0 / x
         q = 25.0 / (x * x)
 
@@ -40,8 +60,9 @@ function besselj0(x::Float64)
         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 = cos(x + xn)
+        xn = fma(xinv, evalpoly(x2, q), - PIO4(T))
+        b = cos(x)*cos(xn) - sin(x)*sin(xn)#cos(x + xn)
+        #b = cos(x + xn)
         return a * b
     end
 end
@@ -187,53 +208,3 @@ function besselj(n::Int, x::Float64)
 
     return sign * ans
 end
-
-
-function b2(v, x)
-
-    a = (x/2)^v
-    out = 1 / gamma(v + 1)
-    for k in 1:100
-        out += (-x^2 / 4)^k / factorial(BigFloat(k)) / factorial(v + BigFloat(k))
-    end
-    return out*a
-end
-
-function b(nu, x::T) where T
-    coef =  (x/2)^nu / gamma(nu + 1)
-
-    x2 = (x/2)^2
-
-
-    maxiter = 10000
-    b = one(T)
-    for k in 1:maxiter
-        z = -x2 / (k * (k + nu))
-        b += z*b
-    end
-    return b * coef
-end
-
-
-### for this function dd can be calculated accurately but there is a build up in error for val as we sum 
-## looks like only works for small arguments and smallish orders.... 
-function b3(v, x::T) where T
-    MaxIter = 5000
-    if v < 20
-        coef = (x / 2)^v / factorial(v)
-    else
-        vinv = inv(v)
-        coef = sqrt(vinv / (2 * π)) * MathConstants.e^(v * (log(x / (2 * v)) + 1)) 
-        coef *= evalpoly(vinv, (1, -1/12, 1/288,  139/51840, -571/2488320, -163879/209018880, 5246819/75246796800, 534703531/902961561600))
-    end
-
-    x2 = (x / 2)^2
-    val = zero(T)
-
-    for k in 1:MaxIter
-        val += coef
-        coef = -coef * x2 / (k * (v + k))
-        abs(coef) < eps(T) * abs(val) && break
-    end
-    return val
-end

src/math_constants.jl

@@ -2,7 +2,8 @@ const ONEOSQPI(::Type{BigFloat}) = big"5.641895835477562869480794515607725858440
 const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1"
 #const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924"
 #const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008"
-
+const SQ2OPI(::Type{BigFloat}) = big"0.7978845608028653558798921198687637369517172623298693153318516593413158517986017"
+const PIO4(::Type{BigFloat}) = big"0.7853981633974483096156608458198757210492923498437764552437361480769541015715495"
 
 const PIO4(::Type{Float64}) = .78539816339744830962
 const THPIO4(::Type{Float64}) = 2.35619449019234492885