Return sequence of Bessel functions

[heltonmc]

Here's a first draft to close #52. I still need to think about this a bit more but these will return right answers....

I think things to discuss:

1. Supporting negative indexes (I've returned an error here as the Amos routine only supports positive orders)
2. Handling overflow in forward recurrence and starting underflow in downward recurrence.

I think for 2 we may want to have a check that if the start value is zero we just error here instead of iterating through. That could create a very silent performance hit but want to think about that a little bit....

[codecov]

Codecov Report

Base: 94.17% // Head: 94.43% // Increases project coverage by +0.25% 🎉

Coverage data is based on head (4a59336) compared to base (cc97c8c).
Patch coverage: 97.18% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #53      +/-   ##
==========================================
+ Coverage   94.17%   94.43%   +0.25%     
==========================================
  Files          18       18              
  Lines        1820     1939     +119     
==========================================
+ Hits         1714     1831     +117     
- Misses        106      108       +2     

Impacted Files	Coverage Δ
src/hankel.jl	93.75% <92.30%> (+2.18%)	⬆️
src/besselk.jl	97.26% <96.00%> (-0.21%)	⬇️
src/besseli.jl	98.08% <96.15%> (-0.43%)	⬇️
src/besselj.jl	92.52% <96.42%> (+0.58%)	⬆️
src/bessely.jl	91.08% <100.00%> (+0.22%)	⬆️
src/modifiedsphericalbessel.jl	97.91% <100.00%> (ø)
src/recurrence.jl	100.00% <100.00%> (ø)
src/sphericalbessel.jl	100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[heltonmc]

I've added this functionality for all the functions now besselj, bessely, besselk, besseli and the hankel functions.

I'm still unsure the best approach if the starting values of the recurrence underflow as I can see advantages either way. Presumably you are using this method because it is efficient, so if we throw an error it will alert the user that this approach may not be efficient and they can adjust their range of nu or just use broadcasting for large orders. Alternatively, we could loop through until it is finite but that to my knowledge would require evaluation of each value to check if nonzero which would be quite inefficient.

An example would be something like besselj(0:1000, 5.0). This function returns 0 in Float64 precision at besselj(210, 5.0) so starting downward recurrence at 1000 would require almost 800 individual evaluations of besselj to reach a point where we can start recurrence.

@lucifer1004 @krcools Do you have a preference?

EDIT: I'm going to go with returning the array instead of returning an error even if it requires evaluation until finite values. I also need to take care not to start with subnormal numbers....

[krcools]

Hi, sorry for not replying sooner. I feel a bit out of my depth here. It actually never occurred to me this could happen. Any idea how AMOS/matlab deals with this?

[heltonmc]

No worries I appreciate the thoughts and feedback as I keep going back and forth on this! Matlab will return the correct results (zeros if underflow in that region and values for lower orders). It does look like there is some performance hit but it is able to account for this fairly well.

The plan for me forward is to have a good underflow check that is much faster than computing the besselj function and will work for most cases where underflow will occur. This will reduce the need to compute besselj many times while always returning the correct result. I'm working on finding as simple of a check for this now so hopefully that is the last piece.

Thank you!

[heltonmc]

Ok I think everything looks good now. This will properly handle starting at underflow.

julia> besselj(0:500, 20.0)
501-element Vector{Float64}:
  0.16702466434057608
  0.06683312417584697
 -0.1603413519229914
 -0.09890139456044525
  0.13067093355485782
  0.15116976798238838
 -0.05508604956366364
 -0.18422139772058657
 -0.07386892884074697
  0.125126254647989
  0.18648255802393707
  0.061356303375948065
 -0.11899062431039419
  ⋮
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0
  0.0

So we get appropriate results even when a lot of values are zero. The performance cost here is fairly small...

julia> @btime besselj(0:500, 20.0)
  7.859 μs (2 allocations: 6.50 KiB)

So for 500 Bessel functions it takes just several microseconds to perform even when the majority of the values underflow. Of course if there is no overflow it will be faster...

julia> @btime besselj(0:500, 400.0)
  1.820 μs (2 allocations: 8.12 KiB)

So those are worst case scenarios (2 allocations have to check for underflow). The best case we get when x > nu.

julia> @btime besselj(0:500, 600.0)
  1.332 μs (1 allocation: 4.06 KiB)

So in that case only 1 allocation is needed and very fast.

Anyway I think this is ready to go now.