Skip to content

Commit e6c7269

Browse files
matthew-wingateararslan
matthew-wingate
authored and
ararslan
committed
Correct besseli(nu, 0) & besselix(nu, 0) with nu<0 and nu an Int. (#34)
* Correct besseli(nu, 0) & besselix(nu, 0) with nu<0 and nu an Int. An AmosException is thrown for besseli(nu, 0) or besselix(nu, 0) due to the call to _besselk. However, for integer nu, besseli(nu, 0) = besseli(-nu,0) = 0 and similar for besselix. A test is added before calling besselk. * Correct test for AmosException in besselix(nu, z) besselix(nu, z) is divergent for z=0 and negative nu, only if nu is not an integer. * Correct besselj and besseljx for nu<0 and nu an Int. * besselj(nu, 0) and besseljx(nu, 0) should equal 0 for integer nu!=0. AmosException was thrown for z=complex(0). * Tests corrected and more added to runtest.jl * Changed `if sinpi(nu) == 0` to `if isinteger(nu)` * Additions to runtests, parens in sincosint, hope to increase coverage * Attempt to get coverage test past at-horner macro.
1 parent c7e0993 commit e6c7269

File tree

3 files changed

+70
-24
lines changed

3 files changed

+70
-24
lines changed

‎src/bessel.jl

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -272,23 +272,35 @@ end
272272

273273
function besseli(nu::Float64, z::Complex128)
274274
if nu < 0
275-
return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi
275+
if isinteger(nu)
276+
return _besseli(-nu,z,Int32(1))
277+
else
278+
return _besseli(-nu,z,Int32(1)) - 2_besselk(-nu,z,Int32(1))*sinpi(nu)/pi
279+
end
276280
else
277281
return _besseli(nu,z,Int32(1))
278282
end
279283
end
280284

281285
function besselix(nu::Float64, z::Complex128)
282286
if nu < 0
283-
return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
287+
if isinteger(nu)
288+
return _besseli(-nu,z,Int32(2))
289+
else
290+
return _besseli(-nu,z,Int32(2)) - 2_besselk(-nu,z,Int32(2))*exp(-abs(real(z))-z)*sinpi(nu)/pi
291+
end
284292
else
285293
return _besseli(nu,z,Int32(2))
286294
end
287295
end
288296

289297
function besselj(nu::Float64, z::Complex128)
290298
if nu < 0
291-
return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu)
299+
if isinteger(nu)
300+
return _besselj(-nu,z,Int32(1))*cospi(nu)
301+
else
302+
return _besselj(-nu,z,Int32(1))*cospi(nu) + _bessely(-nu,z,Int32(1))*sinpi(nu)
303+
end
292304
else
293305
return _besselj(nu,z,Int32(1))
294306
end
@@ -300,7 +312,11 @@ besselj(nu::Cint, x::Float32) = ccall((:jnf, libm), Float32, (Cint, Float32), nu
300312

301313
function besseljx(nu::Float64, z::Complex128)
302314
if nu < 0
303-
return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu)
315+
if isinteger(nu)
316+
return _besselj(-nu,z,Int32(2))*cospi(nu)
317+
else
318+
return _besselj(-nu,z,Int32(2))*cospi(nu) + _bessely(-nu,z,Int32(2))*sinpi(nu)
319+
end
304320
else
305321
return _besselj(nu,z,Int32(2))
306322
end

‎src/sincosint.jl

Lines changed: 21 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ using Base.Math.@horner
1818
function sinint(x::Float64)
1919
t = x*x
2020
if t ≤ 36.0
21-
return x * @horner(t, 1.00000000000000000000E0,
21+
y = (x * @horner(t, 1.00000000000000000000E0,
2222
-0.44663998931312457298E-1,
2323
0.11209146443112369449E-2,
2424
-0.13276124407928422367E-4,
@@ -33,10 +33,11 @@ function sinint(x::Float64)
3333
0.54747121846510390750E-9,
3434
0.10378561511331814674E-11,
3535
0.13754880327250272679E-14,
36-
0.10223981202236205703E-17)
36+
0.10223981202236205703E-17))
37+
return y
3738
elseif t ≤ 144.0
3839
invt = inv(t)
39-
return copysign(π/2, x) - cos(x) *
40+
y = (copysign(π/2, x) - cos(x) *
4041
@horner(invt, 0.99999999962173909991E0,
4142
0.36451060338631902917E3,
4243
0.44218548041288440874E5,
@@ -70,10 +71,11 @@ function sinint(x::Float64)
7071
0.62273134702439012114E10,
7172
0.54570971054996441467E11,
7273
0.18241750166645704670E12,
73-
0.15407148148861454434E12)
74+
0.15407148148861454434E12))
75+
return y
7476
elseif t < Inf
7577
invt = inv(t)
76-
return copysign(π/2, x) - cos(x) / x * (1.0 -
78+
y = (copysign(π/2, x) - cos(x) / x * (1.0 -
7779
@horner(invt, 0.19999999999999978257E1,
7880
0.22206119380434958727E4,
7981
0.84749007623988236808E6,
@@ -107,7 +109,8 @@ function sinint(x::Float64)
107109
0.26028585666152144496E13,
108110
0.85134283716950697226E14,
109111
0.11304079361627952930E16,
110-
0.42519841479489798424E16)*invt)
112+
0.42519841479489798424E16)*invt))
113+
return y
111114
elseif isnan(x)
112115
return NaN
113116
else
@@ -122,7 +125,7 @@ function cosint(x::Float64)
122125
if x < 0.0
123126
throw(DomainErrorNoArgs)
124127
elseif x ≤ 3.0
125-
return log(x/r0) + ((x - r01) - r02) * (x + r0) *
128+
y = (log(x/r0) + ((x - r01) - r02) * (x + r0) *
126129
@horner(t, -0.24607411378767540707E0,
127130
0.72113492241301534559E-2,
128131
-0.11867127836204767056E-3,
@@ -134,9 +137,10 @@ function cosint(x::Float64)
134137
0.78168450570724148921E-4,
135138
0.29959200177005821677E-6,
136139
0.73191677761328838216E-9,
137-
0.94351174530907529061E-12)
140+
0.94351174530907529061E-12))
141+
return y
138142
elseif x ≤ 6.0
139-
return log(x/r1) + ((x - r11) - r12) * (x + r1) *
143+
y = (log(x/r1) + ((x - r11) - r12) * (x + r1) *
140144
@horner(t, -0.15684781827145408780E0,
141145
0.66253165609605468916E-2,
142146
-0.12822297297864512864E-3,
@@ -151,10 +155,11 @@ function cosint(x::Float64)
151155
0.13544922659627723233E-6,
152156
0.27715365686570002081E-9,
153157
0.37718676301688932926E-12,
154-
0.27706844497155995398E-15)
158+
0.27706844497155995398E-15))
159+
return y
155160
elseif x ≤ 12.0
156161
invt = inv(t)
157-
return sin(x) * @horner(invt, 0.99999999962173909991E0,
162+
y = (sin(x) * @horner(invt, 0.99999999962173909991E0,
158163
0.36451060338631902917E3,
159164
0.44218548041288440874E5,
160165
0.22467569405961151887E7,
@@ -187,10 +192,11 @@ function cosint(x::Float64)
187192
0.62273134702439012114E10,
188193
0.54570971054996441467E11,
189194
0.18241750166645704670E12,
190-
0.15407148148861454434E12)
195+
0.15407148148861454434E12))
196+
return y
191197
elseif x < Inf
192198
invt = inv(t)
193-
return sin(x)/x * (1.0 - @horner(invt, 0.19999999999999978257E1,
199+
y = (sin(x)/x * (1.0 - @horner(invt, 0.19999999999999978257E1,
194200
0.22206119380434958727E4,
195201
0.84749007623988236808E6,
196202
0.13959267954823943232E9,
@@ -223,7 +229,8 @@ function cosint(x::Float64)
223229
0.26028585666152144496E13,
224230
0.85134283716950697226E14,
225231
0.11304079361627952930E16,
226-
0.42519841479489798424E16)*invt)
232+
0.42519841479489798424E16)*invt))
233+
return y
227234
elseif isnan(x)
228235
return NaN
229236
else

‎test/runtests.jl

Lines changed: 29 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,9 @@ end
8686
@test_throws ErrorException SF.sinint(big(1.0))
8787
@test_throws ErrorException SF.cosint(big(1.0))
8888
@test_throws DomainError SF.cosint(-1.0)
89+
@test_throws DomainError SF.cosint(Float32(-1.0))
90+
@test_throws DomainError SF.cosint(Float16(-1.0))
91+
@test_throws DomainError SF.cosint(-1//2)
8992
end
9093

9194
@testset "airy" begin
@@ -161,6 +164,13 @@ end
161164
@test SF.besseli(3,-3) ≈ -true_i33
162165
@test SF.besseli(-3,-3) ≈ -true_i33
163166
@test SF.besseli(Float32(-3),Complex64(-3,0)) ≈ -true_i33
167+
true_im3p1_3 = 0.84371226532586351965
168+
@test SF.besseli(-3.1,3) ≈ true_im3p1_3
169+
for i in [-5 -3 -1 1 3 5]
170+
@test SF.besseli(i,0) == 0.0
171+
@test SF.besseli(i,Float32(0)) == 0
172+
@test SF.besseli(i,Complex64(0)) == 0
173+
end
164174
@testset "Error throwing" begin
165175
@test_throws SF.AmosException SF.besseli(1,1000)
166176
@test_throws DomainError SF.besseli(0.4,-1.0)
@@ -172,11 +182,10 @@ end
172182
end
173183
@testset "besselj" begin
174184
@test SF.besselj(0,0) == 1
175-
for i = 1:5
185+
for i in [-5 -3 -1 1 3 5]
176186
@test SF.besselj(i,0) == 0
177-
@test SF.besselj(-i,0) == 0
178-
@test SF.besselj(-i,Float32(0)) == 0
179-
@test SF.besselj(-i,Float32(0)) == 0
187+
@test SF.besselj(i,Float32(0)) == 0
188+
@test SF.besselj(i,Complex64(0)) == 0
180189
end
181190

182191
j33 = SF.besselj(3,3.)
@@ -205,6 +214,10 @@ end
205214
@test SF.besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im
206215
@test SF.besselj(1, 3im) ≈ 3.953370217402609396im
207216
@test SF.besselj(1.0,3im) ≈ SF.besselj(1,3im)
217+
218+
true_jm3p1_3 = -0.45024252862270713882
219+
@test SF.besselj(-3.1,3) ≈ true_jm3p1_3
220+
208221
@testset "Error throwing" begin
209222
@test_throws DomainError SF.besselj(0.1, -0.4)
210223
@test_throws SF.AmosException SF.besselj(20,1000im)
@@ -245,6 +258,8 @@ end
245258
@test_throws DomainError SF.bessely(0.4,-1.0)
246259
@test_throws DomainError SF.bessely(0.4,Float32(-1.0))
247260
@test_throws DomainError SF.bessely(1,Float32(-1.0))
261+
@test_throws DomainError SF.bessely(0.4,BigFloat(-1.0))
262+
@test_throws DomainError SF.bessely(1,BigFloat(-1.0))
248263
@test_throws DomainError SF.bessely(Cint(3),Float32(-3.))
249264
@test_throws DomainError SF.bessely(Cint(3),Float64(-3.))
250265

@@ -279,11 +294,19 @@ end
279294
z == zero(z) || @test SF.besselyx(nu, z) ≈ SF.bessely(nu, z) * exp(-abs(imag(z)))
280295
end
281296
@test SF.besselkx(1, 0) == Inf
297+
for i = [-5 -3 -1 1 3 5]
298+
@test SF.besseljx(i,0) == 0
299+
@test SF.besselix(i,0) == 0
300+
@test SF.besseljx(i,Float32(0)) == 0
301+
@test SF.besselix(i,Float32(0)) == 0
302+
@test SF.besseljx(i,Complex64(0)) == 0
303+
@test SF.besselix(i,Complex64(0)) == 0
304+
end
282305
@testset "Error throwing" begin
283306
@test_throws SF.AmosException SF.hankelh1x(1, 0)
284307
@test_throws SF.AmosException SF.hankelh2x(1, 0)
285-
@test_throws SF.AmosException SF.besselix(-1, 0)
286-
@test_throws SF.AmosException SF.besseljx(-1, 0)
308+
@test_throws SF.AmosException SF.besselix(-1.01, 0)
309+
@test_throws SF.AmosException SF.besseljx(-1.01, 0)
287310
@test_throws SF.AmosException SF.besselyx(1, 0)
288311
@test_throws DomainError SF.besselix(0.4,-1.0)
289312
@test_throws DomainError SF.besseljx(0.4, -1.0)

0 commit comments

Comments
 (0)