Add specialised Bessel functions from SLATEC

[eprovst]

Currently openspecfun only includes the Z- versions provided by SLATEC which expect a double precision complex number as input.

However there are also D- versions specialised for double precision reals, and further specialised -0 and -1 versions hereof. Given Julia's multiple dispatch SpecialFunctions.jl could automatically select the fastest option, if these became available here. Is this feasible?

[eprovst]

For the use cases of DBESJ0, DBESY0, DBESJ1 and DBESY1 SpecialFunctions.jl already uses libm's, the same goes for other integer values for n in Jn and Yn.

I do not know how these methods compare to the SLATEC version, though.

This leaves DBESK0, DBESK1, DBESK, DBESI0, DBESI1, DBESI, DBESJ, DBESY as candidates for inclusion.

[ViralBShah]

I would have generally expected libm to be more robust than slatec. But I have no direct experience with these.

[pdubya]

I did a quick check against netlib/specfun equivalents of dbesk0, dbesk1 and the scaled versions. Seems to be about 10x faster than the generic order complex versions

function besselk0(x::Real)
    ccall((:besk0_, "libbesktest.so"), Cdouble, (Ref{Cdouble},), x)
end

function besselk0x(x::Real)
    ccall((:besek0_, "libbesktest.so"), Cdouble, (Ref{Cdouble},), x)
end

function besselk1(x::Real)
    ccall((:besk1_, "libbesktest.so"), Cdouble, (Ref{Cdouble},), x)
end

function besselk1x(x::Real)
    ccall((:besek1_, "libbesktest.so"), Cdouble, (Ref{Cdouble},), x)
end


using SpecialFunctions: besselk, besselkx

println(besselk(0, 0.5))
println(besselk0(0.5))
println(besselkx(0, 0.5))
println(besselk0x(0.5))
println(besselk(1, 0.5))
println(besselk1(0.5))
println(besselkx(1, 0.5))
println(besselk1x(0.5))


println("k0")
@time begin
    for j = 1:10000
        besselk0(0.5)
    end
end

@time begin
    for j = 1:10000
        besselk(0, 0.5)
    end
end

println("k0x")
@time begin
    for j = 1:10000
        besselk0x(0.5)
    end
end

@time begin
    for j = 1:10000
        besselkx(0, 0.5)
    end
end

println("k1")
@time begin
    for j = 1:10000
        besselk1(0.5)
    end
end

@time begin
    for j = 1:10000
        besselk(1, 0.5)
    end
end

println("k1x")
@time begin
    for j = 1:10000
        besselk1x(0.5)
    end
end

@time begin
    for j = 1:10000
        besselkx(1, 0.5)
    end
end

0.9244190712276656
0.9244190712276659
1.5241093857739092
1.5241093857739096
1.6564411200033007
1.656441120003301
2.7310097082117855
2.731009708211786
k0
  0.000115 seconds
  0.001876 seconds
k0x
  0.000154 seconds
  0.001865 seconds
k1
  0.000145 seconds
  0.001732 seconds
k1x
  0.000282 seconds
  0.001888 seconds

[ViralBShah]

@stevengj Would it be useful to pull in the double precision variants? If so, I don't mind putting it together.

[ViralBShah]

Would the right way here be to package SLATEC up in BinaryBuilder and make it available for use in SpecialFunctions.jl?

[stevengj]

SLATEC has > 1400 functions. Do we want them all? A lot of this seems useless.

[ViralBShah]

I tried to pull just the Bessel functions, but there were dependencies all over the place. I suppose I can find a set that works.

[inkydragon]

Maybe a better choice is "switch to Bessels.jl for real bessel functions" (JuliaMath/SpecialFunctions.jl#409)

Benchmark

julia v1.6

Bessels.jl is 200~300 times faster than SpecialFunctions.jl.

julia> # besselk0
julia> @btime SpecialFunctions.besselk(0, 0.5)
  181.277 ns (1 allocation: 16 bytes)
0.9244190712276656
julia> @btime Bessels.besselk0(0.5)
  0.700 ns (0 allocations: 0 bytes)
0.9244190712276659

julia> # besselk0x
julia> @btime SpecialFunctions.besselkx(0, 0.5)
  261.032 ns (1 allocation: 16 bytes)
1.5241093857739092
julia> @btime Bessels.besselk0x(0.5)
  0.700 ns (0 allocations: 0 bytes)
1.5241093857739096

julia> # besselk1
julia> @btime SpecialFunctions.besselk(1, 0.5)
  189.198 ns (1 allocation: 16 bytes)
1.6564411200033007
julia> @btime Bessels.besselk1(0.5)
  0.700 ns (0 allocations: 0 bytes)
1.656441120003301

julia> # besselk1x
julia> @btime SpecialFunctions.besselkx(1, 0.5)
  270.820 ns (1 allocation: 16 bytes)
2.7310097082117855
julia> @btime Bessels.besselk1x(0.5)
  0.900 ns (0 allocations: 0 bytes)
2.731009708211786

[ViralBShah]

Should we close this and note so in the README? @oscardssmith?

[oscardssmith]

Yeah. These benchmarks are wrong (the functions are getting constant folded), but the idea is correct.

julia> @btime besselk0(x) setup=x=rand()
  11.271 ns (0 allocations: 0 bytes)
julia> @btime besselk0x(x) setup=x=rand()
  14.870 ns (0 allocations: 0 bytes)
julia> @btime besselk1(x) setup=x=rand()
  11.401 ns (0 allocations: 0 bytes)
julia> @btime besselk1x(x) setup=x=rand()
  16.012 ns (0 allocations: 0 bytes)