use taylor expansions around roots

heltonmc · heltonmc · commit 741bf234c835 · 2022-12-01T14:13:19.000-05:00
diff --git a/src/besseli.jl b/src/besseli.jl
@@ -143,72 +143,59 @@ function besseli0(z::ComplexF64)
         z = conj(z)
         isconj = true
     end
-
     if abs(z) > 17.5
+        # use asymptotic expansion for large arguments
         zinv = 1 / z
         e = exp(z)
         sinv =sqrt(zinv)/sqrt(2*pi)
         p = evalpoly(zinv*zinv, (1.0, 0.0703125, 0.112152099609375, 0.5725014209747314, 6.074042001273483, 110.01714026924674, 3038.090510922384, 118838.42625678325, 6.252951493434797e6, 4.259392165047669e8, 3.646840080706556e10, 3.8335346613939443e12, 4.8540146868529006e14, 7.286857349377656e16))
         p2 = evalpoly(zinv*zinv, (0.125, 0.0732421875, 0.22710800170898438, 1.7277275025844574, 24.380529699556064, 551.3358961220206, 18257.755474293175, 832859.3040162893, 5.0069589531988926e7, 3.8362551802304335e9, 3.6490108188498334e11, 4.218971570284097e13, 5.827244631566907e15, 9.47628809926011e17))
         r = e * sinv * muladd(p2, zinv, p) + im * muladd(p2, -zinv, p) * sinv / e
-    elseif abs(z) < 5.2
+    elseif real(z) < 4.5
+        # use taylor series around the roots (slight offset) of J0(z)
+        # use relation I0(z) = J0(im*z)
+        _z = imag(z) + abs(real(z))*im
+        if real(_z) < 2.2
+            # use power series for I0
+            r = evalpoly(z*z, (1.0, 0.25, 0.015625, 0.00043402777777777775, 6.781684027777777e-6, 6.781684027777778e-8, 4.709502797067901e-10, 2.4028075495244395e-12, 9.385966990329842e-15, 2.896903392077112e-17, 7.242258480192779e-20, 1.4963343967340453e-22, 2.5978027721077174e-25, 3.842903509035085e-28, 4.9016626390753635e-31, 5.4462918211948485e-34))
+            r = conj(r) # the following taylor series compute J0 then use relation formula to I0
+        elseif real(_z) < 5.5
+            r = evalpoly(_z - 4.0, (-0.39714980986384735, 0.06604332802354913, 0.19031948892898004, -0.02617922741442789, -0.012327254937700382, 0.0013954187466492201, 0.0003383564874916576, -3.2352699991643884e-5, -5.19447498672009e-6, 4.288219153645963e-7, 5.110006887219985e-8, -3.706408095359726e-9, -3.498992638539183e-10, 2.261387420545399e-11, 1.764093969853949e-12, -1.0276029888008532e-13, -6.822065455052696e-15, 3.615771894851861e-16, 2.087643213976906e-17, -1.0147714299511874e-18, -5.180949462623928e-20, 2.325268708649643e-21, 1.0636683330175965e-22, -4.433363402376956e-24, -1.8364438496343838e-25, 7.144077519453622e-27, 2.703432663739215e-28, -9.858964808813052e-30, -3.433416796708309e-31, 1.1783395850758042e-32, 3.8002747615197184e-34))
+        elseif real(_z) < 8.5
+            r = evalpoly(_z - 7.0, (0.3000792705195556, 0.004682823482345833, -0.15037412265137393, 0.0063961298978456975, 0.0117901293571031, -0.0005931489739085397, -0.00035284910023739957, 1.7226126084364155e-5, 5.6607461669298285e-6, -2.5797924015210036e-7, -5.7071477461194966e-8, 2.405527610719961e-9, 3.9654842151877684e-10, -1.5448890579256987e-11, -2.0176640921954225e-12, 7.282719360974481e-14, 7.849059918114502e-15, -2.6338529522589484e-16, -2.4114034881541882e-17, 7.550478070536054e-19, 6.00042409963671e-20, -1.7595225055714948e-21, -1.2341622420258064e-22, 3.400866475630264e-24, 2.1334826847338247e-25, -5.542486611309409e-27, -3.14340708324299e-28, 7.721501146318583e-30, 3.994514803313029e-31, -9.303265355325922e-33, -4.42301370247718e-34))
+        elseif real(_z) < 11.5
+            r = evalpoly(_z - 10.0, (-0.24593576445134835, -0.04347274616886144, 0.12514151953411723, 0.003001619133391562, -0.010291308511440292, 4.751612657505903e-5, 0.0003290785427221163, -4.8349530769202e-6, -5.5381943804055925e-6, 1.0238253943478287e-7, 5.769323465195413e-8, -1.1408677083069533e-9, -4.100529494311991e-10, 8.181453300774406e-12, 2.1201821949384996e-12, -4.157966370690044e-14, -8.344937881711169e-15, 1.5879358082667785e-16, 2.58619927010138e-17, -4.743519186210765e-19, -6.478222768805454e-20, 1.1415279905758999e-21, 1.339308120471104e-22, -2.2639457012857495e-24, -2.3246536867152943e-25, 3.768116220091341e-27, 3.4361950006830084e-28, -5.342247232985093e-30, -4.378059507841556e-31, 6.532369171961593e-33, 4.8581428358700575e-34))
+        elseif real(_z) < 14.5
+            r = evalpoly(_z - 13.0, (0.20692610237706782, 0.07031805212177837, -0.10616759165475616, -0.00892808222232259, 0.008911624226865098, 0.00030633335736580117, -0.00029379837605402224, -4.243985960944521e-6, 5.111264894811495e-6, 2.334320542560391e-8, -5.478056180434088e-8, 4.428644411870052e-11, 3.982782274431613e-10, -1.5518869433024987e-12, -2.0962106963067093e-12, 1.1997655419738472e-14, 8.3663953608489e-15, -5.700114155181254e-17, -2.6216054600879916e-17, 1.9537138326805023e-19, 6.625117044604972e-20, -5.172603589788702e-22, -1.3794951394555679e-22, 1.100756912484797e-24, 2.408449879552199e-25, -1.93423576653241e-27, -3.5773306530412485e-28, 2.8629880731192636e-30, 4.576360714471254e-31, -3.62565930854413e-33, -5.095559260340188e-34, 3.978315384156111e-36))
+        else
+            r = evalpoly(_z - 16.0, (-0.1748990739836292, -0.09039717566130419, 0.09027444873123035, 0.01312662593374557, -0.007667362695010893, -0.0005550708142218287, 0.0002571415573791134, 1.0850297108500103e-5, -4.565690471161262e-6, -1.2026225777846727e-7, 4.9959717576483296e-8, 8.48815248845687e-10, -3.701703923085197e-10, -4.101051708969331e-12, 1.980421966657433e-12, 1.4173962555705988e-14, -8.014281597026939e-15, -3.5742021762369616e-17, 2.540522616136639e-17, 6.485026844702886e-20, -6.482772100803784e-20, -7.615209126543415e-23, 1.3608987282597849e-22, 2.2067196209386653e-26, -2.3923906484884365e-25, 1.4153681120888267e-28, 3.5743243599678463e-28, -4.139824487468731e-31, -4.595455176935322e-31, 7.2929256888043685e-34, 5.138919311361309e-34))
+        end
+        r = conj(r) 
+    elseif abs(z) < 5.5
+        # use power series for I0
         r = evalpoly(z*z, (1.0, 0.25, 0.015625, 0.00043402777777777775, 6.781684027777777e-6, 6.781684027777778e-8, 4.709502797067901e-10, 2.4028075495244395e-12, 9.385966990329842e-15, 2.896903392077112e-17, 7.242258480192779e-20, 1.4963343967340453e-22, 2.5978027721077174e-25, 3.842903509035085e-28, 4.9016626390753635e-31, 5.4462918211948485e-34, 5.318644356635594e-37))
     else
         if angle(z) <= π / 4.4
+            # use power series but evaluated using second order horner scheme
             zz = z*z
             z4 = zz*zz
             r = evalpoly(z4, (1.0, 0.015625, 6.781684027777777e-6, 4.709502797067901e-10, 9.385966990329842e-15, 7.242258480192779e-20, 2.5978027721077174e-25, 4.9016626390753635e-31, 5.318644356635594e-37, 3.5500798014623073e-43, 1.5365650110207356e-49, 4.499321282809354e-56, 9.228877211181495e-63, 1.3652185223641266e-69, 1.4929270885431173e-76, 1.232764473958843e-83, 7.829549665715613e-91, 3.887148093924665e-98))
             r += zz*evalpoly(z4, (0.25, 0.00043402777777777775, 6.781684027777778e-8, 2.4028075495244395e-12, 2.896903392077112e-17, 1.4963343967340453e-22, 3.842903509035085e-28, 5.4462918211948485e-34, 4.60090342269515e-40, 2.458504017633177e-46, 8.71068600351891e-53, 2.1263333094562166e-59, 3.691550884472598e-66, 4.681819349671216e-73, 4.4379521062518346e-80, 3.206983543077115e-87, 1.7974172786307652e-94, 7.932955293723807e-102))
+        
         else
-            zz = (-1.3144540816651462 + 20.054693865816724im, -1.2152720220184496 + 18.54146805523844im, 
-                  -1.0758582968968085 + 16.41442564500392im, -0.9186088936209557 + 14.015263371275323im, 
-                  -0.7533475204960969 + 11.49386205943562im, -0.5838248556039682 + 8.907445998843897im, 
-                  -0.41092161465730614 + 6.269452314652046im, -0.32055227293065935 + 4.890682596894067im, 
-                  16.906716750808087 + 10.865287107782516im, 0.26715371754140166 + 20.098224520867603im, 
-                  4.121994532242284 + 2.6493699394695174im, 4.053028754723365 + 19.68712670537236im, 
-                  -0.5018982035099 + 7.657486833198153im, 9.416081350634247 + 17.758023876497013im, 
-                  1.9674895590143007 + 4.487648029332259im, -1.275759132467892 + 19.46432302583941im, 
-                  7.709313889965882 + 4.954475197822861im, 13.72022229845232 + 14.68895844098729im, 
-                  -0.3541218766747901 + 5.402855776372861im, -0.9973883749616251 + 15.217205990064667im, 
-                  5.067285183070908 + 3.256546447342955im, 2.1309472708925585 + 19.986722185708082im, 
-                  -0.8162780120038937 + 12.453995821138065im, -0.6642234573493526 + 10.134091621337502im
-                  )
-            w  = (0.023435912479136792 + 0.0im, -0.07518141047674767 + 0.0475948197731169im, 
-                  0.02099862365431344 + 0.10227799317706024im, -0.044790606134137025 - 0.0019900581421525287im, 
-                  0.05092855107420865 - 0.030284027691088355im, 0.11876822116350058 + 0.3545737241031229im, 
-                  -0.2800428715703842 + 0.1385502878368521im, 0.07056156978707125 - 0.015640049217372977im, 
-                  -0.06946091028316168 + 0.016066454632345423im, 0.08524646892263366 + 0.0953622756709394im, 
-                  -0.02790219158045262 + 0.07064926802304639im, 0.1422314233425249 - 0.2536760604980481im, 
-                  -0.20586363442913913 + 0.3014383126557577im, 0.20547145538054265 - 0.17388658146749558im, 
-                  -0.13839690210608482 - 0.07997812132261209im, -0.005116954351018287 + 0.07859682420061695im, 
-                  0.0409490985507271 - 0.2642881024594653im, -0.04174610452138103 - 0.23810783641445965im, 
-                  -0.14078210159561536 - 0.17647153499232213im, -0.07672294420555978 + 0.08789101702567167im, 
-                  -0.15544491331665283 - 0.17872134246439533im, 0.2746212867658749 + 0.0633913060138901im, 
-                  0.003561860807958012 + 0.005781000943183686im, 0.22467614820711904 + 0.050873741433048215im
-                  )
-            f = (4.11757504005866 - 4.095512734321109im, -2.247112620034342 + 3.680759784557969im, 
-                 3.7834318298173315 + 0.562243565983836im, 1.025044899237516 - 2.4295763517152302im, 
-                 -1.1024728966418114 - 0.9985456859791928im, -0.7134052191604513 + 0.634168698676205im, 
-                 0.35489432819051814 + 0.8984965319188546im, 0.15105140195420463 - 0.7254114647779558im, 
-                 0.401057852147248 - 0.0014084695668244449im, 0.5404571564588408 - 0.1885021348797114im,
-                 0.407770977742567 - 0.00680290179701475im, 0.3994965367795384 - 0.0024661141779376626im, 
-                 0.8312271064319111 - 1.0050503936637845im, 0.40006191058465196 - 0.002248719754950748im, 
-                 0.4055260568863696 - 0.01689191195791546im, 5.210832430161148 + 1.7413278357835171im, 
-                 0.40364747476765517 - 0.003288380266379im, 0.40062798907705705 - 0.0018845654535970921im, 
-                 -0.39383201186515915 - 0.17987628384333437im, -2.053184226968296 + 1.6069066509951657im, 
-                 0.40617427977861376 - 0.005320233415147144im, 0.4034434765080089 - 0.006086635360116373im,
-                 -0.07645988937623527 + 1.981923553555365im, 1.8844134122620073 + 0.24179792442042025im
-                ).*w
-
+            # use rational approximation based on AAA algorithm
+            zz = (-1.3144540816651462 + 20.054693865816724im, -1.2152720220184496 + 18.54146805523844im, -1.0758582968968085 + 16.41442564500392im, -0.9186088936209557 + 14.015263371275323im, -0.7533475204960969 + 11.49386205943562im, -0.5838248556039682 + 8.907445998843897im, -0.41092161465730614 + 6.269452314652046im, -0.32055227293065935 + 4.890682596894067im, 16.906716750808087 + 10.865287107782516im, 0.26715371754140166 + 20.098224520867603im, 4.121994532242284 + 2.6493699394695174im, 4.053028754723365 + 19.68712670537236im, -0.5018982035099 + 7.657486833198153im, 9.416081350634247 + 17.758023876497013im, 1.9674895590143007 + 4.487648029332259im, -1.275759132467892 + 19.46432302583941im, 7.709313889965882 + 4.954475197822861im, 13.72022229845232 + 14.68895844098729im, -0.3541218766747901 + 5.402855776372861im, -0.9973883749616251 + 15.217205990064667im, 5.067285183070908 + 3.256546447342955im, 2.1309472708925585 + 19.986722185708082im, -0.8162780120038937 + 12.453995821138065im, -0.6642234573493526 + 10.134091621337502im)
+            w  = (0.023435912479136792 + 0.0im, -0.07518141047674767 + 0.0475948197731169im, 0.02099862365431344 + 0.10227799317706024im, -0.044790606134137025 - 0.0019900581421525287im, 0.05092855107420865 - 0.030284027691088355im, 0.11876822116350058 + 0.3545737241031229im, -0.2800428715703842 + 0.1385502878368521im, 0.07056156978707125 - 0.015640049217372977im, -0.06946091028316168 + 0.016066454632345423im, 0.08524646892263366 + 0.0953622756709394im, -0.02790219158045262 + 0.07064926802304639im, 0.1422314233425249 - 0.2536760604980481im, -0.20586363442913913 + 0.3014383126557577im, 0.20547145538054265 - 0.17388658146749558im, -0.13839690210608482 - 0.07997812132261209im, -0.005116954351018287 + 0.07859682420061695im, 0.0409490985507271 - 0.2642881024594653im, -0.04174610452138103 - 0.23810783641445965im, -0.14078210159561536 - 0.17647153499232213im, -0.07672294420555978 + 0.08789101702567167im, -0.15544491331665283 - 0.17872134246439533im, 0.2746212867658749 + 0.0633913060138901im, 0.003561860807958012 + 0.005781000943183686im, 0.22467614820711904 + 0.050873741433048215im)
+            f = (4.11757504005866 - 4.095512734321109im, -2.247112620034342 + 3.680759784557969im, 3.7834318298173315 + 0.562243565983836im, 1.025044899237516 - 2.4295763517152302im, -1.1024728966418114 - 0.9985456859791928im, -0.7134052191604513 + 0.634168698676205im, 0.35489432819051814 + 0.8984965319188546im, 0.15105140195420463 - 0.7254114647779558im, 0.401057852147248 - 0.0014084695668244449im, 0.5404571564588408 - 0.1885021348797114im, 0.407770977742567 - 0.00680290179701475im, 0.3994965367795384 - 0.0024661141779376626im, 0.8312271064319111 - 1.0050503936637845im, 0.40006191058465196 - 0.002248719754950748im, 0.4055260568863696 - 0.01689191195791546im, 5.210832430161148 + 1.7413278357835171im, 0.40364747476765517 - 0.003288380266379im, 0.40062798907705705 - 0.0018845654535970921im, -0.39383201186515915 - 0.17987628384333437im, -2.053184226968296 + 1.6069066509951657im, 0.40617427977861376 - 0.005320233415147144im, 0.4034434765080089 - 0.006086635360116373im, -0.07645988937623527 + 1.981923553555365im, 1.8844134122620073 + 0.24179792442042025im)
+            f = f.*w
             s1 = 0.0 + 0.0im
             s2 = 0.0 + 0.0im
             @fastmath for ind in eachindex(f)
                 C = inv(z - zz[ind])
                 s1 += C*f[ind]
                 s2 += C*w[ind]
             end
-
-        r = s1 / (s2 * sqrt(z) * exp(-z))
+            r = s1 / (s2 * sqrt(z) * exp(-z))
         end
     end
     isconj && (r = conj(r))
