Skip to content

Commit 4a07d0b

Browse files
heltonmcheltonmc
committed
add docs update coefficients
1 parent 868d87f commit 4a07d0b

File tree

2 files changed

+100
-52
lines changed

2 files changed

+100
-52
lines changed

‎src/besseli.jl

Lines changed: 60 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,49 @@
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
1+
# Modified Bessel functions of the first kind of order zero and one
2+
# besseli0, besseli1
3+
#
4+
# Scaled modified Bessel functions of the first kind of order zero and one
5+
# besseli0x, besselix
6+
#
7+
# Calculation of besseli0 is done in two branches using polynomial approximations [1]
8+
#
9+
# Branch 1: x < 7.75
10+
# besseli0 = [x/2]^2 P16([x/2]^2)
11+
# Branch 2: x >=
12+
# sqrt(x) * exp(-x) * besseli0(x) = P22(1/x)
13+
# where P16 and P22 are a 16 and 22 degree polynomial respectively.
14+
#
15+
# Remez.jl is then used to approximate the polynomial coefficients of
16+
# P22(y) = sqrt(1/y) * exp(-inv(y)) * besseli0(inv(y))
17+
# N,D,E,X = ratfn_minimax(g, (1/1e6, 1/7.75), 21, 0)
18+
#
19+
# A third branch is used for scaled functions for large values
20+
#
21+
#
22+
# Calculation of besseli1 is done in two branches using polynomial approximations [2]
23+
#
24+
# Branch 1: x < 7.75
25+
# besseli1 = x / 2 * (1 + 1/2 * (x/2)^2 + (x/2)^4 * P13([x/2]^2)
26+
# Branch 2: x >=
27+
# sqrt(x) * exp(-x) * besseli0(x) = P22(1/x)
28+
# where P13 and P22 are a 16 and 22 degree polynomial respectively.
29+
#
30+
# Remez.jl is then used to approximate the polynomial coefficients of
31+
# P13(y) = (besseli1(2 * sqrt(y)) / sqrt(y) - 1 - y/2) / y^2
32+
# N,D,E,X = ratfn_minimax(g, (1/1e6, 1/7.75), 21, 0)
33+
#
34+
# A third branch is used for scaled functions for large values
35+
#
36+
# Horner's scheme is then used to evaluate all polynomials.
37+
#
38+
# [1] "Rational Approximations for the Modified Bessel Function of the First Kind
39+
# - I0(x) for Computations with Double Precision" by Pavel Holoborodko
40+
# [2] "Rational Approximations for the Modified Bessel Function of the First Kind
41+
# - I1(x) for Computations with Double Precision" by Pavel Holoborodko
42+
"""
43+
besseli0(x::T) where T <: Union{Float32, Float64}
544
6-
#=
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
13-
=#
45+
Modified Bessel function of the first kind of order zero, ``I_0(x)``.
46+
"""
1447
function besseli0(x::T) where T <: Union{Float32, Float64}
1548
x = abs(x)
1649
if x < 7.75
@@ -22,6 +55,11 @@ function besseli0(x::T) where T <: Union{Float32, Float64}
2255
return a * s
2356
end
2457
end
58+
"""
59+
besseli0x(x::T) where T <: Union{Float32, Float64}
60+
61+
Scaled modified Bessel function of the first kind of order zero, ``I_0(x)*e^{-x}``.
62+
"""
2563
function besseli0x(x::T) where T <: Union{Float32, Float64}
2664
T == Float32 ? branch = 50 : branch = 500
2765
x = abs(x)
@@ -34,6 +72,11 @@ function besseli0x(x::T) where T <: Union{Float32, Float64}
3472
return evalpoly(inv(x), besseli0_large_coefs(T)) / sqrt(x)
3573
end
3674
end
75+
"""
76+
besseli1(x::T) where T <: Union{Float32, Float64}
77+
78+
Modified Bessel function of the first kind of order one, ``I_1(x)``.
79+
"""
3780
function besseli1(x::T) where T <: Union{Float32, Float64}
3881
z = abs(x)
3982
if z < 7.75
@@ -50,6 +93,11 @@ function besseli1(x::T) where T <: Union{Float32, Float64}
5093
end
5194
return z
5295
end
96+
"""
97+
besseli1x(x::T) where T <: Union{Float32, Float64}
98+
99+
Scaled modified Bessel function of the first kind of order one, ``I_1(x)*e^{-x}``.
100+
"""
53101
function besseli1x(x::T) where T <: Union{Float32, Float64}
54102
T == Float32 ? branch = 50 : branch = 500
55103
z = abs(x)

‎src/constants.jl

Lines changed: 40 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -105,70 +105,70 @@ const QQ_j0(::Type{Float64}) = (
105105
6.43178256118178023184E1, 1.00000000000000000000E0
106106
)
107107
const besseli0_small_coefs(::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
108+
1.0000000416172086f0, 0.2499995567758189f0, 0.027778554423089f0,
109+
0.0017355879591086688f0, 6.96204574858338f-5, 1.8959194799355527f-6,
110+
4.298760429567388f-8, 3.8884953805853935f-10, 1.483799986601676f-11
111111
)
112112
const besseli0_small_coefs(::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
113+
0.9999999999999998, 0.2500000000000052, 0.027777777777755364,
114+
0.001736111111149161, 6.94444444107536e-5, 1.9290123635366806e-6,
115+
3.9367592765038015e-8, 6.151201574092085e-10, 7.593827956729909e-12,
116+
7.596677643342155e-14, 6.255282299620455e-16, 4.470993793303175e-18,
117+
2.1859737023077178e-20, 2.0941557335286373e-22
118118
)
119119
const besseli0_med_coefs(::Type{Float32}) = (
120-
3.98942651588301770f-1, 4.98327234176892844f-2, 2.91866904423115499f-2,
121-
1.35614940793742178f-2, 1.31409251787866793f-1
120+
0.3989423094522082f0, 0.049857687814784155f0, 0.028606362029333827f0,
121+
0.018926572151825607f0, 0.11423960556773063f0
122122
)
123123
const besseli0_med_coefs(::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
124+
0.3989422804014326, 0.04986778505064754, 0.028050628512954097,
125+
0.02921968830978531, 0.04466889626137549, 0.10220642174207666,
126+
-0.9937439085650689, 91.25330271974727, -4901.408890977662,
127+
199209.2752981982, -6.181516298413396e6, 1.4830278710991925e8,
128+
-2.7695254643719645e9, 4.0351394830842026e10, -4.5768930327229974e11,
129+
4.0134844243070063e12, -2.6862476523182016e13, 1.3437999451218112e14,
130+
-4.856333741437621e14, 1.1962791200680235e15, -1.796269414464399e15,
131+
1.239942074380968e15
132132
)
133133
const besseli0_large_coefs(::Type{Float32}) = (
134-
3.98942391532752700f-1, 4.98455950638200020f-2, 2.94835666900682535f-2
134+
0.398942288181429f0, 0.04986085645388595f0, 0.028961207830173957f0
135135
)
136136
const besseli0_large_coefs(::Type{Float64}) = (
137-
3.98942280401432905e-01, 4.98677850491434560e-02, 2.80506308916506102e-02,
138-
2.92179096853915176e-02, 4.53371208762579442e-02
137+
0.3989422804014327, 0.04986778505003549, 0.028050629664438713,
138+
0.029218603406797997, 0.045199090067282434
139139
)
140140
const besseli1_small_coefs(::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
141+
0.08333336088599108f0, 0.0069442679189687774f0, 0.00034740666242346436f0,
142+
1.1502079282185252f-5, 2.8884784673820464f-7, 3.686308746029452f-9,
143+
1.2271327751884605f-10
143144
)
144145
const besseli1_med_coefs(::Type{Float32}) = (
145146
3.98942115977513013f-1, -1.49581264836620262f-1, -4.76475741878486795f-2,
146147
-2.65157315524784407f-2, -1.47148600683672014f-1
147148
)
148149
const besseli1_small_coefs(::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
150+
0.08333333333333334, 0.006944444444444374, 0.00034722222222248526,
151+
1.1574074073690356e-5, 2.7557319253050506e-7, 4.920949730519126e-9,
152+
6.834656365321179e-11, 7.593985414952446e-13, 6.904652315442046e-15,
153+
5.2213850252454655e-17, 3.405120412140281e-19, 1.6398527256182257e-21,
154+
1.3161876924566675e-23
154155
)
155156
const besseli1_med_coefs(::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
157+
0.39894228040143276, -0.149603355151029, -0.04675104787903509,
158+
-0.04090746353279043, -0.05744911840910781,-0.12283724006390319,
159+
1.0023632527650936, -94.90954045770921, 5084.06899084327,
160+
-206253.5613716743, 6.387439538535799e6, -1.529244018923123e8,
161+
2.849523551208316e9, -4.141884344471782e10, 4.6860149658304974e11,
162+
-4.097852944580042e12, 2.7345051110005453e13, -1.3634833112030364e14,
163+
4.909983186948892e14, -1.2048200837913132e15, 1.8014682382937435e15,
164+
-1.2377987428989558e15
164165
)
165166
const besseli1_large_coefs(::Type{Float64}) = (
166-
3.989422804014314820e-01, -1.496033551467584157e-01, -4.675105322571775911e-02,
167-
-4.090421597376992892e-02, -5.843630344778927582e-02
167+
0.39894228040143265, -0.14960335515036183, -0.0467510491854453,
168+
-0.04090618769936314, -0.05808395099511361
168169
)
169170
const besseli1_large_coefs(::Type{Float32}) = (
170-
3.989422804014314820f-01, -1.496033551467584157f-01, -4.675105322571775911f-02,
171-
-4.090421597376992892f-02, -5.843630344778927582f-02
171+
0.3989422695655856f0, -0.1495936968314777f0, -0.04802187089304704f0
172172
)
173173
const A_k0(::Type{Float64}) = (
174174
1.37446543561352307156E-16, 4.25981614279661018399E-14, 1.03496952576338420167E-11,

0 commit comments

Comments
 (0)