Skip to content

Commit 51cb088

Browse files
committed
extend to 26 (gives better matching)
1 parent 73179d2 commit 51cb088

File tree

1 file changed

+5
-3
lines changed

1 file changed

+5
-3
lines changed

src/besselj.jl

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,8 @@ const PTS = (( 2.404825557695773 , -1.176691651530894e-16),
4545
(19.615858510468243 , -1.004445634526616e-15),
4646
(21.21163662987926 , 4.947077428784068e-16),
4747
(22.760084380592772 , -4.925749373614922e-16),
48-
(24.352471530749302 , 9.169067133951066e-16))
48+
(24.352471530749302 , 9.169067133951066e-16),
49+
(25.903672087618382 , 4.894530726419825e-16))
4950
const POLYS = ((0.0, -0.5191474972894669, 0.10793870175491979, 0.05660177443794914, -0.008657669593292222, -0.0021942003590739974, 0.0002643770365942964, 4.37291931443113e-5, -4.338825419759815e-6, -5.304927679598406e-7, 4.469819175606667e-8, 4.3284827345621115e-9, -3.135111000732148e-10, -2.628876489517534e-11),
5051
(-0.402759395702553, 2.476919088072758e-16, 0.20137969785127532, -0.017518715285670765, -0.013352611033152278, 0.0010359438492839443, 0.00037218755624680334, -2.4952042421113972e-5, -5.776086353844158e-6, 3.374317129436161e-7, 5.727482259215452e-8, -2.9561880489355444e-9, -3.905845672635605e-10, 1.971332566705736e-11),
5152
(0.0, 0.34026480655836816, -0.030820651425593214, -0.05298855286760721, 0.004631042145890388, 0.002257440229081131, -0.00017518572879518415, -4.6521091062878115e-5, 3.199785909661533e-6, 5.716500268232186e-7, -3.5112898510841466e-8, -4.684643389757727e-9, 2.562685034682206e-10, 2.7958958795750104e-11),
@@ -60,7 +61,8 @@ const POLYS = ((0.0, -0.5191474972894669, 0.10793870175491979, 0.056601774437949
6061
(0.18006337534431555, -5.638484737644332e-17, -0.09003168767215539, 0.0015299132863046024, 0.0074441453680234426, -0.00015060569680378768, -0.0002443398416353605, 5.227001519013193e-6, 4.267152607972633e-6, -9.295966007495808e-8, -4.610438011417262e-8, 1.0049092632275165e-9, 3.339442682325105e-10, -7.4991079301099e-12),
6162
(0.0, -0.17326589422922986, 0.004084217951979124, 0.028749284970146657, -0.0006761643016121907, -0.0014215899173758441, 3.320978125391802e-5, 3.3264379323815026e-5, -7.684444100376941e-7, -4.515479818833484e-7, 1.0271599456634145e-8, 3.993903280527723e-9, -8.793975528824604e-11, -2.4610374225652004e-11),
6263
(-0.16718460047381803, 4.662597138876655e-17, 0.08359230023690671, -0.0012242529339116864, -0.006925682915280748, 0.00012100729852042988, 0.00022821854128498174, -4.230799709959437e-6, -4.007618360337058e-6, 7.601217108010105e-8, 4.358483157250522e-8, -8.317621452517008e-10, -3.178805429572483e-10, 6.284538399593043e-12),
63-
(0.0, 0.16170155068925002, -0.003320023400603752, -0.026859370386562005, 0.0005505380905962079, 0.001331699465921981, -2.7156832110124053e-5, -3.128954487528851e-5, 6.326902787990595e-7, 4.269717284901734e-7, -8.532904814909196e-9, -3.799560194537855e-9, 7.410932094856478e-11, 2.3903214947728713e-11))
64+
(0.0, 0.16170155068925002, -0.0033200234006036146, -0.02685937038656159, 0.0005505380905909195, 0.0013316994659124243, -2.715683205388875e-5, -3.1289544794362e-5, 6.326900360777323e-7, 4.2697141781883964e-7, -8.532447325590315e-9, -3.799009445641295e-9, 7.379552811590934e-11, 2.353654960834717e-11),
65+
(0.15672498625285222, 1.1464880342445208e-19, -0.07836249312642612, 0.0010083833270351796, 0.006501011610557426, -9.993664895655577e-5, -0.00021478298462967253, 3.511086890959515e-6, 3.7857022791122945e-6, -6.350944888510364e-8, -4.135588012575766e-8, 7.135988717747828e-10, 3.2075281131621564e-10, 2.208506585533542e-12))
6466

6567
"""
6668
besselj0(x::T) where T <: Union{Float32, Float64}
@@ -72,7 +74,7 @@ function besselj0(x::Float64)
7274
x = abs(x)
7375
isinf(x) && return zero(x)
7476

75-
if x < 25.0
77+
if x < 26.0
7678
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))
7779
n = round(Int, fma(2/pi,x,-1/2))
7880
root = @inbounds PTS[n]

0 commit comments

Comments
 (0)