Skip to content

Commit 2c874b8

Browse files
heltonmcheltonmc
authored
Merge pull request #11 from heltonmc/benchmark
benchmark suite
2 parents 47bf072 + a39f417 commit 2c874b8

File tree

2 files changed

+73
-149
lines changed

2 files changed

+73
-149
lines changed

‎README.md

Lines changed: 47 additions & 149 deletions
Original file line numberDiff line numberDiff line change
@@ -2,173 +2,71 @@
22
[![Build Status](https://github.com/heltonmc/Bessels.jl/actions/workflows/CI.yml/badge.svg?branch=master)](https://github.com/heltonmc/Bessels.jl/actions/workflows/CI.yml?query=branch%3Amaster)
33
[![Coverage](https://codecov.io/gh/heltonmc/Bessels.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/heltonmc/Bessels.jl)
44

5+
A Julia implementation of Bessel's functions and modified Bessel's functions of the first and second kind.
56

6-
Implementations of Bessel's functions `besselj0`, `besselj1`, `bessely0`, `bessely1`, `besseli0`, `besseli1`, `besselk0`, `besselk1` in Julia.
7-
Most implementations are ported to Julia from the [Cephes](https://www.netlib.org/cephes/) math library originally written by Stephen L. Moshier in the C programming language
8-
which are (partly) used in [SciPy](https://docs.scipy.org/doc/scipy/reference/special.html#faster-versions-of-common-bessel-functions) and [GCC libquadmath](https://gcc.gnu.org/onlinedocs/libquadmath/).
7+
The library currently supports Bessel's function of the first and second kind for orders 0 and 1 and for any integer order for modified Bessel's functions for real arguments.
98

109
# Quick start
1110

12-
Only implemented for real arguments so far.
13-
1411
```julia
15-
using Bessels
12+
# add the package
13+
pkg> add https://github.com/heltonmc/Bessels.jl.git
14+
15+
julia> using Bessels
1616

1717
# Bessel function of the first kind of order 0
1818
julia> besselj0(1.0)
1919
0.7651976865579665
20+
2021
julia> besselj0(1.0f0)
2122
0.7651977f0
22-
23-
# Bessel function of the first kind of order 1
24-
julia> besselj1(1.0)
25-
0.44005058574493355
26-
julia> besselj1(1.0f0)
27-
0.4400506f0
28-
29-
# Bessel function of the second kind of order 0
30-
julia> bessely0(1.0)
31-
0.08825696421567697
32-
julia> bessely0(1.0f0)
33-
0.08825697f0
34-
35-
# Bessel function of the second kind of order 1
36-
julia> bessely1(1.0)
37-
-0.7812128213002888
38-
julia> bessely1(1.0f0)
39-
-0.7812128f0
40-
41-
# Modified Bessel function of the first kind of order 0
42-
julia> besseli0(1.0)
43-
1.2660658777520082
44-
julia> besseli0(1.0f0)
45-
1.266066f0
46-
47-
# Scaled modified Bessel function of the first kind of order 0
48-
julia> besseli0x(1.0)
49-
0.46575960759364043
50-
julia> besseli0x(1.0f0)
51-
0.46575963f0
52-
53-
# Scaled modified Bessel function of the first kind of order 1
54-
julia> besseli1x(1.0)
55-
0.2079104153497085
56-
julia> besseli1x(1.0f0)
57-
0.20791042f0
58-
59-
# Modified Bessel function of the second kind of order 0
60-
julia> besselk0(1.0)
61-
0.42102443824070823
62-
julia> besselk0(1.0f0)
63-
0.42102447f0
64-
65-
# Scaled modified Bessel function of the second kind of order 0
66-
julia> besselk0x(1.0)
67-
1.1444630798068947
68-
julia> besselk0x(1.0f0)
69-
1.1444632f0
70-
71-
# Modified Bessel function of the second kind of order 1
72-
julia> besselk1(1.0)
73-
0.6019072301972346
74-
julia> besselk1(1.0f0)
75-
0.6019073f0
76-
77-
# Scaled modified Bessel function of the second kind of order 1
78-
julia> besselk1x(1.0)
79-
1.636153486263258
80-
julia> besselk1x(1.0f0)
81-
1.6361537f0
8223
```
8324

25+
# Supported functions
26+
27+
- `besselj0(x)`
28+
- `besselj1(x)`
29+
- `bessely0(x)`
30+
- `bessely1(x)`
31+
- `besseli0(x)`
32+
- `besseli1(x)`
33+
- `besselk0(x)`
34+
- `besselk1(x)`
35+
- `besseli(nu, x)`
36+
- `besselk(nu, x)`
37+
38+
Exponentially scaled versions of Modified Bessel's function can also be called with `besselix(nu, x)` and `besselkx(nu, x)`.
39+
8440
# Benchmarks
8541

86-
```julia
87-
# Bessels.jl
88-
julia> @benchmark besselj0(x) setup=(x=100.0*rand())
89-
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
90-
Range (min … max): 7.465 ns … 28.946 ns ┊ GC (min … max): 0.00% … 0.00%
91-
Time (median): 25.860 ns ┊ GC (median): 0.00%
92-
Time (mean ± σ): 24.915 ns ± 4.090 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
93-
Memory estimate: 0 bytes, allocs estimate: 0.
94-
95-
# SpecialFunctions.jl
96-
julia> @benchmark besselj0(x) setup=(x=100.0*rand())
97-
BenchmarkTools.Trial: 10000 samples with 999 evaluations.
98-
Range (min … max): 9.050 ns … 79.412 ns ┊ GC (min … max): 0.00% … 0.00%
99-
Time (median): 51.677 ns ┊ GC (median): 0.00%
100-
Time (mean ± σ): 51.029 ns ± 5.653 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
101-
Memory estimate: 0 bytes, allocs estimate: 0.
102-
```
103-
104-
```julia
105-
# Bessels.jl
106-
julia> @benchmark besselk1(x) setup=(x=100.0*rand())
107-
BenchmarkTools.Trial: 10000 samples with 988 evaluations.
108-
Range (min … max): 47.781 ns … 97.672 ns ┊ GC (min … max): 0.00% … 0.00%
109-
Time (median): 47.950 ns ┊ GC (median): 0.00%
110-
Time (mean ± σ): 48.495 ns ± 2.656 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
111-
Memory estimate: 0 bytes, allocs estimate: 0.
112-
113-
# SpecialFunctions.jl
114-
julia> @benchmark besselk(1, x) setup=(x=10.0*rand())
115-
BenchmarkTools.Trial: 10000 samples with 681 evaluations.
116-
Range (min … max): 184.962 ns … 2.122 μs ┊ GC (min … max): 0.00% … 81.41%
117-
Time (median): 377.081 ns ┊ GC (median): 0.00%
118-
Time (mean ± σ): 388.243 ns ± 96.971 ns ┊ GC (mean ± σ): 0.04% ± 0.81%
119-
Memory estimate: 16 bytes, allocs estimate: 1.
120-
```
121-
122-
Comparing the relative speed (`SpecialFunctions.jl / Bessels.jl`) for a vector of values between 0 and 100.
123-
124-
## Float32
125-
126-
| function | Relative speed |
127-
| ------------- | ------------- |
128-
| besselj0 | 1.7x |
129-
| besselj1 | 1.7x |
130-
| bessely0 | 1.6x |
131-
| bessely1 | 1.6x |
132-
| besseli0 | 24x |
133-
| besseli1 | 22x |
134-
| besselk0 | 13x |
135-
| besselk1 | 15x |
136-
137-
* it looks like SpecialFunctions.jl doesn't always preserve the correct input type so some of the calculations may be done in Float64. This might skew the benchmarks for `Bessels.jl` as it should have all calculations done in the lower precision.
42+
These results compare the median value from BenchmarkTools obtained on one machine for arguments between 0 and 100.
13843

139-
```julia
140-
# SpecialFunctions.jl
141-
julia> @btime besselk(1, x) setup=(x=Float32(10.0*rand()))
142-
179.272 ns (1 allocation: 16 bytes)
143-
0.021948646168061196 # notice incorrect Float64 return type
44+
We compare the relative sped to implementations provided by [SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl).
14445

145-
# Bessels.jl
146-
julia> @btime besselk1(x) setup=(x=Float32(10.0*rand()))
147-
22.967 ns (0 allocations: 0 bytes)
148-
0.0027777348f0 # notice correct Float32 return type
149-
```
46+
| function | `Float32` | `Float64`
47+
| ------------- | ------------- | ------------- |
48+
| besselj0 | 1.7x | 1.8x
49+
| besselj1 | 1.7x | 1.9x
50+
| bessely0 | 1.9x | 1.8x
51+
| bessely1 | 1.7x | 1.7x
52+
| besseli0 | 26x | 13.2x
53+
| besseli1 | 22x | 13.9x
54+
| besseli(20, x) | 5.4x | 2.1x |
55+
| besseli(120, x) | 6x | 4.5x |
56+
| besselk0 | 16x | 12.7x
57+
| besselk1 | 15x | 14.3x
58+
| besselk(20, x) | 5.4x | 5.7x |
59+
| besselk(120, x) | 2.9x | 3.4x |
15060

151-
## Float64
61+
* SpecialFunctions.jl doesn't always preserve the correct input type so some of the calculations may be done in Float64. This might skew the benchmarks for `Bessels.jl`.
15262

153-
| function | Relative speed |
154-
| ------------- | ------------- |
155-
| besselj0 | 1.7x |
156-
| besselj1 | 1.9x |
157-
| bessely0 | 1.5x |
158-
| bessely1 | 1.5x |
159-
| besseli0 | 9x |
160-
| besseli1 | 7x |
161-
| besselk0 | 5x |
162-
| besselk1 | 5x |
63+
# Current Development Plans
16364

164-
```julia
165-
# SpecialFunctions.jl
166-
julia> @btime besselk(1, x) setup=(x=10.0*rand())
167-
184.726 ns (1 allocation: 16 bytes) # notice the small difference in Float32 and Float64 implementations
168-
0.0007517428778913419
65+
- More robust implementations of `besselj(nu, x)` and `bessely(nu, x)`
66+
- Support non-integer orders
67+
- Support complex arguments
16968

170-
# Bessels.jl
171-
julia> @btime besselk1(x) setup=(x=10.0*rand())
172-
47.824 ns (0 allocations: 0 bytes)
173-
0.0057366790518002045
174-
```
69+
# Contributing
70+
71+
Contributions are very welcome, as are feature requests and suggestions. Please open an [issue](https://github.com/heltonmc/Bessels.jl/issues/new) for discussion on newer implementations, share papers, new features, or if you encounter any problems.
72+

‎benchmark/benchmark.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
using BenchmarkTools
2+
using Bessels
3+
4+
suite = BenchmarkGroup()
5+
6+
suite["besseli"] = BenchmarkGroup()
7+
suite["besseli"]["besseli0"] = @benchmarkable besseli0(x) setup=(x = rand()*100)
8+
suite["besseli"]["besseli1"] = @benchmarkable besseli1(x) setup=(x = rand()*100)
9+
suite["besseli"]["besseli20"] = @benchmarkable besseli(20, x) setup=(x = rand()*100)
10+
suite["besseli"]["besseli120"] = @benchmarkable besseli(120, x) setup=(x = rand()*100 + 1.0)
11+
12+
suite["besselk"] = BenchmarkGroup()
13+
suite["besselk"]["besselk0"] = @benchmarkable besselk0(x) setup=(x = rand()*100)
14+
suite["besselk"]["besselk1"] = @benchmarkable besselk1(x) setup=(x = rand()*100)
15+
suite["besselk"]["besselk20"] = @benchmarkable besselk(20, x) setup=(x = rand()*100)
16+
suite["besselk"]["besselk120"] = @benchmarkable besselk(120, x) setup=(x = rand()*100 + 1.0)
17+
18+
suite["besselj"] = BenchmarkGroup()
19+
suite["besselj"]["besselj0"] = @benchmarkable besselj0(x) setup=(x = rand()*100)
20+
suite["besselj"]["besselj1"] = @benchmarkable besselj1(x) setup=(x = rand()*100)
21+
22+
suite["bessely"] = BenchmarkGroup()
23+
suite["bessely"]["bessely0"] = @benchmarkable besselj0(x) setup=(x = rand()*100)
24+
suite["bessely"]["bessely1"] = @benchmarkable besselj1(x) setup=(x = rand()*100)
25+
26+
#tune!(suite)

0 commit comments

Comments
 (0)