Skip to content

Commit 15f9f3e

Browse files
heltonmcheltonmc
committed
fix ordering
1 parent ba4405d commit 15f9f3e

File tree

1 file changed

+60
-112
lines changed

1 file changed

+60
-112
lines changed

‎src/i0.jl

Lines changed: 60 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,4 @@
1-
"""
2-
besseli0(x::T) where T <: Union{Float32, Float64}
3-
4-
Modified Bessel function of the first kind of order zero, ``I_0(x)``.
5-
"""
6-
function besseli0(z::Complex{T}) where T <: Union{Float64}
1+
function besseli0(z::ComplexF64)
72
check = false
83
if real(z) < 0.0
94
z = -z
@@ -13,115 +8,71 @@ function besseli0(z::Complex{T}) where T <: Union{Float64}
138
check = true
149
end
1510

16-
if abs(z) > 20.0
11+
if abs(z) > 17.5
1712
zinv = 1 / z
1813
e = exp(z)
1914
s = sqrt(2 * z * π)
20-
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))
21-
p2 = zinv*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))
15+
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))
16+
p2 = zinv*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))
2217
r = e / s * (p + p2) + im * (p - p2) / (e * s)
23-
elseif abs(z) < 5.0
18+
elseif abs(z) < 5.2
2419
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))
2520
else
26-
if angle(z) <= pi/5
27-
#r = besseli_power_series(0.0, z)
28-
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, 4.60090342269515e-40, 3.5500798014623073e-43, 2.458504017633177e-46, 1.5365650110207356e-49, 8.71068600351891e-53, 4.499321282809354e-56, 2.1263333094562166e-59, 9.228877211181495e-63, 3.691550884472598e-66, 1.3652185223641266e-69, 4.681819349671216e-73, 1.4929270885431173e-76, 4.4379521062518346e-80, 1.232764473958843e-83, 3.206983543077115e-87, 7.829549665715613e-91, 1.7974172786307652e-94, 3.887148093924665e-98, 7.932955293723807e-102))
21+
if angle(z) <= π / 4.4
22+
zz = z*z
23+
zz2 = zz*zz
24+
r = evalpoly(zz2, (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))
25+
r += zz*evalpoly(zz2, (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))
2926
else
30-
zz = (
31-
-1.31445408166514621228770920424722135067 + 20.05469386581672353031535749323666095734im,
32-
-1.21527202201844963802557231247192248702 + 18.54146805523843966057029319927096366882im,
33-
-1.075858296896808452558502722240518778563 + 16.41442564500392009563256578985601663589im,
34-
-0.9186088936209556576883983325387816876173 + 14.01526337127532251258799078641459345818im,
35-
-0.7533475204960968785172781281289644539356 + 11.49386205943561911624328786274418234825im,
36-
-0.5838248556039682402030166485928930342197 + 8.907445998843897427832416724413633346558im,
37-
-0.4109216146573061445579355677182320505381 + 6.269452314652045998855101061053574085236im,
38-
-0.3205522729306593543441294968943111598492 + 4.890682596894066591630689799785614013672im,
39-
16.90671675080808711300051072612404823303 + 10.86528710778251571866803715238347649574im,
40-
0.2671537175414016584973353474197210744023 + 20.09822452086760335987491998821496963501im,
41-
4.121994532242283959533324377844110131264 + 2.649369939469517376551266352180391550064im,
42-
4.053028754723364990297795884544029831886 + 19.68712670537236064660646661650389432907im,
43-
-0.5018982035098999983091516696731559932232 + 7.657486833198152709201167454011738300323im,
44-
9.416081350634247115749531076289713382721 + 17.75802387649701330474272253923118114471im,
45-
1.967489559014300670725106101599521934986 + 4.487648029332259369539315230213105678558im,
46-
-1.275759132467892076334692319505847990513 + 19.46432302583940909812554309610277414322im,
47-
7.709313889965882182764289609622210264206 + 4.954475197822860721430515695828944444656im,
48-
13.72022229845232033085267175920307636261 + 14.68895844098729064342023775679990649223im,
49-
-0.3541218766747901147695642976032104343176 + 5.402855776372860852063695347169414162636im,
50-
-0.9973883749616251348513173979881685227156 + 15.21720599006466656533120840322226285934im,
51-
5.067285183070907805813476443290710449219 + 3.256546447342954841985829261830076575279im,
52-
2.130947270892558531585336822899989783764 + 19.98672218570808212234624079428613185883im,
53-
-0.8162780120038937162418819612565136081901 + 12.45399582113806502547959098592400550842im,
54-
-0.6642234573493526195164804448722861707211 + 10.13409162133750207601678994251415133476im
55-
)
56-
57-
58-
w = ( 0.02343591247913679245784557281240267911926 + 0.0im,
59-
- 0.07518141047674767318831356988084735348821 + 0.04759481977311690037435454314618255011737im,
60-
0.02099862365431344121691203952195792226121 + 0.1022779931770602390717073149062343873084im,
61-
- 0.04479060613413702457430431991269870195538 - 0.001990058142152528741775086018606089055538im,
62-
0.05092855107420864863021492396910616662353 - 0.03028402769108835476674634890059678582475im,
63-
0.1187682211635005813388232809302280656993 + 0.3545737241031228781373840774904238060117im,
64-
- 0.2800428715703842108553089929046109318733 + 0.1385502878368521095797660791504313237965im,
65-
0.07056156978707124605154632490666699595749 - 0.01564004921737297687522882938537804875523im,
66-
- 0.06946091028316167537148828614590456709266 + 0.01606645463234542339781008024601760553196im,
67-
0.08524646892263365582920187080162577331066 + 0.09536227567093939760933807292531128041446im,
68-
- 0.02790219158045261979572693178397457813844 + 0.0706492680230463854229583375854417681694im,
69-
0.142231423342524893049798606625699903816 - 0.2536760604980480837689071904605953022838im,
70-
- 0.2058636344291391306882132994360290467739 + 0.301438312655757711944204402243485674262im,
71-
0.2054714553805426502375297559410682879388 - 0.1738865814674955823093682738544885069132im,
72-
- 0.1383969021060848236803764166324981488287 - 0.07997812132261208906136573659750865772367im,
73-
- 0.005116954351018287300290054986362520139664 + 0.07859682420061694929636075812595663592219im,
74-
0.04094909855072709908840877801594615448266 - 0.264288102459465323867959796189097687602im,
75-
- 0.04174610452138102778540940107632195577025 - 0.2381078364144596504203832409984897822142im,
76-
- 0.1407821015956153554160579233212047256529 - 0.1764715349923221265893147347014746628702im,
77-
- 0.07672294420555977878528608471242478117347 + 0.08789101702567167495594446791074005886912im,
78-
- 0.1554449133166528329574873623641906306148 - 0.1787213424643953330051004968481720425189im,
79-
0.2746212867658748835175686053844401612878 + 0.06339130601389010577495497500422061420977im,
80-
0.003561860807958012152540927530708358972333 + 0.005781000943183686215098848748539239750244im,
81-
0.2246761482071190363374313392341719008982 + 0.05087374143304821544342431138829851988703im)
82-
83-
84-
f = ( 4.117575040058659929798068333184346556664 - 4.095512734321109071800037781940773129463im,
85-
- 2.247112620034342089780921014607883989811 + 3.680759784557968927742876985576003789902im,
86-
3.78343182981733150427317013964056968689 + 0.5622435659838359578444055841828230768442im,
87-
1.02504489923751607172164312942186370492 - 2.429576351715230231320674647577106952667im,
88-
- 1.102472896641811361817531178530771285295 - 0.9985456859791927985980919402209110558033im,
89-
- 0.713405219160451298243685869238106533885 + 0.6341686986762049560439891138230450451374im,
90-
0.3548943281905181379443092737346887588501 + 0.8984965319188545906925469353154767304659im,
91-
0.1510514019542046337818419488030485808849 - 0.7254114647779558167073332697327714413404im,
92-
0.4010578521472479840426217378990259021521 - 0.00140846956682444488656580361407577584032im,
93-
0.540457156458840848323177397105609998107 - 0.1885021348797114026929477859084727242589im,
94-
0.4077709777425669868122781736019533127546 - 0.006802901797014750291670015514000624534674im,
95-
0.3994965367795383914817364257032750174403 - 0.002466114177937662587519751511422327894252im,
96-
0.8312271064319111113505300636461470276117 - 1.005050393663784458198051652288995683193im,
97-
0.400061910584651958533441984400269575417 - 0.002248719754950747901772745152015886560548im,
98-
0.4055260568863696124530804354435531422496 - 0.01689191195791546126758753132435231236741im,
99-
5.210832430161148387526282022008672356606 + 1.741327835783517130607833678368479013443im,
100-
0.4036474747676551677599832146370317786932 - 0.00328838026637900005672010550483719271142im,
101-
0.4006279890770570450975185394781874492764 - 0.001884565453597092116871936084976368874777im,
102-
- 0.393832011865159148378268127999035641551 - 0.1798762838433343724808821662008995190263im,
103-
- 2.053184226968296055559903834364376962185 + 1.606906650995165719564283790532499551773im,
104-
0.406174279778613755986782507534371688962 - 0.005320233415147143825330022792741146986373im,
105-
0.4034434765080088802768898403883213177323 - 0.006086635360116373315297888524355585104786im,
106-
- 0.07645988937623526826570241610170342028141 + 1.981923553555364980738318081421311944723im,
107-
1.884413412262007314623701859090942889452 + 0.2417979244204202515788892924319952726364im
108-
)
109-
d = w.*f
110-
111-
s1 = 0.0 + 0.0im
112-
s2 = 0.0 + 0.0im
113-
114-
@fastmath for ind in eachindex(f)
115-
C = inv(z - zz[ind])
116-
s1 += C*d[ind]
117-
s2 += C*w[ind]
118-
end
119-
120-
121-
#C = 1 ./ (z .- transpose(zz))
122-
#r = (C*d) ./ (C*w)
123-
124-
r = (s1/s2) / (sqrt(z) * exp(-z))
27+
zz = (-1.3144540816651462 + 20.054693865816724im, -1.2152720220184496 + 18.54146805523844im,
28+
-1.0758582968968085 + 16.41442564500392im, -0.9186088936209557 + 14.015263371275323im,
29+
-0.7533475204960969 + 11.49386205943562im, -0.5838248556039682 + 8.907445998843897im,
30+
-0.41092161465730614 + 6.269452314652046im, -0.32055227293065935 + 4.890682596894067im,
31+
16.906716750808087 + 10.865287107782516im, 0.26715371754140166 + 20.098224520867603im,
32+
4.121994532242284 + 2.6493699394695174im, 4.053028754723365 + 19.68712670537236im,
33+
-0.5018982035099 + 7.657486833198153im, 9.416081350634247 + 17.758023876497013im,
34+
1.9674895590143007 + 4.487648029332259im, -1.275759132467892 + 19.46432302583941im,
35+
7.709313889965882 + 4.954475197822861im, 13.72022229845232 + 14.68895844098729im,
36+
-0.3541218766747901 + 5.402855776372861im, -0.9973883749616251 + 15.217205990064667im,
37+
5.067285183070908 + 3.256546447342955im, 2.1309472708925585 + 19.986722185708082im,
38+
-0.8162780120038937 + 12.453995821138065im, -0.6642234573493526 + 10.134091621337502im
39+
)
40+
w = (0.023435912479136792 + 0.0im, -0.07518141047674767 + 0.0475948197731169im,
41+
0.02099862365431344 + 0.10227799317706024im, -0.044790606134137025 - 0.0019900581421525287im,
42+
0.05092855107420865 - 0.030284027691088355im, 0.11876822116350058 + 0.3545737241031229im,
43+
-0.2800428715703842 + 0.1385502878368521im, 0.07056156978707125 - 0.015640049217372977im,
44+
-0.06946091028316168 + 0.016066454632345423im, 0.08524646892263366 + 0.0953622756709394im,
45+
-0.02790219158045262 + 0.07064926802304639im, 0.1422314233425249 - 0.2536760604980481im,
46+
-0.20586363442913913 + 0.3014383126557577im, 0.20547145538054265 - 0.17388658146749558im,
47+
-0.13839690210608482 - 0.07997812132261209im, -0.005116954351018287 + 0.07859682420061695im,
48+
0.0409490985507271 - 0.2642881024594653im, -0.04174610452138103 - 0.23810783641445965im,
49+
-0.14078210159561536 - 0.17647153499232213im, -0.07672294420555978 + 0.08789101702567167im,
50+
-0.15544491331665283 - 0.17872134246439533im, 0.2746212867658749 + 0.0633913060138901im,
51+
0.003561860807958012 + 0.005781000943183686im, 0.22467614820711904 + 0.050873741433048215im
52+
)
53+
f = (4.11757504005866 - 4.095512734321109im, -2.247112620034342 + 3.680759784557969im,
54+
3.7834318298173315 + 0.562243565983836im, 1.025044899237516 - 2.4295763517152302im,
55+
-1.1024728966418114 - 0.9985456859791928im, -0.7134052191604513 + 0.634168698676205im,
56+
0.35489432819051814 + 0.8984965319188546im, 0.15105140195420463 - 0.7254114647779558im,
57+
0.401057852147248 - 0.0014084695668244449im, 0.5404571564588408 - 0.1885021348797114im,
58+
0.407770977742567 - 0.00680290179701475im, 0.3994965367795384 - 0.0024661141779376626im,
59+
0.8312271064319111 - 1.0050503936637845im, 0.40006191058465196 - 0.002248719754950748im,
60+
0.4055260568863696 - 0.01689191195791546im, 5.210832430161148 + 1.7413278357835171im,
61+
0.40364747476765517 - 0.003288380266379im, 0.40062798907705705 - 0.0018845654535970921im,
62+
-0.39383201186515915 - 0.17987628384333437im, -2.053184226968296 + 1.6069066509951657im,
63+
0.40617427977861376 - 0.005320233415147144im, 0.4034434765080089 - 0.006086635360116373im,
64+
-0.07645988937623527 + 1.981923553555365im, 1.8844134122620073 + 0.24179792442042025im
65+
).*w
66+
67+
s1 = 0.0 + 0.0im
68+
s2 = 0.0 + 0.0im
69+
@fastmath for ind in eachindex(f)
70+
C = inv(z - zz[ind])
71+
s1 += C*f[ind]
72+
s2 += C*w[ind]
73+
end
74+
75+
r = s1 / (s2 * sqrt(z) * exp(-z))
12576
end
12677
end
12778
check && (r = conj(r))
@@ -152,7 +103,4 @@ function k0_large_argument(z::T) where T
152103
p = evalpoly(1/z, (1.0, -1/8, 9/128, -75/1024, 3675/32768, -59535/262144, 2401245/4194304, -57972915/33554432, 13043905875/2147483648, -418854310875/17179869184, 30241281245175/274877906944, -1212400457192925/2199023255552, 213786613951685775/70368744177664, -10278202593831046875/562949953421312, 1070401384414690453125/9007199254740992, -60013837619516978071875/72057594037927936, 57673297952355815927071875/9223372036854775808, -3694483615889146090857721875/73786976294838206464, 502860269940467106811189921875/1180591620717411303424, -36232405765710498380237842265625/9444732965739290427392, 11021897833929133607268351617203125/302231454903657293676544, -882276678992136837800861860405640625/2417851639229258349412352, T(148302689041496455735799416353639046875)/T(38685626227668133590597632)))
153104
return c*p
154105
end
155-
=#
156-
157-
158-
## these are pretty good but still need some work
106+
=#

0 commit comments

Comments
 (0)