`besselk(x, 0im)` fails

[theogf]

besselk(x, 0im) should just return T(Inf) for any x but returns an AmosException instead. Here is a MWE:

using SpecialFunctions
besselk(2, 0im)
ERROR: AmosException with id 1: input error.
Stacktrace:
 [1] _besselk(nu::Float64, z::ComplexF64, kode::Int32)
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/bessel.jl:274
 [2] besselk
   @ ~/.julia/packages/SpecialFunctions/9pXme/src/bessel.jl:408 [inlined]
 [3] besselk(nu::Int64, z::Complex{Int64})
   @ SpecialFunctions ~/.julia/packages/SpecialFunctions/9pXme/src/bessel.jl:586
 [4] top-level scope
   @ REPL[14]:1

Error noted in https://github.com/JuliaStats/Distributions.jl/pull/1504

[simonbyrne]

That seems reasonable, but it is complicated by the fact that besselk has a branch cut along the negative real axis, so the exact limit isn't well-defined (i.e. it's a singularity, but the arg depends on which direction you approach from).

Other functions have similar issues (e.g. log), and there are some conventions for dealing with this (https://people.freebsd.org/~das/kahan86branch.pdf).

[simonbyrne]

I think it might make sense to define it as complex(T(Inf), -zero(T)*real(z)*imag(z))?

[ViralBShah]

How does https://github.com/heltonmc/Bessels.jl do on some of these cases?

[theogf]

It looks like besselk in Bessels.jl does not allow for complex arguments...

[oscardssmith]

yeah. They're in scope, but not implemented yet. If anyone wants to start working on it, it would be a good addition.