JuliaMath
return sequence of bessel functions
Moved from: https://github.com/JuliaLang/julia/issues/10395
I'm very new to julia, so not sure I managed to find the actual source code being used for those bessel(jhy) functions. There's a related thread suggesting that the chosen algorithm for integer arguments is different.
Nonetheless, for the general case, the Amos implementation returns a sequence of bessel functions, nu+I-1, I=1,...,N, and it would be very nice if the julia wrapper could return the whole sequence. Currently, I believe julia implements vectorisation over nu as a for loop and just takes the one value from the Amos sequence (or openlibm's). Since it is not uncommon to need a sequence of values (Bessel functions tend to be used in recurrence relations), it is a bit wasteful to iterate over nu when it was already done at the Fortran level.
Related question on StackOverflow
The actual source code is in https://github.com/JuliaLang/openspecfun ... the functions indeed have the option of returning a sequence of Bessel functions (see the N argument of zbesj, for example. It seems like it would be straightforward to add support for this for the case where the nu argument is a range with step(nu) == 1.
This change can be implemented by, for example, changing the following lines in bessel.jl:
ccall((:zbesh_,openspecfun), Cvoid,
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int},
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
real(z), imag(z), nu, kode, k, 1,
ai1, ai2, ae1, ae2)
to
ccall((:zbesh_,openspecfun), Cvoid,
(Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int},
Ptr{Float64}, Ptr{Float64}, Ref{Int32}, Ref{Int32}),
real(z), imag(z), nu, kode, k, N,
ai1, ai2, ae1, ae2)
where ai1,ai2 are arrays instead of references. Limited testing shows a drastic cut in time and allocated memory:
julia> @btime besselh.(0:5, 1, 1.1+0.1im)
5.658 μs (61 allocations: 1.11 KiB)
6-element Array{Complex{Float64},1}:
0.651368+0.11899im
0.393305-0.665955im
-0.0513027-1.38437im
-1.03222-4.31003im
-7.65254-21.4245im
-68.2153-145.21im
julia> @btime besselhseq(0, 5, 1, 1.1+0.1im)
1.460 μs (11 allocations: 560 bytes)
6-element Array{Complex{Float64},1}:
0.651368+0.11899im
0.393305-0.665955im
-0.0513027-1.38437im
-1.03222-4.31003im
-7.65254-21.4245im
-68.2153-145.21im
However, I'm unsure what would be the best way to wrap these functions. If using the existing signatures, one would have to handle the case where some Bessel orders are negative, since zbesh only works for nu >= 0.
