Scope of Bessels.jl and inclusion of other special functions (e.g. gamma, airy, scorer)

[heltonmc]

EDIT: This issue started with discussion on the gamma implementation but has moved to discussing the scope of this package and inclusion of other special functions or splitting into smaller packages. Please leave comments on what functions you would like to be included here or if you would like smaller packages.

SPLIT FROM https://github.com/heltonmc/Bessels.jl/pull/26#issuecomment-1200481143

A small note: I was perusing your gamma implementation and I think there are a couple small improvements that could be made, assuming they don't cost precision somehow. For one, if I change the v = x^(0.5*x - 0.25) to v = exp((0.5*x-0.25)*log(x)) in large_gamma, I see a big speedup. And it doesn't seem like the log would make an issue, because you'll never hit that branch when x is zero.

Second, if there's a reason not to do that, then an intermediate method like this seems to be faster when x > 11.5 but it isn't too too large:

function smallish_gamma(x)
  _x = min(x, x-ceil(x-11.5))
  _g = small_gamma(_x)
  while _x < x # use the fact that gamma(x+1) == x*gamma(x)
    _g *= _x
    _x += one(_x)
  end
  _g
end

You could then branch on x < 30.0 or something instead and at least hit the large_gamma branch less often.

I haven't checked very hard to see if these tricks introduce new issues, but I wouldn't think that they do.

Originally posted by @cgeoga in https://github.com/heltonmc/Bessels.jl/issues/26#issuecomment-1200481143

[heltonmc]

Thank you for looking at this @cgeoga !!! Let's discuss here to try and keep the other thread focused on bessely.

I'm very happy to take PRs on this in general which we can iron out potentially a better implementation of this. As mentioned https://github.com/heltonmc/Bessels.jl/pull/26#issuecomment-1200514696 we would definitely need to test accuracy for this. The good thing about 1d functions is we can exhaustively test all the ranges. I would say for a function as important as gamma here keeping the max error below ~5 ULP would be ideal. I need to move the gamma function into its own file as well to keep it separate and add tests as well. But as mentioned hopefully the powers are faster in 1.9, though I just did a quick test on a very old 1.9 branch and found it was much slower than 1.7, so keeping track of this by julia version would also be very helpful.

[heltonmc]

On a quick first test it looks like the first suggestion increased the mean ULP from 0.26595 to 1.98555 and the maximum ULP from 3 to 571.0 at x=162.81202133604663.

[oscardssmith]

One thing worth noting is that although these techniques won't work for Float64, for Float32, there are a lot of options opened since you can do the internals in Float64 to avoid loosing precision.

[cgeoga]

Makes sense about moving the issues and apologies for the clutter. Also, wow about the additional error. I assumed it would be a little bit worse, but that is certainly more than I expected.

You've mentioned these back-of-the-envelope errors in terms of ULPs before---is there some easy demo you can point me to so I can figure out how to do that kind of error testing myself? I should learn how to do this anyway, and then I can sort of pre-screen some of these ideas before posting.

[heltonmc]

I haven't tried it out but see Oscar's https://github.com/JuliaMath/FunctionAccuracyTests.jl.
Though in this case I just tested against SpecialFunctions.jl

import SpecialFunctions
julia> x = rand(10000)*170
julia> maximum(abs.(1 .- SpecialFunctions.gamma.(BigFloat.(x)) ./ gamma.(x)))/ eps(Float64)

[oscardssmith]

If you want more details on my thought process here, my Juliacon talk was on floating point math implementation and testing. https://www.youtube.com/watch?v=-tFS4PUHqT0

[cgeoga]

Nice, thank you both! Will do some studying.

[juliohm]

Can you guys comment a bit more on the current gamma implementation defined in Bessels.jl? Is it ready to use? Why is it not exported? I have a package that depends on both besselk and gamma from SpecialFunctions.jl and would like to replace this dependency by Bessels.jl that is pure Julia with zero dependencies. However, I would like to understand the future plans better before making this migration.

[heltonmc]

The function is definitely ready to use. We use it in practically every function in this library which is why we rolled our own implementation that is much faster. It has a max ULP around 2.5-3. But it does have some limitations as it only works for real arguments and for single and double precision. So if you need it for complex or higher precision it won't work.

The reason it is not exported is because this is probably not the best package for a native gamma function implementation. There have been some discussion about moving it over to SpecialFunctions.jl but my plan was to actually create a lightweight Gamma.jl package that just has julia implementations of gamma and similar functions. So I'd be hesitant to export it here as it is likely that we will move it out of the package in the future.

[cgeoga]

+1 for a Gamma.jl package, that would be great.

[juliohm]

Awesome @heltonmc , thank you for the prompt reply. I am also +1 for Gamma.jl or maybe a new package that wraps all these special functions in pure Julia, maybe SpecialFuns.jl?

[oscardssmith]

the argument for splitting them up is that a lot of libraries want to use gamma functions without the precompilation time of loading everything.

[juliohm]

Thank you @oscardssmith for adding the comment. Is the compilation time a real issue for Bessels.jl or is this another case of premature modularization? I feel that many of these package splits in Julia are more harmful than helpful, specially when it comes to maintenance. The collateral effect of having many packages is that some of them die after a while without maintainers. Working on the same repo with the same code style, multiple maintainers, etc. pays off in my opinion.

[oscardssmith]

Well using SpecialFunctions is already .35 seconds (and most of the code is just calling C functions).

[heltonmc]

@juliohm that was my general idea to have a higher level package that wraps all of the mathematical functions outside of Base. I'd like to see more of a canonicalization of the math functions in Julia (particularly the special functions). In my opinion there are many packages that do similar things but the quality is variable and sometimes hard to discern for new users. And they are often put together quickly for a specific application, registered, and then never maintained after…. A python user is immediately going to SciPy or mpmath and they have some sense of quality and stability which I think is a good thing.

Part of the problem with SpecialFunctions.jl is it is actually hard to find maintainers and get things merged where the scope is not so clearly defined (and maybe nobody feels any kind of authority?). One example that I think is concerning is the most recent PR and issue about the fresnel integrals which existed in a different package where a first time contributor to the language submitted a PR to add that functionality that was carefully reviewed by a primary contributor to the language but was decided that it is better in a separate package. So that contributor submitted to the other package but that package is unmaintained…. now we are point people on SpecialFuncitons.jl to an unmaintained package where the existing code base can be improved. So to your point that is definitely an issue with all of these smaller packages. Though, ideally these can be brought under a higher level package under an organization so we can continue to improve them.

So yes… what to do with gamma 🤷‍♂️ ? I’m happy for this package to have special functions but I start to become a bit worried about scope. Which is something I think SpecialFunctions.jl is not so clearly defined. Should we also include loggamma, beta, polygama, digamma, and their incomplete forms here? But I feel like then why not just start making Bessels.jl (or rename it) to just general special functions. I still think the right thing to do is try to improve the scattered special functions already registered in the general repository and improve them and then bring them under some higher organization and repository to standardize in some way. For example, I’ve been spending some time trying to improve Struve.jl so that hopefully we can have several similar quality libraries to piece together.

But I think I agree with your overall sentiment…. Though maybe I was hoping to improve the already existing registered with very defined scope that a higher level package can reexport. So new users can be a little more confident in the quality of the package they are using and having an easy interface with a single package. Though say if you wanted just the gamma function, you can then just depend on that package which has a very defined scope and very lightweight.

Really appreciate your thoughts and definitely welcome to other suggestions!