Skip to content

Commit 377e08e

Browse files
committed
fix coverage remove 0 branch
1 parent 891deb8 commit 377e08e

File tree

3 files changed

+5
-122
lines changed

3 files changed

+5
-122
lines changed

src/besselj.jl

Lines changed: 3 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,8 @@
88
# where r1 and r2 are zeros of J0
99
# and P3 and Q8 are a 3 and 8 degree polynomial respectively
1010
# Polynomial coefficients are from [1] which is based on [2]
11-
# See [1] for more details and [2] for coefficients of polynomials
11+
# See [1] for more details and [2] for coefficients of polynomials.
12+
# For tiny arugments the power series expansion is used.
1213
#
1314
# Branch 2: 5.0 < x < 80.0
1415
# besselj0 = sqrt(2/(pi*x))*(cos(x - pi/4)*R7(x) - sin(x - pi/4)*R8(x))
@@ -75,7 +76,6 @@ end
7576
function besselj0(x::Float32)
7677
T = Float32
7778
x = abs(x)
78-
iszero(x) && return one(x)
7979
isinf(x) && return zero(x)
8080

8181
if x <= 2.0f0
@@ -100,15 +100,14 @@ end
100100
function besselj1(x::Float64)
101101
T = Float64
102102
x = abs(x)
103-
iszero(x) && return zero(x)
104103
isinf(x) && return zero(x)
105104

106105
if x <= 5.0
107106
z = x * x
108107
w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
109108
w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
110109
return w
111-
elseif x < 100.0
110+
elseif x < 75.0
112111
w = 5.0 / x
113112
z = w * w
114113
p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
@@ -140,7 +139,6 @@ end
140139

141140
function besselj1(x::Float32)
142141
x = abs(x)
143-
iszero(x) && return zero(x)
144142
isinf(x) && return zero(x)
145143

146144
if x <= 2.0f0
@@ -158,66 +156,3 @@ function besselj1(x::Float32)
158156
return p
159157
end
160158
end
161-
162-
function besselj(n::Int, x::Float64)
163-
if n < 0
164-
n = -n
165-
if (n & 1) == 0
166-
sign = 1
167-
else
168-
sign = -1
169-
end
170-
else
171-
sign = 1
172-
end
173-
174-
if x < zero(x)
175-
if (n & 1)
176-
sign = -sign
177-
x = -x
178-
end
179-
end
180-
181-
if n == 0
182-
return sign * besselj0(x)
183-
elseif n == 1
184-
return sign * besselj1(x)
185-
elseif n == 2
186-
return sign * (2.0 * besselj1(x) / x - besselj0(x))
187-
end
188-
189-
#if x < MACHEP
190-
# return 0.0
191-
#end
192-
193-
k = 40 # or 53
194-
pk = 2 * (n + k)
195-
ans = pk
196-
xk = x * x
197-
198-
for _ in 1:k
199-
pk -= 2.0
200-
ans = pk - (xk / ans)
201-
end
202-
203-
ans = x / ans
204-
205-
pk = 1.0
206-
pkm1 = inv(ans)
207-
k = n - 1
208-
r = 2 * k
209-
210-
for _ in 1:k
211-
pkm2 = (pkm1 * r - pk * x) / x
212-
pk = pkm1
213-
pkm1 = pkm2
214-
r -= 2.0
215-
end
216-
if abs(pk) > abs(pkm1)
217-
ans = besselj1(x) / pk
218-
else
219-
ans = besselj0(x) / pkm1
220-
end
221-
222-
return sign * ans
223-
end

src/bessely.jl

Lines changed: 0 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -138,55 +138,3 @@ function _bessely1_compute(x::Float32)
138138
return p
139139
end
140140
end
141-
142-
# Bessel function of second kind, order zero
143-
#
144-
# Bessel function of second kind, order one
145-
#=
146-
Ported to Julia from:
147-
Cephes Math Library Release 2.2: June, 1992
148-
Copyright 1984, 1987, 1992 by Stephen L. Moshier
149-
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
150-
https://github.com/jeremybarnes/cephes/blob/master/single/j0f.c
151-
https://github.com/jeremybarnes/cephes/blob/master/single/j1f.c
152-
=#
153-
154-
#=
155-
function bessely(n::Int, x)
156-
if n < 0
157-
n = -n
158-
if n & 1 == 0
159-
sign = 1
160-
else sign = -1
161-
end
162-
else
163-
sign = 1
164-
end
165-
166-
if n == 0
167-
return sign * bessely0(x)
168-
elseif n == 1
169-
return sign * bessely1(x)
170-
end
171-
172-
if x <= 0.0
173-
return NaN
174-
end
175-
176-
anm2 = bessely0(x)
177-
anm1 = bessely1(x)
178-
an = zero(x)
179-
180-
k = 1
181-
r = 2 * k
182-
183-
for _ in 1:n
184-
an = r * anm1 / x - anm2
185-
anm2 = anm1
186-
anm1 = an
187-
r += 2.0
188-
end
189-
190-
return sign * an
191-
end
192-
=#

src/math_constants.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@ const ONEOSQPI(::Type{BigFloat}) = big"5.641895835477562869480794515607725858440
22
const TWOOPI(::Type{BigFloat}) = big"6.3661977236758134307553505349005744813784E-1"
33
#const SQPIO2(::Type{BigFloat}) = big"1.253314137315500251207882642405522626503493370304969158314961788171146827303924"
44
#const SQ1O2PI(::Type{BigFloat}) = big"0.3989422804014326779399460599343818684758586311649346576659258296706579258993008"
5-
const SQ2OPI(::Type{BigFloat}) = big"0.7978845608028653558798921198687637369517172623298693153318516593413158517986017"
6-
const PIO4(::Type{BigFloat}) = big"0.7853981633974483096156608458198757210492923498437764552437361480769541015715495"
5+
#const SQ2OPI(::Type{BigFloat}) = big"0.7978845608028653558798921198687637369517172623298693153318516593413158517986017"
6+
#const PIO4(::Type{BigFloat}) = big"0.7853981633974483096156608458198757210492923498437764552437361480769541015715495"
77

88
const PIO4(::Type{Float64}) = .78539816339744830962
99
const THPIO4(::Type{Float64}) = 2.35619449019234492885

0 commit comments

Comments
 (0)