Skip to content

Commit a12f2a9

Browse files
heltonmcheltonmc
committed
add evalpoly routines of besseli
1 parent 547bfad commit a12f2a9

File tree

2 files changed

+129
-102
lines changed

2 files changed

+129
-102
lines changed

‎src/besseli.jl

Lines changed: 68 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,49 +1,86 @@
1+
# Modified Bessel functions of the first kind of order zero and one
2+
# besseli0, besseli1
3+
# Scaled modified Bessel functions of the first kind of order zero and one
4+
# besseli0x, besselix
5+
16
#=
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
7+
Approximate forms used here are given in
8+
"Rational Approximations for the Modified Bessel Function of the First Kind - I1(x) for Computations with Double Precision"
9+
by Pavel Holoborodko,
10+
http://www.advanpix.com/2015/11/12/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision/
11+
Actual coefficients used are from the boost math library.
12+
https://github.com/boostorg/math/tree/develop/include/boost/math/special_functions/detail
613
=#
14+
715
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))
16+
T == Float32 ? branch = 50 : branch = 500
17+
if x < 7.75
18+
a = x * x / 4
19+
return muladd(a, evalpoly(a, P1_i0(T)), 1)
20+
elseif x < branch
21+
return exp(x) * evalpoly(inv(x), P2_i0(T)) / sqrt(x)
1122
else
12-
return exp(x) * chbevl(T(32) / x - T(2), B_i0(T)) / sqrt(x)
23+
a = exp(x / 2)
24+
s = a * evalpoly(inv(x), P3_i0(T)) / sqrt(x)
25+
return a * s
1326
end
1427
end
1528
function besseli0x(x::T) where T <: Union{Float32, Float64}
16-
x = abs(x)
17-
if x <= 8
18-
y = muladd(x, T(.5), T(-2))
19-
return chbevl(y, A_i0(T))
29+
T == Float32 ? branch = 50 : branch = 500
30+
if x < 7.75
31+
a = x * x / 4
32+
return muladd(a, evalpoly(a, P1_i0(T)), 1) * exp(-x)
33+
elseif x < branch
34+
return evalpoly(inv(x), P2_i0(T)) / sqrt(x)
2035
else
21-
return chbevl(T(32) / x - T(2), B_i0(T)) / sqrt(x)
36+
return evalpoly(inv(x), P3_i0(T)) / sqrt(x)
2237
end
2338
end
24-
function besseli1(x::T) where T <: Union{Float32, Float64}
25-
z = abs(x)
26-
if x <= 8
27-
y = muladd(z, T(.5), T(-2))
28-
z = chbevl(y, A_i1(T)) * z * exp(z)
39+
function besseli1(x::Float32)
40+
T = Float32
41+
if x < 7.75
42+
a = x * x / 4
43+
inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
44+
return x * evalpoly(a, inner) / 2
2945
else
30-
z = exp(z) * chbevl(T(32) / z - T(2), B_i1(T)) / sqrt(z)
46+
a = exp(x / 2)
47+
s = a * evalpoly(inv(x), P2_i1(T)) / sqrt(x)
48+
return a * s
3149
end
32-
if x < zero(x)
33-
z = -z
50+
end
51+
function besseli1(x::Float64)
52+
T = Float64
53+
if x < 7.75
54+
a = x * x / 4
55+
inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
56+
return x * evalpoly(a, inner) / 2
57+
elseif x < 500
58+
return exp(x) * evalpoly(inv(x), P2_i1(T)) / sqrt(x)
59+
else
60+
a = exp(x / 2)
61+
s = a * evalpoly(inv(x), P3_i1(T)) / sqrt(x)
62+
return a * s
3463
end
35-
return z
3664
end
37-
function besseli1x(x::T) where T <: Union{Float32, Float64}
38-
z = abs(x)
39-
if z <= 8
40-
y = muladd(z, T(.5), T(-2))
41-
z = chbevl(y, A_i1(T)) * z
65+
function besseli1x(x::Float32)
66+
T = Float32
67+
if x < 7.75
68+
a = x * x / 4
69+
inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
70+
return x * evalpoly(a, inner) / 2 * exp(-x)
4271
else
43-
z = chbevl(T(32) / z - T(2), B_i1(T)) / sqrt(z)
72+
return evalpoly(inv(x), P2_i1(T)) / sqrt(x)
4473
end
45-
if x < zero(x)
46-
z = -z
74+
end
75+
function besseli1x(x::Float64)
76+
T = Float64
77+
if x < 7.75
78+
a = x * x / 4
79+
inner = (one(T), T(0.5), evalpoly(a, P1_i1(T)))
80+
return x * evalpoly(a, inner) / 2 * exp(-x)
81+
elseif x < 500
82+
return evalpoly(inv(x), P2_i1(T)) / sqrt(x)
83+
else
84+
return evalpoly(inv(x), P3_i1(T)) / sqrt(x)
4785
end
48-
return z
4986
end

‎src/constants.jl

Lines changed: 61 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -104,51 +104,67 @@ const QQ_j0(::Type{Float64}) = (
104104
7.24046774195652478189E3, 3.88240183605401609683E3, 8.56430025976980587198E2,
105105
6.43178256118178023184E1, 1.00000000000000000000E0
106106
)
107-
const A_i0(::Type{Float64}) = (
108-
-4.41534164647933937950E-18, 3.33079451882223809783E-17, -2.43127984654795469359E-16,
109-
1.71539128555513303061E-15, -1.16853328779934516808E-14, 7.67618549860493561688E-14,
110-
-4.85644678311192946090E-13, 2.95505266312963983461E-12, -1.72682629144155570723E-11,
111-
9.67580903537323691224E-11, -5.18979560163526290666E-10, 2.65982372468238665035E-9,
112-
-1.30002500998624804212E-8, 6.04699502254191894932E-8, -2.67079385394061173391E-7,
113-
1.11738753912010371815E-6, -4.41673835845875056359E-6, 1.64484480707288970893E-5,
114-
-5.75419501008210370398E-5, 1.88502885095841655729E-4, -5.76375574538582365885E-4,
115-
1.63947561694133579842E-3, -4.32430999505057594430E-3, 1.05464603945949983183E-2,
116-
-2.37374148058994688156E-2, 4.93052842396707084878E-2, -9.49010970480476444210E-2,
117-
1.71620901522208775349E-1, -3.04682672343198398683E-1, 6.76795274409476084995E-1
118-
)
119-
const B_i0(::Type{Float64}) = (
120-
-7.23318048787475395456E-18, -4.83050448594418207126E-18, 4.46562142029675999901E-17,
121-
3.46122286769746109310E-17, -2.82762398051658348494E-16, -3.42548561967721913462E-16,
122-
1.77256013305652638360E-15, 3.81168066935262242075E-15, -9.55484669882830764870E-15,
123-
-4.15056934728722208663E-14, 1.54008621752140982691E-14, 3.85277838274214270114E-13,
124-
7.18012445138366623367E-13, -1.79417853150680611778E-12, -1.32158118404477131188E-11,
125-
-3.14991652796324136454E-11, 1.18891471078464383424E-11, 4.94060238822496958910E-10,
126-
3.39623202570838634515E-9, 2.26666899049817806459E-8, 2.04891858946906374183E-7,
127-
2.89137052083475648297E-6, 6.88975834691682398426E-5, 3.36911647825569408990E-3,
128-
8.04490411014108831608E-1
129-
)
130-
const A_i1(::Type{Float64}) = (
131-
2.77791411276104639959E-18, -2.11142121435816608115E-17, 1.55363195773620046921E-16,
132-
-1.10559694773538630805E-15, 7.60068429473540693410E-15, -5.04218550472791168711E-14,
133-
3.22379336594557470981E-13, -1.98397439776494371520E-12, 1.17361862988909016308E-11,
134-
-6.66348972350202774223E-11, 3.62559028155211703701E-10, -1.88724975172282928790E-9,
135-
9.38153738649577178388E-9, -4.44505912879632808065E-8, 2.00329475355213526229E-7,
136-
-8.56872026469545474066E-7, 3.47025130813767847674E-6, -1.32731636560394358279E-5,
137-
4.78156510755005422638E-5, -1.61760815825896745588E-4, 5.12285956168575772895E-4,
138-
-1.51357245063125314899E-3, 4.15642294431288815669E-3, -1.05640848946261981558E-2,
139-
2.47264490306265168283E-2, -5.29459812080949914269E-2, 1.02643658689847095384E-1,
140-
-1.76416518357834055153E-1, 2.52587186443633654823E-1
141-
)
142-
const B_i1(::Type{Float64}) = (
143-
7.51729631084210481353E-18, 4.41434832307170791151E-18, -4.65030536848935832153E-17,
144-
-3.20952592199342395980E-17, 2.96262899764595013876E-16, 3.30820231092092828324E-16,
145-
-1.88035477551078244854E-15, -3.81440307243700780478E-15, 1.04202769841288027642E-14,
146-
4.27244001671195135429E-14, -2.10154184277266431302E-14, -4.08355111109219731823E-13,
147-
-7.19855177624590851209E-13, 2.03562854414708950722E-12, 1.41258074366137813316E-11,
148-
3.25260358301548823856E-11, -1.89749581235054123450E-11, -5.58974346219658380687E-10,
149-
-3.83538038596423702205E-9, -2.63146884688951950684E-8, -2.51223623787020892529E-7,
150-
-3.88256480887769039346E-6, -1.10588938762623716291E-4, -9.76109749136146840777E-3,
151-
7.78576235018280120474E-1
107+
const P1_i0(::Type{Float32}) = (
108+
1.00000003928615375f0, 2.49999576572179639f-1, 2.77785268558399407f-2,
109+
1.73560257755821695f-3, 6.96166518788906424f-5, 1.89645733877137904f-6,
110+
4.29455004657565361f-8, 3.90565476357034480f-10, 1.48095934745267240f-11
111+
)
112+
const P1_i0(::Type{Float64}) = (
113+
1.00000000000000000e0, 2.49999999999999909e-1, 2.77777777777782257e-2,
114+
1.73611111111023792e-3, 6.94444444453352521e-5, 1.92901234513219920e-6,
115+
3.93675991102510739e-8, 6.15118672704439289e-10, 7.59407002058973446e-12,
116+
7.59389793369836367e-14, 6.27767773636292611e-16, 4.34709704153272287e-18,
117+
2.63417742690109154e-20, 1.13943037744822825e-22, 9.07926920085624812e-25
118+
)
119+
const P2_i0(::Type{Float32}) = (
120+
3.98942651588301770f-1, 4.98327234176892844f-2, 2.91866904423115499f-2,
121+
1.35614940793742178f-2, 1.31409251787866793f-1
122+
)
123+
const P2_i0(::Type{Float64}) = (
124+
3.98942280401425088e-01, 4.98677850604961985e-02, 2.80506233928312623e-02,
125+
2.92211225166047873e-02, 4.44207299493659561e-02, 1.30970574605856719e-01,
126+
-3.35052280231727022e+00, 2.33025711583514727e+02, -1.13366350697172355e+04,
127+
4.24057674317867331e+05, -1.23157028595698731e+07, 2.80231938155267516e+08,
128+
-5.01883999713777929e+09, 7.08029243015109113e+10, -7.84261082124811106e+11,
129+
6.76825737854096565e+12, -4.49034849696138065e+13, 2.24155239966958995e+14,
130+
-8.13426467865659318e+14, 2.02391097391687777e+15, -3.08675715295370878e+15,
131+
2.17587543863819074e+15
132+
)
133+
const P3_i0(::Type{Float32}) = (
134+
3.98942391532752700f-1, 4.98455950638200020f-2, 2.94835666900682535f-2
135+
)
136+
const P3_i0(::Type{Float64}) = (
137+
3.98942280401432905e-01, 4.98677850491434560e-02, 2.80506308916506102e-02,
138+
2.92179096853915176e-02, 4.53371208762579442e-02
139+
)
140+
const P1_i1(::Type{Float32}) = (
141+
8.333333221f-2, 6.944453712f-3, 3.472097211f-4, 1.158047174f-5,
142+
2.739745142f-7, 5.135884609f-9, 5.262251502f-11, 1.331933703f-12
143+
)
144+
const P2_i1(::Type{Float32}) = (
145+
3.98942115977513013f-1, -1.49581264836620262f-1, -4.76475741878486795f-2,
146+
-2.65157315524784407f-2, -1.47148600683672014f-1
147+
)
148+
const P1_i1(::Type{Float64}) = (
149+
8.333333333333333803e-02, 6.944444444444341983e-03, 3.472222222225921045e-04,
150+
1.157407407354987232e-05, 2.755731926254790268e-07, 4.920949692800671435e-09,
151+
6.834657311305621830e-11, 7.593969849687574339e-13, 6.904822652741917551e-15,
152+
5.220157095351373194e-17, 3.410720494727771276e-19, 1.625212890947171108e-21,
153+
1.332898928162290861e-23
154+
)
155+
const P2_i1(::Type{Float64}) = (
156+
3.989422804014406054e-01, -1.496033551613111533e-01, -4.675104253598537322e-02,
157+
-4.090895951581637791e-02, -5.719036414430205390e-02, -1.528189554374492735e-01,
158+
3.458284470977172076e+00, -2.426181371595021021e+02, 1.178785865993440669e+04,
159+
-4.404655582443487334e+05, 1.277677779341446497e+07, -2.903390398236656519e+08,
160+
5.192386898222206474e+09, -7.313784438967834057e+10, 8.087824484994859552e+11,
161+
-6.967602516005787001e+12, 4.614040809616582764e+13, -2.298849639457172489e+14,
162+
8.325554073334618015e+14, -2.067285045778906105e+15, 3.146401654361325073e+15,
163+
-2.213318202179221945e+15
164+
)
165+
const P3_i1(::Type{Float64}) = (
166+
3.989422804014314820e-01, -1.496033551467584157e-01, -4.675105322571775911e-02,
167+
-4.090421597376992892e-02, -5.843630344778927582e-02
152168
)
153169
const A_k0(::Type{Float64}) = (
154170
1.37446543561352307156E-16, 4.25981614279661018399E-14, 1.03496952576338420167E-11,
@@ -184,32 +200,6 @@ const B_k1(::Type{Float64}) = (
184200
1.95215518471351631108E-4, -2.85781685962277938680E-3, 1.03923736576817238437E-1,
185201
2.72062619048444266945E0
186202
)
187-
const A_i0(::Type{Float32}) = (
188-
-1.30002500998624804212f-8, 6.04699502254191894932f-8, -2.67079385394061173391f-7,
189-
1.11738753912010371815f-6, -4.41673835845875056359f-6, 1.64484480707288970893f-5,
190-
-5.75419501008210370398f-5, 1.88502885095841655729f-4, -5.76375574538582365885f-4,
191-
1.63947561694133579842f-3, -4.32430999505057594430f-3, 1.05464603945949983183f-2,
192-
-2.37374148058994688156f-2, 4.93052842396707084878f-2, -9.49010970480476444210f-2,
193-
1.71620901522208775349f-1, -3.04682672343198398683f-1, 6.76795274409476084995f-1
194-
)
195-
const B_i0(::Type{Float32}) = (
196-
3.39623202570838634515f-9, 2.26666899049817806459f-8, 2.04891858946906374183f-7,
197-
2.89137052083475648297f-6, 6.88975834691682398426f-5, 3.36911647825569408990f-3,
198-
8.04490411014108831608f-1
199-
)
200-
const A_i1(::Type{Float32}) = (
201-
9.38153738649577178388f-9, -4.44505912879632808065f-8, 2.00329475355213526229f-7,
202-
-8.56872026469545474066f-7, 3.47025130813767847674f-6, -1.32731636560394358279f-5,
203-
4.78156510755005422638f-5, -1.61760815825896745588f-4, 5.12285956168575772895f-4,
204-
-1.51357245063125314899f-3, 4.15642294431288815669f-3, -1.05640848946261981558f-2,
205-
2.47264490306265168283f-2, -5.29459812080949914269f-2, 1.02643658689847095384f-1,
206-
-1.76416518357834055153f-1, 2.52587186443633654823f-1
207-
)
208-
const B_i1(::Type{Float32}) = (
209-
-3.83538038596423702205f-9, -2.63146884688951950684f-8, -2.51223623787020892529f-7,
210-
-3.88256480887769039346f-6, -1.10588938762623716291f-4, -9.76109749136146840777f-3,
211-
7.78576235018280120474f-1
212-
)
213203
const JP_j0(::Type{Float32}) = (
214204
-1.729150680240724f-1, 1.332913422519003f-2, -3.969646342510940f-4,
215205
6.388945720783375f-6, -6.068350350393235f-8

0 commit comments

Comments
 (0)