increase range for large argument approx

Author: heltonmc

src/U_polynomials.jl

@@ -43,7 +43,7 @@ function Uk_poly_In(p, v, p2, ::Type{Float32})
     u13 = -1 / 955134445051714068660879360000 * evalpoly(p2, (17438611142828905996129258798828125, -3698121486504259988094897605296209375, 136735019134677724428035696765082813750, -2069933923586966756183324291232117362250, 16843538631795795357786827345307534156063, -83924867223075156862785921508524155665245, 274983827478138958623041508409195988431140, -616410216242554698436702353237166008656700, 962926533925253374359704288384340809260875, -1049095945162229046324321461816274931715625, 782463969315283937856703223178540650343750, -381190503845282445953419057314665534156250, 109361210755577700442544717509565392265625, -14020668045586884672121117629431460546875))
     u14 = 1 / 45846453362482275295722209280000 * evalpoly(p2, (5448320367052402487647812713291015625, -1338184074771428116079233795614103631250, 57170953417612444837142230812990944671875, -1000503839668383458999731491914516564625300, 9440449669103391509091075981237243128469201, -54857705817658080981995319669299096598482382, 211477117385619365164298957821904115470535555, -564850830044980230366140582063618983657685400, 1070439683260179398514645952824098811025619475, -1451823699927947602004297385351260623500372750, 1401302601668131482630653233972052729323190625, -940627071986145750405920450097257805227812500, 417630985812040040477569028872405152769921875, -110320224449895843354117801955418504167031250, 13133360053559028970728309756597440972265625))   
     u15 = -1 / 9820310310243703368343697227776000000 * evalpoly(p2, (8178936810213560828419581728001773291015625, -2303431987527333128955769182299845911305390625, 112597271053778779753048514469995937998172890625, -2254933495791765108580529087615802250458013685625, 24403480234538299231733883413666768614198435948125, -163359140754958502896104062604202934925448173291477, 730367145705123976114617970888707594104468177381925, -2284251621937242886581917667066122422330060024456125, 5136561256208409671660362298619778869859994724706875, -8420533422834140468835467666391400380550043688871875, 10085018700249896522602267572484630409640764997271875, -8735135969643867540297524795790262235822823374296875, 5329871927856528282113994744458999865006055974609375, -2173722139119126857509156976742601742357422740234375, 532039967451707060045861997017872377315039521484375, -59115551939078562227317999668652486368337724609375))
-    @show evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
+
     return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
 end
 
@@ -64,15 +64,13 @@ function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
     u13 = -evalpoly(p2, (18257.755474293175, -3.8718334425726123e6, 1.43157876718889e8, -2.167164983223795e9, 1.763473060683497e10, -8.786707217802327e10, 2.879006499061506e11, -6.453648692453765e11, 1.0081581068653821e12, -1.0983751560812233e12, 8.192186695485773e11, -3.990961752244665e11, 1.144982377320258e11, -1.4679261247695616e10))
     u14 = evalpoly(p2, (118838.42625678326, -2.9188388122220814e7, 1.2470092935127103e9, -2.1822927757529224e10, 2.0591450323241e11, -1.1965528801961816e12, 4.612725780849132e12, -1.2320491305598287e13, 2.334836404458184e13, -3.166708858478516e13, 3.056512551993532e13, -2.0516899410934438e13, 9.109341185239898e12, -2.406297900028504e12, 2.86464035717679e11))   
     u15 = -evalpoly(p2, (832859.3040162893, -2.3455796352225152e8, 1.1465754899448236e10, -2.2961937296824646e11, 2.4850009280340854e12, -1.663482472489248e13, 7.437312290867914e13, -2.3260483118893994e14, 5.230548825784446e14, -8.57461032982895e14, 1.0269551960827625e15, -8.894969398810265e14, 5.4273966498765975e14, -2.213496387025252e14, 5.417751075510605e13, -6.019723417234006e12))
-    @show evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
-
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
+    u16 = evalpoly(p2, (6.252951493434797e6, -2.0016469281917763e9, 1.1099740513917902e11, -2.5215584749128545e12, 3.100743647289646e13, -2.3665253045164925e14, 1.2126758042503475e15, -4.3793258383640155e15, 1.1486706978449752e16, -2.2268225133911144e16, 3.213827526858624e16, -3.4447226006485144e16, 2.705471130619708e16, -1.5129826322457682e16, 5.705782159023671e15, -1.3010127235496995e15, 1.3552215870309369e14))
+    u17 = -evalpoly(p2, (5.0069589531988926e7, -1.8078220384658062e10, 1.128709145410874e12, -2.886383763141476e13, 4.0004445704303625e14, -3.4503855118462725e15, 2.0064271476309532e16, -8.270945651585064e16, 2.4960365126160426e17, -5.62631788074636e17, 9.575335098169139e17, -1.2336116931960694e18, 1.1961991142756308e18, -8.592577980317548e17, 4.4347954614171904e17, -1.5552983504313904e17, 3.3192764720355224e16, -3.254192619642669e15))
+    
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17))
 end
 
 
-
-    
-   
 #=
 function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat
     u0 = one(T)

src/besselj.jl

@@ -174,7 +174,7 @@ function _besselj(nu, x)
             jnup1 = besselj_debye(v+1, x)
             return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
         end
-    elseif x > 1.85*nu
+    elseif x > 1.65*nu
         return besselj_large_argument(nu, x)
     end
     if nu <= 100
@@ -269,24 +269,14 @@ function besselj_debye(v, x)
     b = sqrt(vmx)
 
     n = v * log(vdx + sqrt(vdx^2 - 1)) - b
+
     coef = SQ1O2PI(S) * exp(-n) / sqrt(b)
     p = v / b
     p2  = v^2 / vmx
 
     return coef * Uk_poly_Jn(p, v, p2, T)
 end
 
-function besselj_debye_large_order(v, x)
-    T = eltype(x)
-    α = asech(x/v)
-    coef = exp(v*(tanh(α) - α))
-    coef /= sqrt(2*T(pi)*v*tanh(α))
-    p = coth(α)
-    p2 = p^2
-
-    return coef * Uk_poly_Jn(p, v, p2, T)
-end
-
 # For 0.0 <= x < 171.5
 # Mean ULP = 0.55
 # Max ULP = 2.4

test/besselj_test.jl

@@ -103,7 +103,7 @@ end
 x = 0.1:3.0:800.0
 nu = 2:4:500
 for v in nu, xx in x
-    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(BigFloat(v), BigFloat(xx)), rtol=1e-14)
+    @test isapprox(Bessels._besselj(BigFloat(v), BigFloat(xx)), SpecialFunctions.besselj(BigFloat(v), BigFloat(xx)), rtol=5e-13)
 end
 
 # test half orders (SpecialFunctions does not give big float precision)