Skip to content

Commit 1d41167

Browse files
heltonmcheltonmc
authored
Merge pull request #13 from heltonmc/besselj0_largeargs
Better large order approximations for J0,J1, Y0,Y1
2 parents 1a793d9 + c953d2b commit 1d41167

File tree

7 files changed

+181
-150
lines changed

7 files changed

+181
-150
lines changed

‎src/Bessels.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@ include("Float128/constants.jl")
3636
include("math_constants.jl")
3737
include("U_polynomials.jl")
3838
include("recurrence.jl")
39+
include("misc.jl")
3940
#include("hankel.jl")
4041

4142
end

‎src/Float128/besselj.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,4 @@ function besselj0(x::BigFloat)
6767
return z
6868
end
6969
end
70+

‎src/besselj.jl

Lines changed: 69 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -1,42 +1,77 @@
1-
#=
2-
Cephes Math Library Release 2.8: June, 2000
3-
Copyright 1984, 1987, 2000 by Stephen L. Moshier
4-
https://github.com/jeremybarnes/cephes/blob/master/bessel/j0.c
5-
https://github.com/jeremybarnes/cephes/blob/master/bessel/j1.c
6-
=#
7-
function besselj0(x::Float64)
8-
T = Float64
1+
# Bessel functions of the first kind of order zero and one
2+
# besselj0, besselj1
3+
#
4+
# Calculation of besselj0 is done in three branches using polynomial approximations
5+
#
6+
# Branch 1: x <= 5.0
7+
# besselj0 = (x^2 - r1^2)*(x^2 - r2^2)*P3(x^2) / Q8(x^2)
8+
# where r1 and r2 are zeros of J0
9+
# and P3 and Q8 are a 3 and 8 degree polynomial respectively
10+
# Polynomial coefficients are from [1] which is based on [2]
11+
# For tiny arugments the power series expansion is used.
12+
#
13+
# Branch 2: 5.0 < x < 75.0
14+
# besselj0 = sqrt(2/(pi*x))*(cos(x - pi/4)*R7(x) - sin(x - pi/4)*R8(x))
15+
# Hankel's asymptotic expansion is used
16+
# where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
17+
# See section 4 of [3] for more details and [1] for coefficients of polynomials
18+
#
19+
# Branch 3: x >= 75.0
20+
# besselj0 = sqrt(2/(pi*x))*beta(x)*(cos(x - pi/4 - alpha(x))
21+
# See modified expansions given in [3]. Exact coefficients are used
22+
#
23+
# Calculation of besselj1 is done in a similar way as besselj0.
24+
# See [3] for details on similarities.
25+
#
26+
# [1] https://github.com/deepmind/torch-cephes
27+
# [2] Cephes Math Library Release 2.8: June, 2000 by Stephen L. Moshier
28+
# [3] Harrison, John. "Fast and accurate Bessel function computation."
29+
# 2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
30+
#
31+
function besselj0(x::T) where T
932
x = abs(x)
10-
iszero(x) && return one(x)
1133
isinf(x) && return zero(x)
1234

1335
if x <= 5
1436
z = x * x
1537
if x < 1.0e-5
1638
return 1.0 - z / 4.0
1739
end
18-
19-
DR1 = 5.78318596294678452118E0
20-
DR2 = 3.04712623436620863991E1
21-
40+
DR1 = 5.78318596294678452118e0
41+
DR2 = 3.04712623436620863991e1
2242
p = (z - DR1) * (z - DR2)
2343
p = p * evalpoly(z, RP_j0(T)) / evalpoly(z, RQ_j0(T))
2444
return p
25-
else
45+
elseif x < 75.0
2646
w = 5.0 / x
2747
q = 25.0 / (x * x)
2848

2949
p = evalpoly(q, PP_j0(T)) / evalpoly(q, PQ_j0(T))
3050
q = evalpoly(q, QP_j0(T)) / evalpoly(q, QQ_j0(T))
3151
xn = x - PIO4(T)
32-
p = p * cos(xn) - w * q * sin(xn)
52+
sc = sincos(xn)
53+
p = p * sc[2] - w * q * sc[1]
3354
return p * SQ2OPI(T) / sqrt(x)
55+
else
56+
xinv = inv(x)
57+
x2 = xinv*xinv
58+
59+
p = (one(T), -1/16, 53/512, -4447/8192, 5066403/524288)
60+
p = evalpoly(x2, p)
61+
a = SQ2OPI(T) * sqrt(xinv) * p
62+
63+
q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296)
64+
xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))
65+
66+
# the following computes b = cos(x + xn) more accurately
67+
# see src/misc.jl
68+
b = cos_sum(x, xn)
69+
return a * b
3470
end
3571
end
3672
function besselj0(x::Float32)
3773
T = Float32
3874
x = abs(x)
39-
iszero(x) && return one(x)
4075
isinf(x) && return zero(x)
4176

4277
if x <= 2.0f0
@@ -61,28 +96,42 @@ end
6196
function besselj1(x::Float64)
6297
T = Float64
6398
x = abs(x)
64-
iszero(x) && return zero(x)
6599
isinf(x) && return zero(x)
66100

67101
if x <= 5.0
68102
z = x * x
69103
w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
70104
w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
71105
return w
72-
else
106+
elseif x < 75.0
73107
w = 5.0 / x
74108
z = w * w
75109
p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
76110
q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
77111
xn = x - THPIO4(T)
78-
p = p * cos(xn) - w * q * sin(xn)
112+
sc = sincos(xn)
113+
p = p * sc[2] - w * q * sc[1]
79114
return p * SQ2OPI(T) / sqrt(x)
115+
else
116+
xinv = inv(x)
117+
x2 = xinv*xinv
118+
119+
p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288)
120+
p = evalpoly(x2, p)
121+
a = SQ2OPI(T) * sqrt(xinv) * p
122+
123+
q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144)
124+
xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
125+
126+
# the following computes b = cos(x + xn) more accurately
127+
# see src/misc.jl
128+
b = cos_sum(x, xn)
129+
return a * b
80130
end
81131
end
82132

83133
function besselj1(x::Float32)
84134
x = abs(x)
85-
iszero(x) && return zero(x)
86135
isinf(x) && return zero(x)
87136

88137
if x <= 2.0f0
@@ -100,66 +149,3 @@ function besselj1(x::Float32)
100149
return p
101150
end
102151
end
103-
104-
function besselj(n::Int, x::Float64)
105-
if n < 0
106-
n = -n
107-
if (n & 1) == 0
108-
sign = 1
109-
else
110-
sign = -1
111-
end
112-
else
113-
sign = 1
114-
end
115-
116-
if x < zero(x)
117-
if (n & 1)
118-
sign = -sign
119-
x = -x
120-
end
121-
end
122-
123-
if n == 0
124-
return sign * besselj0(x)
125-
elseif n == 1
126-
return sign * besselj1(x)
127-
elseif n == 2
128-
return sign * (2.0 * besselj1(x) / x - besselj0(x))
129-
end
130-
131-
#if x < MACHEP
132-
# return 0.0
133-
#end
134-
135-
k = 40 # or 53
136-
pk = 2 * (n + k)
137-
ans = pk
138-
xk = x * x
139-
140-
for _ in 1:k
141-
pk -= 2.0
142-
ans = pk - (xk / ans)
143-
end
144-
145-
ans = x / ans
146-
147-
pk = 1.0
148-
pkm1 = inv(ans)
149-
k = n - 1
150-
r = 2 * k
151-
152-
for _ in 1:k
153-
pkm2 = (pkm1 * r - pk * x) / x
154-
pk = pkm1
155-
pkm1 = pkm2
156-
r -= 2.0
157-
end
158-
if abs(pk) > abs(pkm1)
159-
ans = besselj1(x) / pk
160-
else
161-
ans = besselj0(x) / pkm1
162-
end
163-
164-
return sign * ans
165-
end

‎src/bessely.jl

Lines changed: 67 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,33 @@
1-
#=
2-
Cephes Math Library Release 2.8: June, 2000
3-
Copyright 1984, 1987, 2000 by Stephen L. Moshier
4-
https://github.com/jeremybarnes/cephes/blob/master/bessel/j0.c
5-
https://github.com/jeremybarnes/cephes/blob/master/bessel/j1.c
6-
=#
1+
# Bessel functions of the second kind of order zero and one
2+
# bessely0, bessely1
3+
#
4+
# Calculation of bessely0 is done in three branches using polynomial approximations
5+
#
6+
# Branch 1: x <= 5.0
7+
# bessely0 = R(x^2) + 2*log(x)*besselj0(x) / pi
8+
# where r1 and r2 are zeros of J0
9+
# and P3 and Q8 are a 3 and 8 degree polynomial respectively
10+
# Polynomial coefficients are from [1] which is based on [2].
11+
# For tiny arugments the power series expansion is used.
12+
#
13+
# Branch 2: 5.0 < x < 75.0
14+
# bessely0 = sqrt(2/(pi*x))*(sin(x - pi/4)*R7(x) - cos(x - pi/4)*R8(x))
15+
# Hankel's asymptotic expansion is used
16+
# where R7 and R8 are rational functions (Pn(x)/Qn(x)) of degree 7 and 8 respectively
17+
# See section 4 of [3] for more details and [1] for coefficients of polynomials
18+
#
19+
# Branch 3: x >= 75.0
20+
# bessely0 = sqrt(2/(pi*x))*beta(x)*(sin(x - pi/4 - alpha(x))
21+
# See modified expansions given in [3]. Exact coefficients are used.
22+
#
23+
# Calculation of bessely1 is done in a similar way as bessely0.
24+
# See [3] for details on similarities.
25+
#
26+
# [1] https://github.com/deepmind/torch-cephes
27+
# [2] Cephes Math Library Release 2.8: June, 2000 by Stephen L. Moshier
28+
# [3] Harrison, John. "Fast and accurate Bessel function computation."
29+
# 2009 19th IEEE Symposium on Computer Arithmetic. IEEE, 2009.
30+
#
731
function bessely0(x::T) where T <: Union{Float32, Float64}
832
if x <= zero(x)
933
if iszero(x)
@@ -23,14 +47,30 @@ function _bessely0_compute(x::Float64)
2347
w = evalpoly(z, YP_y0(T)) / evalpoly(z, YQ_y0(T))
2448
w += TWOOPI(T) * log(x) * besselj0(x)
2549
return w
26-
else
50+
elseif x < 75.0
2751
w = T(5) / x
2852
z = w*w
2953
p = evalpoly(z, PP_y0(T)) / evalpoly(z, PQ_y0(T))
3054
q = evalpoly(z, QP_y0(T)) / evalpoly(z, QQ_y0(T))
3155
xn = x - PIO4(T)
32-
p = p * sin(xn) + w * q * cos(xn);
56+
sc = sincos(xn)
57+
p = p * sc[1] + w * q * sc[2]
3358
return p * SQ2OPI(T) / sqrt(x)
59+
else
60+
xinv = inv(x)
61+
x2 = xinv*xinv
62+
63+
p = (one(T), -1/16, 53/512, -4447/8192, 5066403/524288)
64+
p = evalpoly(x2, p)
65+
a = SQ2OPI(T) * sqrt(xinv) * p
66+
67+
q = (-1/8, 25/384, -1073/5120, 375733/229376, -55384775/2359296)
68+
xn = muladd(xinv, evalpoly(x2, q), - PIO4(T))
69+
70+
# the following computes b = sin(x + xn) more accurately
71+
# see src/misc.jl
72+
b = sin_sum(x, xn)
73+
return a * b
3474
end
3575
end
3676
function _bessely0_compute(x::Float32)
@@ -71,22 +111,38 @@ function _bessely1_compute(x::Float64)
71111
w = x * (evalpoly(z, YP_y1(T)) / evalpoly(z, YQ_y1(T)))
72112
w += TWOOPI(T) * (besselj1(x) * log(x) - inv(x))
73113
return w
74-
else
114+
elseif x < 75.0
75115
w = T(5) / x
76116
z = w * w
77117
p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
78118
q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
79119
xn = x - THPIO4(T)
80-
p = p * sin(xn) + w * q * cos(xn)
120+
sc = sincos(xn)
121+
p = p * sc[1] + w * q * sc[2]
81122
return p * SQ2OPI(T) / sqrt(x)
123+
else
124+
xinv = inv(x)
125+
x2 = xinv*xinv
126+
127+
p = (one(T), 3/16, -99/512, 6597/8192, -4057965/524288)
128+
p = evalpoly(x2, p)
129+
a = SQ2OPI(T) * sqrt(xinv) * p
130+
131+
q = (3/8, -21/128, 1899/5120, -543483/229376, 8027901/262144)
132+
xn = muladd(xinv, evalpoly(x2, q), - 3 * PIO4(T))
133+
134+
# the following computes b = sin(x + xn) more accurately
135+
# see src/misc.jl
136+
b = sin_sum(x, xn)
137+
return a * b
82138
end
83139
end
84140

85141
function _bessely1_compute(x::Float32)
86142
T = Float32
87143
if x <= 2.0f0
88144
z = x * x
89-
YO1 = 4.66539330185668857532f0
145+
YO1 = 4.66539330185668857532f0
90146
w = (z - YO1) * x * evalpoly(z, YP32)
91147
w += TWOOPI(Float32) * (besselj1(x) * log(x) - inv(x))
92148
return w
@@ -100,55 +156,3 @@ function _bessely1_compute(x::Float32)
100156
return p
101157
end
102158
end
103-
104-
# Bessel function of second kind, order zero
105-
#
106-
# Bessel function of second kind, order one
107-
#=
108-
Ported to Julia from:
109-
Cephes Math Library Release 2.2: June, 1992
110-
Copyright 1984, 1987, 1992 by Stephen L. Moshier
111-
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
112-
https://github.com/jeremybarnes/cephes/blob/master/single/j0f.c
113-
https://github.com/jeremybarnes/cephes/blob/master/single/j1f.c
114-
=#
115-
116-
#=
117-
function bessely(n::Int, x)
118-
if n < 0
119-
n = -n
120-
if n & 1 == 0
121-
sign = 1
122-
else sign = -1
123-
end
124-
else
125-
sign = 1
126-
end
127-
128-
if n == 0
129-
return sign * bessely0(x)
130-
elseif n == 1
131-
return sign * bessely1(x)
132-
end
133-
134-
if x <= 0.0
135-
return NaN
136-
end
137-
138-
anm2 = bessely0(x)
139-
anm1 = bessely1(x)
140-
an = zero(x)
141-
142-
k = 1
143-
r = 2 * k
144-
145-
for _ in 1:n
146-
an = r * anm1 / x - anm2
147-
anm2 = anm1
148-
anm1 = an
149-
r += 2.0
150-
end
151-
152-
return sign * an
153-
end
154-
=#

0 commit comments

Comments
 (0)