update evalpoly routine

Author: heltonmc

src/besselj.jl

@@ -151,41 +151,72 @@ function besselj1(x::Float32)
     end
 end
 
-function besselj_large_argument(v, x::T) where T
+function _α_αp_asymptotic(v, x::T) where T
     μ = 4 * v^2
     s0 = 1
     s1 = (1 - μ) / 8
-    s2 = (-μ^2 + 26μ - 25) / 128
-    s3 = (-μ^3 + 115μ^2 - 1187μ + 1073) / 1024
-    s4 = (-5μ^4 + 1540μ^3 - 56238μ^2 + 430436μ - 375733) / 32768
-    s5 = (-7μ^5 + 4515μ^4 - 397190μ^3 + 9716998μ^2 - 64709091μ + 55384775) / 262144
-    s6 = (-21μ^6 + 24486μ^5 - 4238531μ^4 + 235370036μ^3 - 4733751627μ^2 + 29215626566μ - 24713030909) / 4194304
-    s7 = (-33μ^7 + 63063μ^6 - 18939349μ^5 + 1989535379μ^4 - 87480924739μ^3 + 1573356635461μ^2 - 9268603618823μ + 7780757249041) / 33554432
-    s8 = (-429μ^8 + 1252680μ^7 - 598859404μ^6 + 106122595896μ^5 - 8593140373614μ^4 + 329343318168440μ^3 - 5517359285625804μ^2 + 31505470994964360μ - 26308967412122125) / 2147483648
-    s9 = (-715μ^9 + 3026595μ^8 - 2163210764μ^7 + 597489288988μ^6 - 79979851361130μ^5 + 5536596631240042μ^4 - 194026764558396188μ^3 + 3095397360215483916μ^2 - 17285630127691126691μ + 14378802319925055947) / 17179869184
-    s10 = (-2431μ^10 + 14318590μ^9 - 14587179035μ^8 + 5925778483368μ^7 - 1216961874423502μ^6 + 137164402798287604μ^5 - 8556293060689145118μ^4 + 281576004385356401192μ^3 - 4334421752432088249971μ^2 + 23801928703130666089534μ - 19740662615375374580231) / 274877906944
-    s11 = (-4199μ^11 + 33302269μ^10 - 46575001473μ^9 + 26616800349043μ^8 - 7938052610367302μ^7 + 1355921534379626242μ^6 - 136097258378143928354μ^5 + 7882145023568373727078μ^4 - 247560495216085669748963μ^3 + 3708133715826433118900337μ^2 - 20097626064642945009122253μ + 16629305448257355302267575) / 2199023255552
-    s12 = (-29393μ^12 + 305569628μ^11 - 569135407698μ^10 + 441734362333708μ^9 - 183408097342762543μ^8 + 45038004674954790968μ^7 - 6784356890740220944060μ^6 + 626292216963766324300664μ^5 - 34303895772550410316410943μ^4 + 1039206680182498322210164588μ^3 - 15233583205322702897380357650μ^2 + 81699279492428006742420030332μ - 67471218624230362526181277601) / 70368744177664
-
-
-    αp = s0 + s1/x^2 + s2/x^4 + s3/x^6 + s4/x^8 + s5/x^10 + s6/x^12 + s7/x^14 + s8/x^16 + s9/x^18 + s10/x^20 + s11/x^22 + s12/x^24
-    α = s0*x - s1/x - s2/3x^3 - s3/5x^5 - s4/7x^7 - s5/9x^9 - s6/11x^11 - s7/13x^13 - s8/15x^15 - s9/17x^17 - s10/19x^19 - s11/21x^21 - s12/23x^23
-
+    s2 = evalpoly(μ, (-25, 26, -1)) / 128
+    s3 = evalpoly(μ, (1073, -1187, 115, -1)) / 1024
+    s4 = evalpoly(μ, (-375733, 430436, -56238, 1540, -5)) / 32768
+    s5 = evalpoly(μ, (55384775, -64709091, 9716998, -397190, 4515, -7)) / 262144
+    s6 = evalpoly(μ, (-24713030909, 29215626566, -4733751627, 235370036, -4238531, 24486, -21)) / 4194304
+
+    if x < 15*v
+        s7 = evalpoly(μ, (7780757249041, -9268603618823, 1573356635461, -87480924739, 1989535379, -18939349, 63063, -33)) / 33554432
+        s8 = evalpoly(μ, (-26308967412122125, 31505470994964360, -5517359285625804, 329343318168440, -8593140373614, 106122595896, -598859404, 1252680, -429)) / 2147483648
+        s9 = evalpoly(μ, (14378802319925055947, -17285630127691126691, 3095397360215483916, -194026764558396188, 5536596631240042, -79979851361130, 597489288988, -2163210764, 3026595, -715)) / 17179869184
+        s10 = evalpoly(μ, (-19740662615375374580231, 23801928703130666089534, -4334421752432088249971, 281576004385356401192, -8556293060689145118, 137164402798287604, -1216961874423502, 5925778483368, -14587179035, 14318590, -2431)) / 274877906944
+        s11 = evalpoly(μ, (16629305448257355302267575, -20097626064642945009122253, 3708133715826433118900337, -247560495216085669748963, 7882145023568373727078, -136097258378143928354, 1355921534379626242, -7938052610367302, 26616800349043, -46575001473, 33302269, -4199)) / 2199023255552
+        s12 = evalpoly(μ, (-67471218624230362526181277601, 81699279492428006742420030332, -15233583205322702897380357650, 1039206680182498322210164588, -34303895772550410316410943, 626292216963766324300664, -6784356890740220944060, 45038004674954790968, -183408097342762543, 441734362333708, -569135407698, 305569628, -29393)) / 70368744177664
+        s13 = evalpoly(μ, (81121053093143918050102987600099,  -98383136980007985475571939027975, 18503851640550924693860888790250, -1284587641197372537980918233562, 43642612323123398148213493825, -832261743620394930038568309, 9606500048207046619311900, -69893489053616393418780, 325411350127146969525, -959952773579760385, 1717262634951770, -1676314996330, 692939975, -52003)) / 562949953421312
+
+        αp = evalpoly((1/x)^2, (s0, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, s13))
+        α = x * (evalpoly((1/x)^2, (s0, -s1, -s2/3, -s3/5, -s4/7, -s5/9, -s6/11, -s7/13, -s8/15, -s9/17, -s10/19, -s11/21, -s12/23, -s13/25)))
+    else
+        αp = evalpoly((1/x)^2, (s0, s1, s2, s3, s4, s5, s6))
+        α = x * (evalpoly((1/x)^2, (s0, -s1, -s2/3, -s3/5, -s4/7, -s5/9, -s6/11)))
+    end
     α = α - T(pi)/4 - T(pi)/2 * v
-
     b = sqrt(2 / T(pi)) / sqrt(αp * x)
+
+    return α, b
+end
+function besselj_large_argument(v, x::T) where T
+    α, b = _α_αp_asymptotic(v, x)
     return cos(α)*b
 end
+function bessely_large_argument(v, x::T) where T
+    α, b = _α_αp_asymptotic(v, x)
+    return sin(α)*b
+end
 
+
+# only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))
 function besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 200
+    MaxIter = 1_000
     out = zero(T)
     a = (x/2)^v / gamma(v + one(T))
     t2 = (x/2)^2
     for i in 0:MaxIter
         out += a
-        abs(a) < eps(T)*abs(out) && break
-        a = -a / (v + i + one(T)) * inv(i + one(T)) * t2
+        abs(a) < eps(T) * abs(out) && break
+        a = -a * inv((v + i + one(T)) * (i + one(T))) * t2
     end
     return out
 end
+
+# perhaps use when v is small i believe v also has to be positive for this to work
+# need for bessely 
+function log_besselj_small_arguments_orders(v, x::T) where T
+    MaxIter = 200
+    out = zero(T)
+    a   = one(T)
+    x2 = (x/2)^2
+    for i in 0:MaxIter
+        out += a
+        a = -a * x2 * inv((i + one(T)) * (v + i + one(T)))
+        #(abs(a) < eps(T) * abs(out)) && break
+    end
+    logout = -loggamma(v + 1) + v * log(x/2) + log(out)
+    return logout
+end