fix behavior at 0

Author: oscardssmith

src/besselj.jl

@@ -73,7 +73,7 @@ function besselj0(x::Float64)
     isinf(x) && return zero(x)
 
     if x < 25.0
-        x<pi/2 && return evalpoly(x*x, (0.9999999999999997, -0.24999999999998027, 0.015624999999831789, -0.00043402777723176535, 6.78168315716265e-6, -6.781608601544183e-8, 4.705878468811015e-10, -2.311849729995118e-12))
+        x<pi/2 && return evalpoly(x*x, (1.0, -0.25, 0.01562499999999994, -0.00043402777777725544, 6.781684026082576e-6, -6.781683757550061e-8, 4.709479394601058e-10, -2.4016837144506874e-12, 9.104258208703104e-15))
         n = round(Int, fma(2/pi,x,-1/2))
         root = @inbounds PTS[n]
         r = x - root[1] - root[2]