use taylor expansion for besseli1

heltonmc · heltonmc · commit 938bf686f5a6 · 2022-12-01T15:08:18.000-05:00
diff --git a/src/besseli.jl b/src/besseli.jl
@@ -212,7 +212,6 @@ function besseli0(z::ComplexF32)
         z = conj(z)
         isconj = true
     end
-
     if abs(z) > 10.0
         zinv = 1 / z
         e = exp(z)
@@ -243,64 +242,51 @@ function besseli1(z::ComplexF64)
     else
         isconj = false
     end
-
     if abs(z) > 17.5
+        # asymptotic expansion
         zinv = 1 / z
         e = exp(z)
         sinv = sqrt(zinv) * SQ1O2PI(Float64)
         p = evalpoly(zinv*zinv, (1.0, -0.1171875, -0.144195556640625, -0.6765925884246826, -6.883914268109947, -121.59789187653587, -3302.2722944808525, -127641.2726461746, -6.656367718817688e6, -4.502786003050393e8, -3.8338575207427895e10, -4.0118385991331978e12, -5.060568503314726e14, -7.572616461117957e16, -1.3262572853205555e19))
         p2 = evalpoly(zinv*zinv, (0.375, 0.1025390625, 0.2775764465332031, 1.993531733751297, 27.248827311268542, 603.8440767050702, 19718.37591223663, 890297.8767070678, 5.310411010968523e7, 4.043620325107754e9, 3.827011346598606e11, 4.406481417852279e13, 6.065091351222699e15, 9.83388387659068e17))
         r = e * sinv * muladd(-p2, zinv, p) - im * muladd(p2, zinv, p) * sinv / e
+    elseif real(z) < 4.5
+        # taylor series around the roots (slight offset) of J1(z)
+        # use relation I1(z) = J1(im*z) / im
+        _z = imag(z) + abs(real(z))*im
+        if real(_z) < 2.2
+            # power series for I0
+            r = z*evalpoly(z*z, (0.5, 0.0625, 0.0026041666666666665, 5.425347222222222e-5, 6.781684027777778e-7, 5.651403356481481e-9, 3.363930569334215e-11, 1.5017547184527747e-13, 5.214426105738801e-16, 1.4484516960385557e-18, 3.2919356728148996e-21, 6.234726653058522e-24, 9.991549123491221e-27, 1.3724655389411017e-29, 1.6338875463584545e-32, 1.7019661941233902e-35, 1.564307163716351e-38))
+            r = imag(r) + real(r) * im  # the following taylor series compute J1 then use relation formula to I1
+        elseif real(_z) < 5.5
+            r = evalpoly(_z - 4.0, (-0.06604332802354913, -0.3806389778579601, 0.07853768224328367, 0.04930901975080153, -0.006977093733246101, -0.0020301389249499455, 0.00022646889994150718, 4.155579989376072e-5, -3.859397238281367e-6, -5.110006887219984e-7, 4.077048904895698e-8, 4.19879116624702e-9, -2.939803646709019e-10, -2.4697315577955285e-11, 1.5414044832012797e-12, 1.0915304728084314e-13, -6.146812221248163e-15, -3.757757785158431e-16, 1.9280657169072558e-17, 1.0361898925247856e-18, -4.88306428816425e-20, -2.340070332638712e-21, 1.0196735825466999e-22, 4.407465239122521e-24, -1.7860193798634053e-25, -7.02892492572196e-27, 2.661920498379524e-28, 9.613567030783265e-30, -3.417184796719832e-31, -1.1400824284559156e-32, 3.818052689081327e-34))
+        elseif real(_z) < 8.5
+            r = evalpoly(_z - 7.0, (-0.004682823482345833, 0.30074824530274785, -0.019188389693537092, -0.0471605174284124, 0.0029657448695426985, 0.002117094601424397, -0.00012058288259054909, -4.528596933543863e-5, 2.3218131613689033e-6, 5.707147746119497e-7, -2.6460803717919574e-8, -4.758581058225322e-9, 2.0083557753034083e-10, 2.8247297290735913e-11, -1.0924079041461721e-12, -1.2558495868983204e-13, 4.477550018840212e-15, 4.340526278677539e-16, -1.43459083340185e-17, -1.200084819927342e-18, 3.6949972617001396e-20, 2.715156932456774e-21, -7.821992893949606e-23, -5.12035844336118e-24, 1.3856216528273523e-25, 8.172858416431775e-27, -2.0848053095060175e-28, -1.1184641449276482e-29, 2.6979469530445177e-31, 1.326904110743154e-32, -3.035362398996416e-34))
+        elseif real(_z) < 11.5
+            r = evalpoly(_z - 10.0, (0.04347274616886144, -0.25028303906823446, -0.009004857400174687, 0.04116523404576117, -0.00023758063287529515, -0.0019744712563326975, 3.38446715384414e-5, 4.430555504324474e-5, -9.214428549130458e-7, -5.769323465195413e-7, 1.2549544791376487e-8, 4.920635393174389e-9, -1.0635889291006728e-10, -2.9682550729139e-11, 6.236949556035066e-13, 1.335190061073787e-13, -2.6994908740535234e-15, -4.655158686182484e-16, 9.012686453800453e-18, 1.2956445537610907e-18, -2.3972087802093897e-20, -2.9464778650364287e-21, 5.2070751129572234e-23, 5.579168848116706e-24, -9.420290550228352e-26, -8.934107001775822e-27, 1.442406752905975e-28, 1.2258566621956357e-29, -1.8943870598688617e-31, -1.4574428507610173e-32, 2.158353205458848e-34))
+        elseif real(_z) < 14.5
+            r = evalpoly(_z - 13.0, (-0.07031805212177837, 0.2123351833095123, 0.02678424666696777, -0.035646496907460384, -0.0015316667868290057, 0.0017627902563241335, 2.9707901726611656e-5, -4.0890119158491976e-5, -2.10088848830435e-7, 5.47805618043409e-7, -4.871508853057117e-10, -4.779338729317937e-9, 2.017453026293242e-11, 2.934694974829395e-11, -1.799648312960788e-13, -1.3386232577358232e-13, 9.690194063808074e-16, 4.718889828158393e-16, -3.7120562820930424e-18, -1.3250234089209876e-18, 1.086246753855577e-20, 3.034889306802292e-21, -2.531740898715398e-23, -5.7802797109249766e-24, 4.83558941632857e-26, 9.301059697909257e-27, -7.730067797438434e-29, -1.281381000050621e-29, 1.0514411994670651e-31, 1.5286677781106927e-32, -1.2332777691576932e-34))
+        else
+            r = evalpoly(_z - 16.0, (0.09039717566130419, -0.1805488974624607, -0.03937987780123671, 0.030669450780043572, 0.0027753540711091436, -0.0015428493442746806, -7.595207975950073e-5, 3.65255237692901e-5, 1.0823603200062054e-6, -4.995971757648329e-7, -9.336967737302558e-9, 4.442044707702236e-9, 5.331367221660131e-11, -2.772590753320406e-11, -2.126094383355898e-13, 1.2822850555243102e-13, 6.076143699602834e-16, -4.57294070904595e-16, -1.2321551004935482e-18, 1.2965544201607567e-18, 1.599193916574117e-21, -2.993977202171527e-21, -5.0754551281589305e-25, 5.7417375563722476e-24, -3.5384202802220665e-27, -9.2932433359164e-27, 1.1177526116165572e-29, 1.2867274495418905e-29, -2.114948449753267e-32, -1.5416757934083927e-32, 3.0470627981401083e-35))
+        end
+        r = conj(r) * im 
     elseif abs(z) < 5.2
+        # power series for I1
         r = z*evalpoly(z*z, (0.5, 0.0625, 0.0026041666666666665, 5.425347222222222e-5, 6.781684027777778e-7, 5.651403356481481e-9, 3.363930569334215e-11, 1.5017547184527747e-13, 5.214426105738801e-16, 1.4484516960385557e-18, 3.2919356728148996e-21, 6.234726653058522e-24, 9.991549123491221e-27, 1.3724655389411017e-29, 1.6338875463584545e-32, 1.7019661941233902e-35, 1.564307163716351e-38))
     else
         if angle(z) <= π / 4.4
+            # power series for I1 with second order Horner scheme
             zz = z*z
             z4 = zz*zz
             r = evalpoly(z4, (0.5, 0.0026041666666666665, 6.781684027777778e-7, 3.363930569334215e-11, 5.214426105738801e-16, 3.2919356728148996e-21, 9.991549123491221e-27, 1.6338875463584545e-32, 1.564307163716351e-38, 9.342315267006072e-45, 3.658488121477942e-51, 9.781133223498595e-58, 1.845775442236299e-64, 2.5281824488224564e-71, 2.574012221626064e-78, 1.9883297967078112e-85, 1.1862954038963049e-92, 5.553068705606664e-100))
             r += zz*evalpoly(z4, (0.0625, 5.425347222222222e-5, 5.651403356481481e-9, 1.5017547184527747e-13, 1.4484516960385557e-18, 6.234726653058522e-24, 1.3724655389411017e-29, 1.7019661941233902e-35, 1.2780287285264307e-41, 6.146260044082942e-48, 1.9797013644361157e-54, 4.4298610613671176e-61, 7.099136316293458e-68, 8.360391695841456e-75, 7.396586843753058e-82, 5.010911786057992e-89, 2.6432607038687723e-96))
             r *= z    
         else
-            zz = (-1.3144285719254842 + 20.054304662400146im, -1.2128746095534062 + 18.50489060934118im,
-                  -1.0886579984944276 + 16.60971135387317im, -0.9427507983505322 + 14.383597659587654im, 
-                  -0.7876981565498402 + 12.017951489232377im, -0.6255235989509672 + 9.543645881424778im, 
-                  -0.4528832294653219 + 6.9096628407010305im, -0.321046242530308 + 4.898219116608337im, 
-                  16.905702128729647 + 10.872315095446137im, 4.119703326027316 + 2.6529313043347877im, 
-                  -1.2665282079632723 + 19.323486333541947im, -0.3857883819775798 + 5.8859932845642575im, 
-                  5.040563157078601 + 19.45771628582095im, 8.609436758553485 + 5.532949040120932im, 
-                  0.694001756513151 + 20.08801537140881im, -0.5457083486368534 + 8.325900481870493im, 
-                  2.790091405909159 + 4.028075216114001im, 12.479613762034912 + 15.75656181882421im, 
-                  -1.0195312436724646 + 15.555040882512422im, 13.06890672928664 + 8.39887636916513im, 
-                  5.1276393456482685 + 3.295333712441098im, 2.7253844460913528 + 19.914373693917753im, 
-                  14.171480465360736 + 9.107457483790645im, 0.8038847904678548 + 4.833608304740307im
-                 )
-            w  = (0.013404709995345458 + 0.0im, -0.10348199668238728 + 0.03550804014042965im, 
-                  -0.11679300367073918 + 0.16207579437860672im, -0.03568356238246015 - 0.16027459313035955im, 
-                  -0.060434459186566966 - 0.02767438331782018im, -0.07086407749919574 + 0.06788442528741509im,
-                  -0.08650269086297906 - 0.029416238653716714im, -0.00415412382674645 - 0.00607269307464697im, 
-                  -0.07587326685881253 + 0.19888070181762688im, -0.08774917184927954 - 0.014101104643344704im, 
-                  0.007193185558007104 + 0.06995793660348812im, -0.02672800181750328 - 0.03807797543563465im, 
-                  0.19776072630033434 + 0.05342203578372187im, 0.3489984378687633 - 0.29126891724262455im, 
-                  -0.010279561959668174 + 0.09872927973051646im, -0.12150771925043644 - 0.012812175864887374im,
-                  -0.1299759980313551 - 0.07418447933244608im, -0.05894210892668888 + 0.24161474004043546im,
-                  -0.2648591701697078 - 0.07362611929106554im, -0.000770797654486471 + 0.00014803398133387874im,
-                  0.12194676645309753 - 0.25480936217584116im, 0.0718337239451382 + 0.20114630545134515im,
-                  0.5172999110471845 - 0.09099393187858433im, -0.023839427516594167 - 0.056052737377955054im
-                  )
-            f = (-3.215770403317309 + 4.184450606321457im, 3.196349209039463 - 3.528571814265253im,
-                 -3.0032703696390066 + 0.8968536198145782im, 1.7024961674754864 + 2.2913462636360844im, 
-                 2.0860054490489937 - 0.919904256366073im, 0.019539554711925467 - 1.3245943976605632im, 
-                 -0.5518276420937166 - 0.23608328648174784im, 0.7317073656843707 + 0.7151746200203796im, 
-                 0.3926351501115586 + 0.00413652715802052im, 0.372728884764365 + 0.01864352217896087im, 
-                 -3.733141424945208 - 2.8413601177219308im, 0.9786880011963323 - 0.618277150144583im, 
-                 0.3971643533953302 + 0.007251821636383407im, 0.38646933950897416 + 0.008351088369392122im,
-                 0.33879033458506147 + 0.08702715214401817im, 1.3930956061616984 + 0.6707847701397638im, 
-                 0.38110598793307165 + 0.02730656109552859im, 0.39435250888241274 + 0.005949183211885836im, 
-                 1.251416864948188 - 2.932439718874822im, 0.39076282374400617 + 0.005394185668824364im, 
-                 0.3778378309823485 + 0.014609112728729376im, 0.39660674100748977 + 0.008337843914085801im, 
-                 0.39140459884033496 + 0.004960218953458415im, 0.4210441290407154 + 0.10788322118417641im
-                ).*w
-
+            # rational approximation using the AAA algorithm
+            zz = (-1.3144285719254842 + 20.054304662400146im, -1.2128746095534062 + 18.50489060934118im, -1.0886579984944276 + 16.60971135387317im, -0.9427507983505322 + 14.383597659587654im, -0.7876981565498402 + 12.017951489232377im, -0.6255235989509672 + 9.543645881424778im, -0.4528832294653219 + 6.9096628407010305im, -0.321046242530308 + 4.898219116608337im, 16.905702128729647 + 10.872315095446137im, 4.119703326027316 + 2.6529313043347877im, -1.2665282079632723 + 19.323486333541947im, -0.3857883819775798 + 5.8859932845642575im, 5.040563157078601 + 19.45771628582095im, 8.609436758553485 + 5.532949040120932im, 0.694001756513151 + 20.08801537140881im, -0.5457083486368534 + 8.325900481870493im, 2.790091405909159 + 4.028075216114001im, 12.479613762034912 + 15.75656181882421im, -1.0195312436724646 + 15.555040882512422im, 13.06890672928664 + 8.39887636916513im, 5.1276393456482685 + 3.295333712441098im, 2.7253844460913528 + 19.914373693917753im, 14.171480465360736 + 9.107457483790645im, 0.8038847904678548 + 4.833608304740307im)
+            w  = (0.013404709995345458 + 0.0im, -0.10348199668238728 + 0.03550804014042965im, -0.11679300367073918 + 0.16207579437860672im, -0.03568356238246015 - 0.16027459313035955im, -0.060434459186566966 - 0.02767438331782018im, -0.07086407749919574 + 0.06788442528741509im, -0.08650269086297906 - 0.029416238653716714im, -0.00415412382674645 - 0.00607269307464697im, -0.07587326685881253 + 0.19888070181762688im, -0.08774917184927954 - 0.014101104643344704im, 0.007193185558007104 + 0.06995793660348812im, -0.02672800181750328 - 0.03807797543563465im, 0.19776072630033434 + 0.05342203578372187im, 0.3489984378687633 - 0.29126891724262455im, -0.010279561959668174 + 0.09872927973051646im, -0.12150771925043644 - 0.012812175864887374im, -0.1299759980313551 - 0.07418447933244608im, -0.05894210892668888 + 0.24161474004043546im, -0.2648591701697078 - 0.07362611929106554im, -0.000770797654486471 + 0.00014803398133387874im, 0.12194676645309753 - 0.25480936217584116im, 0.0718337239451382 + 0.20114630545134515im, 0.5172999110471845 - 0.09099393187858433im, -0.023839427516594167 - 0.056052737377955054im)
+            f = (-3.215770403317309 + 4.184450606321457im, 3.196349209039463 - 3.528571814265253im, -3.0032703696390066 + 0.8968536198145782im, 1.7024961674754864 + 2.2913462636360844im, 2.0860054490489937 - 0.919904256366073im, 0.019539554711925467 - 1.3245943976605632im, -0.5518276420937166 - 0.23608328648174784im, 0.7317073656843707 + 0.7151746200203796im, 0.3926351501115586 + 0.00413652715802052im, 0.372728884764365 + 0.01864352217896087im, -3.733141424945208 - 2.8413601177219308im, 0.9786880011963323 - 0.618277150144583im, 0.3971643533953302 + 0.007251821636383407im, 0.38646933950897416 + 0.008351088369392122im, 0.33879033458506147 + 0.08702715214401817im, 1.3930956061616984 + 0.6707847701397638im, 0.38110598793307165 + 0.02730656109552859im, 0.39435250888241274 + 0.005949183211885836im, 1.251416864948188 - 2.932439718874822im, 0.39076282374400617 + 0.005394185668824364im, 0.3778378309823485 + 0.014609112728729376im, 0.39660674100748977 + 0.008337843914085801im, 0.39140459884033496 + 0.004960218953458415im, 0.4210441290407154 + 0.10788322118417641im)
+            f = f.*w
             s1 = 0.0 + 0.0im
             s2 = 0.0 + 0.0im
             @fastmath for ind in eachindex(f)
@@ -330,7 +316,6 @@ function besseli1(z::ComplexF32)
     else
         isconj = false
     end
-
     if abs(z) > 10.0
         zinv = 1 / z
         e = exp(z)
