Skip to content

Commit ba4405d

Browse files
heltonmcheltonmc
committed
improve coefficients
1 parent bef8cb8 commit ba4405d

File tree

1 file changed

+97
-84
lines changed

1 file changed

+97
-84
lines changed

‎src/i0.jl

Lines changed: 97 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -27,91 +27,101 @@ function besseli0(z::Complex{T}) where T <: Union{Float64}
2727
#r = besseli_power_series(0.0, z)
2828
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))
2929
else
30-
zz = [
31-
-1.678231235746891147897486007423140108585 + 20.29678265316490737291132973041385412216im,
32-
-1.523025526043160660094599734293296933174 + 18.74439925781758020661982300225645303726im,
33-
-1.602232319163288343233375599083956331015 + 19.57994203993233739424795203376561403275im,
34-
-1.297548314533615521071396869956515729427 + 15.67962995660580283185936423251405358315im,
35-
-1.083284426743612138821504231600556522608 + 13.2089741880721351918737127562053501606im,
36-
-0.8756326508559164611966707525425590574741 + 10.59162940924948692611451406264677643776im,
37-
-0.6632096279193985255417942425992805510759 + 8.104093334608059251422673696652054786682im,
38-
-0.420162057582130654687091464438708499074 + 5.169659732326902457089090603403747081757im,
39-
-0.5248407961081728023700065932644065469503 + 6.352576446680355815033180988393723964691im,
40-
16.81912545211429588221108133438974618912 + 11.61870568962804384227638365700840950012im,
41-
3.956891894434294787430417272844351828098 + 2.733954688042517844337453425396233797073im,
42-
0.1577160301669917608080595528008416295052 + 20.49513202616177309778322523925453424454im,
43-
-0.390590921015831871176970935266581363976 + 4.825951828615894889651372068328782916069im,
44-
5.086291973430167701053505879826843738556 + 19.856103426521318766617696383036673069im,
45-
-1.393155126117511022343364857078995555639 + 16.86197099455173997739620972424745559692im,
46-
12.05655555056696748295053112087771296501 + 16.57204880191947538037311460357159376144im,
47-
2.460344943558517627479886868968605995178 + 4.13151392233745973214809055207297205925im,
48-
-1.582576281040411148026691989798564463854 + 20.43576085524474805765748897101730108261im,
49-
8.550688367994505156843842996750026941299 + 5.903012749617017718151146254967898130417im,
50-
-0.9675432524019281776972434272465761750937 + 11.81279027041785312235333549324423074722im,
51-
4.703843537806423391600674221990630030632 + 3.252480256589709295411694256472401320934im,
52-
16.06349030987121295765973627567291259766 + 12.69219331974861475487159623298794031143im,
53-
1.937834775724029512389279261697083711624 + 20.39705296333508854900173901114612817764im,
54-
-1.594679144973751627745173209405038505793 + 19.93238790107540125973173417150974273682im
55-
]
56-
57-
58-
w = [ -0.006546682597891755918395606528292773873545 + 0.0im,
59-
-0.008478994639378135619867116190562228439376 - 0.006865842132444678960756512253738037543371im,
60-
-0.04084417105131761538405754663472180254757 - 0.06528666329071823593022116938300314359367im,
61-
-0.04712461267953402949126839871496486011893 + 0.01452078721700496635738097950252267764881im,
62-
-0.007118804602080300823752079253381452872418 + 0.06492808208351363852273152588168159127235im,
63-
-0.06116738548641747347245356536404869984835 - 0.02725053730346357894198661142581840977073im,
64-
-0.02494387696913849192248413544348295545205 + 0.02277714123083315542195315117623977130279im,
65-
-0.04245963986977926291066509634219983126968 + 0.004116793595813594310028893374919789494015im,
66-
-0.02406420353735235287406801774068298982456 + 0.02715441509989760873744479852121003204957im,
67-
-0.4904880414270129662668296077754348516464 + 0.1081159025127441941638295475058839656413im,
68-
-0.03348344660847365344968906697431521024555 + 0.05785977178570968909587790562909503933042im,
69-
0.003536996149578506415389611561295168939978 + 0.07752104476498027085806796776523697189987im,
70-
0.01586745997824464543546341133151145186275 - 0.01721172123513132687366855577693058876321im,
71-
0.1214650011644559657320030510163633152843 + 0.0199244295229646650735588764291605912149im,
72-
-0.01413089963393250950152157940920005785301 + 0.03915878901517700488854600848753761965781im,
73-
0.1940642102318209105682456083741271868348 + 0.189125716371402241566812563178245909512im,
74-
-0.0629596993796838894086320692622393835336 + 0.03943410994216315496041502797197608742863im,
75-
0.02664896602239549827650932911637937650084 + 0.007480078058428093340515019349368230905384im,
76-
-0.01891162412416614799215430764434131560847 - 0.224558251516702916950052326683362480253im,
77-
-0.07674844370101535639960843582230154424906 + 0.06062885086104936177564539434570178855211im,
78-
-0.05865831566992357748446806908759754151106 - 0.1487332818764729724936302091009565629065im,
79-
0.5977968397140664968958390090847387909889 - 0.3921463933084315400812158713961252942681im,
80-
0.09468179551658635617616965873821754939854 + 0.06801353116428925094094637415764736942947im,
81-
-0.03593292611744278858276757659950817469507 + 0.08129780087330089333175209276305395178497im]
82-
83-
84-
f = [ 3.267431673967855942919413791969418525696 - 11.08592268019023485692287067649886012077im,
85-
-1.408157475121772916892837201885413378477 + 8.194246359208529284501310030464082956314im,
86-
10.16651031219827316931514360476285219193 + 1.139462837357397617665810685139149427414im,
87-
0.05333588002795774246633797588401648681611 + 5.332849960825335244862799299880862236023im,
88-
3.731512177678290687055095986579544842243 + 1.009153918438497665732711539021693170071im,
89-
2.079126843703456906098381296033039689064 - 1.574054723602581251640231130295433104038im,
90-
-0.3027225201326076420293986757314996793866 - 1.336196815913121005436892119178082793951im,
91-
-0.321207149871566732812766531424131244421 - 0.5904585480867294844387060948065482079983im,
92-
0.5332409552090798809942384650639723986387 + 1.12410228905475118033052694954676553607im,
93-
0.4009721205547844280481228906864998862147 - 0.001452802358581500212844628272534919233294im,
94-
0.4076248719737363135351415621698834002018 - 0.007221143389127381538583616560345035395585im,
95-
0.3572930776018199416910192667273804545403 - 0.2904023883125960159290457340830471366644im,
96-
0.2227466624325098176750969969361904077232 - 0.86310886714297674338070009980583563447im,
97-
0.3994990829482118477322671878937399014831 - 0.002392663621286134686266811044674795994069im,
98-
5.223429579451738469231258932268247008324 - 4.319291202018760600367386359721422195435im,
99-
0.4003498255225810820157050784473540261388 - 0.002032103233714888101957285471144132316113im,
100-
0.4060493902079563288687324984493898227811 - 0.01098448458664161297981820553104626014829im,
101-
0.164972906678356817655739519068447407335 - 9.454457188929774602570432762149721384048im,
102-
0.4029739937090559553922730628983117640018 - 0.002997378292125214404445499027929145086091im,
103-
-2.361384545672124968263005939661525189877 + 0.1425123166686921016843569987031514756382im,
104-
0.406335503077652593351132281895843334496 - 0.005903178413394780457701394027481001103297im,
105-
0.4008673043752389864025076349207665771246 - 0.001578422845174482381375158368541633535642im,
106-
0.3995411619962026539276678249734686687589 - 0.01068906314642052712837738681628252379596im,
107-
8.455251786581097661610328941605985164642 - 5.380704447124407430180781375383958220482im
108-
]
109-
110-
111-
C = 1 ./ (z .- transpose(zz))
112-
r = ((C*(w.*f))./(C*w))
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)
113123

114-
r = r[1] / (sqrt(z) * exp(-z))
124+
r = (s1/s2) / (sqrt(z) * exp(-z))
115125
end
116126
end
117127
check && (r = conj(r))
@@ -143,3 +153,6 @@ function k0_large_argument(z::T) where T
143153
return c*p
144154
end
145155
=#
156+
157+
158+
## these are pretty good but still need some work

0 commit comments

Comments
 (0)