fix downard recurrence and large orders

Author: heltonmc

src/U_polynomials.jl

@@ -48,7 +48,13 @@ function Uk_poly_In(p, v, p2, ::Type{Float32})
     u18 = evalpoly(p2, (4.259392165047669e8, -1.722832387173505e11, 1.2030115826419191e13, -3.4396530474307594e14, 5.335106978708839e15, -5.1605093193485224e16, 3.37667624979061e17, -1.5736434765189599e18, 5.402894876715982e18, -1.3970803516443374e19, 2.757282981650519e19, -4.178861444656839e19, 4.859942729324836e19, -4.301555703831444e19, 2.846521225167657e19, -1.3639420410571592e19, 4.47020096401231e18, -8.966114215270463e17, 8.30195760673191e16))
     u19 = -evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
     u20 = evalpoly(p2, (3.646840080706556e10, -1.818726203851104e13, 1.5613123930484672e15, -5.48403360388329e16, 1.0461721131134344e18, -1.2483700995047234e19, 1.0126774169536592e20, -5.8917941350694964e20, 2.548961114664972e21, -8.405915817108351e21, 2.1487414815055883e22, -4.302534303482379e22, 6.783661642951883e22, -8.423222750084323e22, 8.19433100543513e22, -6.173206302884415e22, 3.528435843903409e22, -1.4787743528433614e22, 4.285296082829494e21, -7.671943936729004e20, 6.393286613940837e19))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20))
+    u21 = -evalpoly(p2, (3.6490108188498334e11, -2.0052440123627112e14, 1.894406984252143e16, -7.319501491566134e17, 1.5365025218443373e19, -2.0197335419300872e20, 1.8081594057131945e21, -1.1640246461465369e22, 5.591591380366263e22, -2.0566149136271542e23, 5.8965434619782445e23, -1.3337178907798302e24, 2.3967237744351682e24, -3.430872898515746e24, 3.905264103536985e24, -3.511096528332644e24, 2.461506085403875e24, -1.3170969618092387e24, 5.194289094766812e23, -1.4228394823321413e23, 2.417461500896379e22, -1.91862023880665e21))
+    u22 = evalpoly(p2, (3.8335346613939443e12, -2.3109159761323565e15, 2.3920280120269997e17, -1.0121818379942089e19, 2.3275346258089414e20, -3.3544689122226785e21, 3.297557757461478e22, -2.336107524486965e23, 1.238524103792452e24, -5.0463598652544e24, 1.6103128541137314e25, -4.077501349206541e25, 8.26258535798955e25, -1.3459193994556415e26, 1.7635713272326644e26, -1.8526731041549917e26, 1.548092083577385e26, -1.0148048982766395e26, 5.103920268388802e25, -1.9006807535664433e25, 4.936185283790662e24, -7.980021228256559e23, 6.04547062746709e22))
+    u23 = -evalpoly(p2, (4.218971570284097e13, -2.778481101311081e16, 3.1385283211499996e18, -1.4486387749510863e20, 3.6341499869780876e21, -5.7179919065432055e22, 6.144339925144987e23, -4.766924608251481e24, 2.774466490672939e25, -1.2449342046124282e26, 4.392130563430048e26, -1.2355529146787609e27, 2.7982068996977173e27, -5.131998439010333e27, 7.641216535678268e27, -9.228395023257356e27, 8.999255845917453e27, -7.02322235515725e27, 4.322773732100187e27, -2.050902994929233e27, 7.234243234844319e26, -1.7860680966743495e26, 2.753863007576946e25, -1.9955529040412654e24))
+    u24 = evalpoly(p2, (4.8540146868529006e14, -3.4792991439250445e17, 4.273207395701127e19, -2.1435653415108537e21, 5.844687629283339e22, -1.0000750138961727e24, 1.1699189691874474e25, -9.896648661695488e25, 6.29370256208713e26, -3.0939194683063286e27, 1.1998211967644424e28, -3.7252346341093444e28, 9.358117764887965e28, -1.9153963148099324e29, 3.206650343980748e29, -4.395132918078325e29, 4.9215508698387624e29, -4.4775348387950634e29, 3.277658265637452e29, -1.9012207767547338e29, 8.536184882279286e28, -2.8599776383548e28, 6.728957650918171e27, -9.916401268407057e26, 6.886389769727123e25))
+    u25 = -evalpoly(p2, (5.827244631566907e15, -4.5305357275125955e18, 6.029638127487473e20, -3.2761234100445222e22, 9.675654883193622e23, -1.7941040647617987e25, 2.2764310713849358e26, -2.0914533474677945e27, 1.4471195817119858e28, -7.757785573404132e28, 3.2900927159291354e29, -1.1210232552135908e30, 3.1034661143911036e30, -7.036055338636485e30, 1.3128796688902614e31, -2.0208792587851872e31, 2.5653099826522344e31, -2.6771355605594045e31, 2.2823085118856488e31, -1.5730388076301427e31, 8.627355824571355e30, -3.676221426681414e30, 1.1728484268744769e30, -2.6355294419807464e29, 3.7195112743738626e28, -2.479674182915908e27))
+
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20, u21, u22, u23, u24, u25))
 end
 
 function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
@@ -73,9 +79,14 @@ function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
     u18 = evalpoly(p2, (4.259392165047669e8, -1.722832387173505e11, 1.2030115826419191e13, -3.4396530474307594e14, 5.335106978708839e15, -5.1605093193485224e16, 3.37667624979061e17, -1.5736434765189599e18, 5.402894876715982e18, -1.3970803516443374e19, 2.757282981650519e19, -4.178861444656839e19, 4.859942729324836e19, -4.301555703831444e19, 2.846521225167657e19, -1.3639420410571592e19, 4.47020096401231e18, -8.966114215270463e17, 8.30195760673191e16))
     u19 = -evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
     u20 = evalpoly(p2, (3.646840080706556e10, -1.818726203851104e13, 1.5613123930484672e15, -5.48403360388329e16, 1.0461721131134344e18, -1.2483700995047234e19, 1.0126774169536592e20, -5.8917941350694964e20, 2.548961114664972e21, -8.405915817108351e21, 2.1487414815055883e22, -4.302534303482379e22, 6.783661642951883e22, -8.423222750084323e22, 8.19433100543513e22, -6.173206302884415e22, 3.528435843903409e22, -1.4787743528433614e22, 4.285296082829494e21, -7.671943936729004e20, 6.393286613940837e19))
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20))
-end
+    u21 = -evalpoly(p2, (3.6490108188498334e11, -2.0052440123627112e14, 1.894406984252143e16, -7.319501491566134e17, 1.5365025218443373e19, -2.0197335419300872e20, 1.8081594057131945e21, -1.1640246461465369e22, 5.591591380366263e22, -2.0566149136271542e23, 5.8965434619782445e23, -1.3337178907798302e24, 2.3967237744351682e24, -3.430872898515746e24, 3.905264103536985e24, -3.511096528332644e24, 2.461506085403875e24, -1.3170969618092387e24, 5.194289094766812e23, -1.4228394823321413e23, 2.417461500896379e22, -1.91862023880665e21))
+    u22 = evalpoly(p2, (3.8335346613939443e12, -2.3109159761323565e15, 2.3920280120269997e17, -1.0121818379942089e19, 2.3275346258089414e20, -3.3544689122226785e21, 3.297557757461478e22, -2.336107524486965e23, 1.238524103792452e24, -5.0463598652544e24, 1.6103128541137314e25, -4.077501349206541e25, 8.26258535798955e25, -1.3459193994556415e26, 1.7635713272326644e26, -1.8526731041549917e26, 1.548092083577385e26, -1.0148048982766395e26, 5.103920268388802e25, -1.9006807535664433e25, 4.936185283790662e24, -7.980021228256559e23, 6.04547062746709e22))
+    u23 = -evalpoly(p2, (4.218971570284097e13, -2.778481101311081e16, 3.1385283211499996e18, -1.4486387749510863e20, 3.6341499869780876e21, -5.7179919065432055e22, 6.144339925144987e23, -4.766924608251481e24, 2.774466490672939e25, -1.2449342046124282e26, 4.392130563430048e26, -1.2355529146787609e27, 2.7982068996977173e27, -5.131998439010333e27, 7.641216535678268e27, -9.228395023257356e27, 8.999255845917453e27, -7.02322235515725e27, 4.322773732100187e27, -2.050902994929233e27, 7.234243234844319e26, -1.7860680966743495e26, 2.753863007576946e25, -1.9955529040412654e24))
+    u24 = evalpoly(p2, (4.8540146868529006e14, -3.4792991439250445e17, 4.273207395701127e19, -2.1435653415108537e21, 5.844687629283339e22, -1.0000750138961727e24, 1.1699189691874474e25, -9.896648661695488e25, 6.29370256208713e26, -3.0939194683063286e27, 1.1998211967644424e28, -3.7252346341093444e28, 9.358117764887965e28, -1.9153963148099324e29, 3.206650343980748e29, -4.395132918078325e29, 4.9215508698387624e29, -4.4775348387950634e29, 3.277658265637452e29, -1.9012207767547338e29, 8.536184882279286e28, -2.8599776383548e28, 6.728957650918171e27, -9.916401268407057e26, 6.886389769727123e25))
+    u25 = -evalpoly(p2, (5.827244631566907e15, -4.5305357275125955e18, 6.029638127487473e20, -3.2761234100445222e22, 9.675654883193622e23, -1.7941040647617987e25, 2.2764310713849358e26, -2.0914533474677945e27, 1.4471195817119858e28, -7.757785573404132e28, 3.2900927159291354e29, -1.1210232552135908e30, 3.1034661143911036e30, -7.036055338636485e30, 1.3128796688902614e31, -2.0208792587851872e31, 2.5653099826522344e31, -2.6771355605594045e31, 2.2823085118856488e31, -1.5730388076301427e31, 8.627355824571355e30, -3.676221426681414e30, 1.1728484268744769e30, -2.6355294419807464e29, 3.7195112743738626e28, -2.479674182915908e27))
 
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20, u21))
+end
 #=
 function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat
     u0 = one(T)

src/besselj.jl

@@ -170,11 +170,12 @@ function _besselj(nu, x)
     nu > debye_cutoff && return besselj_debye(nu, x)
 
     if nu >= x
-        nu_shift = ceil(7.0 + 1.00033*x + (1427.61*x)^(1/3))
+        nu_shift = ceil(Int, 7.2 + 1.00033*x + (1427.61*x)^(1/3) - nu)
         v = nu + nu_shift
+        arr = range(v, stop = nu, length = nu_shift + 1)
         jnu = besselj_debye(v, x)
         jnup1 = besselj_debye(v+1, x)
-        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
+        return besselj_down_recurrence(x, jnu, jnup1, arr)[2]
     end
 
     # at this point x > nu and  x < nu * 1.65
@@ -193,11 +194,12 @@ function _besselj(nu, x)
         jnum1 = besselj_large_argument(v2 - 1, x)
         return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
     else
-        nu_shift = ceil(debye_diff)
+        nu_shift = ceil(Int, debye_diff)
         v = nu + nu_shift
+        arr = range(v, stop = nu, length = nu_shift + 1)
         jnu = besselj_debye(v, x)
         jnup1 = besselj_debye(v+1, x)
-        return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
+        return besselj_down_recurrence(x, jnu, jnup1, arr)[2]
     end
 end
 
@@ -257,7 +259,7 @@ function besselj_debye(v, x)
     S = promote_type(T, Float64)
     x = S(x)
 
-    vmx = (v + x) * (v − x)
+    vmx = (v + x) * (v - x)
     vs = sqrt(vmx)
     n  = fma(v, -log(x / (v + vs)), -vs)
 

src/recurrence.jl

@@ -45,14 +45,18 @@ end
     end
     return jnup1, jnum1
 end
-@inline function besselj_down_recurrence(x, jnu, jnup1, nu_start, nu_end)
+@inline function besselj_down_recurrence(x, jnu, jnup1, arr)
+    # arr is the index of Bessel orders arr = nu_start:-1:nu_end
+    # but needs special care if nu is a decimal
+    # use  v = nu + nu_shift
+    # arr = range(v, stop = nu, length = nu_shift + 1)
     jnum1 = jnup1
     x2 = 2 / x
-    for n in nu_start:-1:nu_end
+    for n in arr
         a = x2 * n
         jnum1 = muladd(a, jnu, -jnup1)
         jnup1 = jnu
         jnu = jnum1
     end
     return jnum1, jnup1
-end
+end

test/besselj_test.jl

@@ -71,49 +71,36 @@ j1_32 = besselj1.(Float32.(x))
 @test besselj1(-80.0) ≈ SpecialFunctions.besselj1(-80.0)
 
 ## Tests for besselj 
-# note this is not complete just a simple test
-# this needs work and removing for now
 
-#@test besselj(3, 1.0) ≈ SpecialFunctions.besselj(3, 1.0)
-#@test besselj(-5, 6.1) ≈ SpecialFunctions.besselj(-5, 6.1)
-
-## test the small value approximation using power series 
-nu = [0.5, 1.5, 3.0, 10.0, 22.2, 35.0, 52.1, 100.0, 250.2, 500.0, 1000.0]
-x = [0.001, 1.0, 5.0, 15.0, 19.9]
-
-for v in nu, x in x
-    @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
-end
-
-## test the large argument asymptotic expansion
-vnu = ((0.5, 21.0), (0.5, 45.0), (3.5, 21.0), (10.0, 21.0), (21.2, 45.0), (43.2, 90.0), (100.1, 205.2))
-
-for (v, x) in vnu
-    @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
-end
-
-## test the debye uniform asymptotic expansion for x < nu
-vnu = ((60.5, 21.0), (100.5, 45.0), (150.5, 61.0))
-
-for (v, x) in vnu
-    @test Bessels._besselj(v, x) ≈ SpecialFunctions.besselj(v, x)
-end
+#=
+Notes
+    - power series shows larger error when nu is large (146) and x is small (1.46)
+    - asymptotic expansion shows larger error when nu is large or x is large
+=#
 
 ## test all numbers and orders for 0<nu<100
-x = 0.1:3.0:800.0
-nu = 2:4:500
+x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05]
+nu = [2, 4, 6, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100]
 for v in nu, xx in x
-    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(BigFloat(v), BigFloat(xx)), rtol=5e-13)
+    xx *= v
+    sf = SpecialFunctions.besselj(BigFloat(v), BigFloat(xx))
+    @test isapprox(Bessels._besselj(v, xx), Float64(sf), rtol=5e-14)
 end
 
 # test half orders (SpecialFunctions does not give big float precision)
-# The SpecialFunctions implementation is actually very inaccurate
-# julia> 1 - SpecialFunctions.besselj(10.5, 29.6) / -0.009263478934797420709865
-#-1.63202784619898e-13
-# with value from https://keisan.casio.com/exec/system/1180573474
-x = 0.1:0.5:100.0
-nu = 2.5:1.0:10.5
+# The SpecialFunctions implementation is only accurate to about 1e-11 - 1e-13
+
+x = [0.05, 0.1, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92, 0.95, 0.97, 0.99, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2, 1.4, 1.5, 1.6, 1.8, 2.0, 2.5, 3.0]
+nu = [0.1, 0.4567, 0.8123, 1.5, 2.5, 4.1234, 6.8, 12.3, 18.9, 28.2345, 38.1235, 51.23, 72.23435, 80.5, 98.5, 104.2]
 for v in nu, xx in x
-    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(v, xx), rtol=1e-12)
+    xx *= v
+    @test isapprox(Bessels._besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=1e-12)
 end
 
+## test large orders
+nu = [150, 165.2, 200.0, 300.0, 500.0, 1000.0, 5000.2, 10000.0, 50000.0]
+x = [0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8, 0.85, 0.9, 0.92,0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,0.995, 0.999, 1.0, 1.01, 1.05, 1.08, 1.1, 1.2]
+for v in nu, xx in x
+    xx *= v
+    @test isapprox(Bessels._besselj(v, xx), SpecialFunctions.besselj(v, xx), rtol=5e-11)
+end