Skip to content

Commit f6c594b

Browse files
heltonmcheltonmc
committed
add tests
1 parent 15f9f3e commit f6c594b

File tree

5 files changed

+136
-108
lines changed

5 files changed

+136
-108
lines changed

‎src/besseli.jl

Lines changed: 114 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@
5151
# [3] https://dlmf.nist.gov/10.41
5252

5353
"""
54-
besseli0(x::T) where T <: Union{Float32, Float64}
54+
besseli0(x::T) where T <: Union{Float32, Float64, ComplexF32, ComplexF64}
5555
5656
Modified Bessel function of the first kind of order zero, ``I_0(x)``.
5757
"""
@@ -130,6 +130,119 @@ function besseli1x(x::T) where T <: Union{Float32, Float64}
130130
return z
131131
end
132132

133+
#####
134+
##### Implementation for complex arguments
135+
#####
136+
137+
function besseli0(z::ComplexF64)
138+
check = false
139+
if real(z) < 0.0
140+
z = -z
141+
end
142+
if imag(z) < 0.0
143+
z = conj(z)
144+
check = true
145+
end
146+
147+
if abs(z) > 17.5
148+
zinv = 1 / z
149+
e = exp(z)
150+
s = sqrt(2 * z * π)
151+
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))
152+
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))
153+
r = e / s * (p + p2) + im * (p - p2) / (e * s)
154+
elseif abs(z) < 5.2
155+
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))
156+
else
157+
if angle(z) <= π / 4.4
158+
zz = z*z
159+
zz2 = zz*zz
160+
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))
161+
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))
162+
else
163+
zz = (-1.3144540816651462 + 20.054693865816724im, -1.2152720220184496 + 18.54146805523844im,
164+
-1.0758582968968085 + 16.41442564500392im, -0.9186088936209557 + 14.015263371275323im,
165+
-0.7533475204960969 + 11.49386205943562im, -0.5838248556039682 + 8.907445998843897im,
166+
-0.41092161465730614 + 6.269452314652046im, -0.32055227293065935 + 4.890682596894067im,
167+
16.906716750808087 + 10.865287107782516im, 0.26715371754140166 + 20.098224520867603im,
168+
4.121994532242284 + 2.6493699394695174im, 4.053028754723365 + 19.68712670537236im,
169+
-0.5018982035099 + 7.657486833198153im, 9.416081350634247 + 17.758023876497013im,
170+
1.9674895590143007 + 4.487648029332259im, -1.275759132467892 + 19.46432302583941im,
171+
7.709313889965882 + 4.954475197822861im, 13.72022229845232 + 14.68895844098729im,
172+
-0.3541218766747901 + 5.402855776372861im, -0.9973883749616251 + 15.217205990064667im,
173+
5.067285183070908 + 3.256546447342955im, 2.1309472708925585 + 19.986722185708082im,
174+
-0.8162780120038937 + 12.453995821138065im, -0.6642234573493526 + 10.134091621337502im
175+
)
176+
w = (0.023435912479136792 + 0.0im, -0.07518141047674767 + 0.0475948197731169im,
177+
0.02099862365431344 + 0.10227799317706024im, -0.044790606134137025 - 0.0019900581421525287im,
178+
0.05092855107420865 - 0.030284027691088355im, 0.11876822116350058 + 0.3545737241031229im,
179+
-0.2800428715703842 + 0.1385502878368521im, 0.07056156978707125 - 0.015640049217372977im,
180+
-0.06946091028316168 + 0.016066454632345423im, 0.08524646892263366 + 0.0953622756709394im,
181+
-0.02790219158045262 + 0.07064926802304639im, 0.1422314233425249 - 0.2536760604980481im,
182+
-0.20586363442913913 + 0.3014383126557577im, 0.20547145538054265 - 0.17388658146749558im,
183+
-0.13839690210608482 - 0.07997812132261209im, -0.005116954351018287 + 0.07859682420061695im,
184+
0.0409490985507271 - 0.2642881024594653im, -0.04174610452138103 - 0.23810783641445965im,
185+
-0.14078210159561536 - 0.17647153499232213im, -0.07672294420555978 + 0.08789101702567167im,
186+
-0.15544491331665283 - 0.17872134246439533im, 0.2746212867658749 + 0.0633913060138901im,
187+
0.003561860807958012 + 0.005781000943183686im, 0.22467614820711904 + 0.050873741433048215im
188+
)
189+
f = (4.11757504005866 - 4.095512734321109im, -2.247112620034342 + 3.680759784557969im,
190+
3.7834318298173315 + 0.562243565983836im, 1.025044899237516 - 2.4295763517152302im,
191+
-1.1024728966418114 - 0.9985456859791928im, -0.7134052191604513 + 0.634168698676205im,
192+
0.35489432819051814 + 0.8984965319188546im, 0.15105140195420463 - 0.7254114647779558im,
193+
0.401057852147248 - 0.0014084695668244449im, 0.5404571564588408 - 0.1885021348797114im,
194+
0.407770977742567 - 0.00680290179701475im, 0.3994965367795384 - 0.0024661141779376626im,
195+
0.8312271064319111 - 1.0050503936637845im, 0.40006191058465196 - 0.002248719754950748im,
196+
0.4055260568863696 - 0.01689191195791546im, 5.210832430161148 + 1.7413278357835171im,
197+
0.40364747476765517 - 0.003288380266379im, 0.40062798907705705 - 0.0018845654535970921im,
198+
-0.39383201186515915 - 0.17987628384333437im, -2.053184226968296 + 1.6069066509951657im,
199+
0.40617427977861376 - 0.005320233415147144im, 0.4034434765080089 - 0.006086635360116373im,
200+
-0.07645988937623527 + 1.981923553555365im, 1.8844134122620073 + 0.24179792442042025im
201+
).*w
202+
203+
s1 = 0.0 + 0.0im
204+
s2 = 0.0 + 0.0im
205+
@fastmath for ind in eachindex(f)
206+
C = inv(z - zz[ind])
207+
s1 += C*f[ind]
208+
s2 += C*w[ind]
209+
end
210+
211+
r = s1 / (s2 * sqrt(z) * exp(-z))
212+
end
213+
end
214+
check && (r = conj(r))
215+
return r
216+
end
217+
218+
function besseli0(z::ComplexF32)
219+
z = ComplexF64(z)
220+
check = false
221+
if real(z) < 0.0
222+
z = -z
223+
end
224+
if imag(z) < 0.0
225+
z = conj(z)
226+
check = true
227+
end
228+
229+
if abs(z) > 10.0
230+
zinv = 1 / z
231+
e = exp(z)
232+
s = sqrt(2 * z * π)
233+
p = evalpoly(zinv*zinv, (1.0, 0.0703125, 0.112152099609375, 0.5725014209747314, 6.074042001273483, 110.01714026924674))
234+
p2 = zinv*evalpoly(zinv*zinv, (0.125, 0.0732421875, 0.22710800170898438, 1.7277275025844574, 24.380529699556064, 551.3358961220206))
235+
r = e / s * (p + p2) + im * (p - p2) / (e * s)
236+
else
237+
zz = z*z
238+
zz2 = zz*zz
239+
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))
240+
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))
241+
end
242+
check && (r = conj(r))
243+
return ComplexF32(r)
244+
end
245+
133246
# Modified Bessel functions of the first kind of order nu
134247
# besseli(nu, x)
135248
#

‎src/besselj.jl

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
#
2525

2626
"""
27-
besselj0(x::T) where T <: Union{Float32, Float64}
27+
besselj0(x::T) where T <: Union{Float32, Float64, ComplexF32, ComplexF64}
2828
2929
Bessel function of the first kind of order zero, ``J_0(x)``.
3030
"""
@@ -173,6 +173,12 @@ function _besselj1(x::Float32)
173173
end
174174
end
175175

176+
#####
177+
##### Implementation for complex arguments
178+
#####
179+
180+
besselj0(z::T) where T <: Union{ComplexF32, ComplexF64} = besseli0(z/im)
181+
176182
# Bessel functions of the first kind of order nu
177183
# besselj(nu, x)
178184
#

‎src/i0.jl

Lines changed: 0 additions & 106 deletions
This file was deleted.

0 commit comments

Comments
 (0)