Skip to content

Commit 608efb8

Browse files
oscardssmithoscardssmith
committed
merge Float32 and Float64 impls
1 parent 104b6aa commit 608efb8

File tree

5 files changed

+710
-0
lines changed

5 files changed

+710
-0
lines changed

‎src/besseli.jl

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
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/i0.c
5+
https://github.com/jeremybarnes/cephes/blob/master/bessel/i1.c
6+
=#
7+
function besseli0(x::T) where T <: Union{Float32, Float64}
8+
if x <= 8
9+
y = muladd(x, T(.5), T(-2))
10+
return exp(x) * chbevl(y, A_i0(T))
11+
else
12+
return exp(x) * chbevl(T(32) / x - T(-2), B_i0(T)) / sqrt(x)
13+
end
14+
end
15+
function besseli0x(x::T) where T <: Union{Float32, Float64}
16+
T = Float64
17+
x = abs(x)
18+
if x <= 8
19+
y = muladd(x, T(.5), T(-2))
20+
return chbevl(y, A_i0(T))
21+
else
22+
return chbevl(T(32) / x - T(-2), B_i0(T)) / sqrt(x)
23+
end
24+
end
25+
function besseli1(x::T) where T <: Union{Float32, Float64}
26+
z = abs(x)
27+
if x <= 8
28+
y = muladd(x, T(.5), T(-2))
29+
z = chbevl(y, A_i1(T)) * z * exp(z)
30+
else
31+
z = exp(z) * chbevl(T(32) / x - T(-2), B_i1(T)) / sqrt(z)
32+
end
33+
if x < zero(x)
34+
z = -z
35+
end
36+
return z
37+
end
38+
function besseli1x(x::Float64)
39+
T = Float64
40+
z = abs(x)
41+
if z <= 8
42+
y = muladd(x, T(.5), T(-2))
43+
z = chbevl(y, A_i1(T)) * z
44+
else
45+
z = chbevl(T(32) / x - T(-2), B_i1(T)) / sqrt(z)
46+
end
47+
if x < zero(x)
48+
z = -z
49+
end
50+
return z
51+
end

‎src/besselj.jl

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
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
9+
x = abs(x)
10+
iszero(x) && return one(x)
11+
isinf(x) && return zero(x)
12+
13+
if x <= 5
14+
z = x * x
15+
if x < 1.0e-5
16+
return 1.0 - z / 4.0
17+
end
18+
19+
DR1 = 5.78318596294678452118E0
20+
DR2 = 3.04712623436620863991E1
21+
22+
p = (z - DR1) * (z - DR2)
23+
p = p * evalpoly(z, RP_j0(T)) / evalpoly(z, RQ_j0(T))
24+
return p
25+
else
26+
w = 5.0 / x
27+
q = 25.0 / (x * x)
28+
29+
p = evalpoly(q, PP_j0(T)) / evalpoly(q, PQ_j0(T))
30+
q = evalpoly(q, QP_j0(T)) / evalpoly(q, QQ_j0(T))
31+
xn = x - PIO4(T)
32+
p = p * cos(xn) - w * q * sin(xn)
33+
return p * SQ2OPI(T) / sqrt(x)
34+
end
35+
end
36+
function besselj0(x::Float32)
37+
T = Float32
38+
x = abs(x)
39+
iszero(x) && return zero(x)
40+
isinf(x) && return zero(x)
41+
42+
if x <= 2.0f0
43+
z = x * x
44+
if x < 1.0f-3
45+
return 1.0f0 - 0.25f0 * z
46+
end
47+
DR1 = 5.78318596294678452118f0
48+
p = (z - DR1) * evalpoly(z, JP_j0(T))
49+
return p
50+
else
51+
q = inv(x)
52+
w = sqrt(q)
53+
p = w * evalpoly(q, MO_j0(T))
54+
w = q * q
55+
xn = q * evalpoly(w, PH_j0(T)) - PIO4(Float32)
56+
p = p * cos(xn + x)
57+
return p
58+
end
59+
end
60+
61+
function besselj1(x::Float64)
62+
T = Float64
63+
x = abs(x)
64+
iszero(x) && return zero(x)
65+
isinf(x) && return zero(x)
66+
67+
if x <= 5.0
68+
z = x * x
69+
w = evalpoly(z, RP_j1(T)) / evalpoly(z, RQ_j1(T))
70+
w = w * x * (z - 1.46819706421238932572e1) * (z - 4.92184563216946036703e1)
71+
return w
72+
else
73+
w = 5.0 / x
74+
z = w * w
75+
p = evalpoly(z, PP_j1(T)) / evalpoly(z, PQ_j1(T))
76+
q = evalpoly(z, QP_j1(T)) / evalpoly(z, QQ_j1(T))
77+
xn = x - THPIO4(T)
78+
p = p * cos(xn) - w * q * sin(xn)
79+
return p * SQ2OPI(T) / sqrt(x)
80+
end
81+
end
82+
83+
function besselj1(x::Float32)
84+
x = abs(x)
85+
iszero(x) && return zero(x)
86+
isinf(x) && return zero(x)
87+
88+
if x <= 2.0f0
89+
z = x * x
90+
Z1 = 1.46819706421238932572f1
91+
p = (z - Z1) * x * evalpoly(z, JP32)
92+
return p
93+
else
94+
q = inv(x)
95+
w = sqrt(q)
96+
p = w * evalpoly(q, MO132)
97+
w = q * q
98+
xn = q * evalpoly(w, PH132) - THPIO4(Float32)
99+
p = p * cos(xn + x)
100+
return p
101+
end
102+
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/besselk.jl

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
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/k0.c
5+
https://github.com/jeremybarnes/cephes/blob/master/bessel/k1.c
6+
=#
7+
function besselk0(x::T) where T <: Union{Float32, Float64}
8+
if x <= zero(x)
9+
return throw(DomainError(x, "NaN result for non-NaN input."))
10+
end
11+
12+
if x <= 2
13+
y = muladd(x, x, T(-2))
14+
y = chbevl(y, A_k0(T)) - log(T(.5) * x) * besseli0(x)
15+
return y
16+
else
17+
z = T(8) / x - T(2)
18+
return exp(-x) * chbevl(z, B_k0(T)) / sqrt(x)
19+
end
20+
end
21+
22+
function besselk0x(x::T) where T <: Union{Float32, Float64}
23+
if x <= zero(x)
24+
return throw(DomainError(x, "NaN result for non-NaN input."))
25+
end
26+
if x <= 2
27+
y = muladd(x, x, T(-2))
28+
y = chbevl(y, A_k0(T)) - log(T(.5) * x) * besseli0(x)
29+
return y * exp(x)
30+
else
31+
z = T(8) / x - T(2)
32+
return chbevl(z, B_k0(T)) / sqrt(x)
33+
end
34+
end
35+
function besselk1(x::T) where T <: Union{Float32, Float64}
36+
z = T(.5) * x
37+
if x <= zero(x)
38+
return throw(DomainError(x, "NaN result for non-NaN input."))
39+
end
40+
if x <= 2
41+
y = muladd(x, x, T(-2))
42+
y = log(z) * besseli1(x) + chbevl(y, A_k1(T)) / x
43+
return y
44+
else
45+
return exp(-x) * chbevl(T(8) / x - T(2), B_k1(T)) / sqrt(x)
46+
end
47+
end
48+
function besselk1x(x::T) where T <: Union{Float32, Float64}
49+
z = T(.5) * x
50+
if x <= zero(x)
51+
return throw(DomainError(x, "NaN result for non-NaN input."))
52+
end
53+
if x <= 2
54+
y = muladd(x, x, T(-2))
55+
y = log(z) * besseli1(x) + chbevl(y, A_k1(T)) / x
56+
return y * exp(x)
57+
else
58+
return chbevl(T(8) / x - T(2), B_k1(T)) / sqrt(x)
59+
end
60+
end

0 commit comments

Comments
 (0)