Consider replacing SpecialFunctions.jl by pure Julia alternatives

[juliohm]

I've noticed that DualNumbers.jl checks if the function is available in SpecialFunctions.jl to wrap the definition with the Dual type. More recently, initiatives such as Bessels.jl have gained traction. They are pure Julia implementations of special functions that are faster and differentiable. What are your thoughts on replacing SpecialFunctions.jl by Bessels.jl and similar efforts to clean up the dependency tree downstream?

DualNumbers.jl is an indirect dependency of Distributions.jl which is pretty much everywhere in the ecosystem. So installing Distributions.jl nowadays comes with the downside of installing external _jll libraries and a ton of dependencies that could be avoided with Bessels.jl and friends.

I've already migrated all my packages from SpecialFunctions.jl to Bessels.jl because they provide all the bessel functions and the gamma (non-exported) functions I need without a single dependency.

cc: @heltonmc @cgeoga @oscardssmith

[dlfivefifty]

Looks as though Bessels.jl has a long way to go before its a full replacement: e.g. complex numbers

I'm surprised it's restricted to specific floating point types (Float32, Float64, etc.), eg, why not support BigFloat? I realise this is much more complicated but I think if one is going through all the effort its worth going that extra mile

[heltonmc]

I'll have to admit I'm not very familiar with this package or its downstream dependents and it is true that you can't propagate Dual numbers. That wasn't really a priority as taking derivatives are best handled by identities (unless wrt to order e.g., BesselK.jl).

If support for Complex numbers is really what's holding back adoption then I can focus on that. I wouldn't say it has a long way necessarily because all the algorithms can be used for complex numbers I just haven't really tested where best to make those branches. Though, the Amos routines that SpecialFunctions.jl use just develop routines for the modified bessel functions and use identities to piece together all the other functions. That is definitely the simplest but far from optimal.

Higher precision is definitely on the radar and we already have cutoffs for most regions. Though, I will say I don't really have intention to support BigFloat precision maybe like Float128 or Double64 at most. The algorithms even for Double64 and Float64 are pretty different.... so it takes some time to have these different implementations.