SpecialFunctions.jl Successor

Currently, the definitive Julia package for special functions is SpecialFunctions.jl, but it’s now being signalled that PRs adding functionality to it should be taken elsewhere because the package is mostly just for wrapping openspecfun.

Is anyone in the community interested in collaborating on a successor to SpecialFunctions.jl? To start, it could depend on SpecialFunctions.jl for a lot of functionality, but it’d be nice for it to eventually be pure-julia. I’m no expert on this sort of stuff, but I do have a copy of Numerical Recipes lying around. I’d be hesitant to go forward with this if we can’t find at least one collaborator who feels familiar with the implementation of special functions and has a good feel for numerical analysis.

What would such a successor package be called? I guess one possible path would be to move the openspecfun stuff out into a new repo called OpenSpecFunctions.jl, and keep the pure julia stuff in SpecialFunctions.jl and then get new maintainers in working on SpecialFunctions.jl, but this would require a lot of coordination and consent from the current maintainers.

Or it could be called something boring like SpecialFunctions2.jl.

5 Likes

I would personally really like efforts to be directed towards SIMD implementations of special functions.
Optimal SIMD and serial definitions are likely to be different, as the former ought to be branch-free.

2 Likes

A SimdSpecialFunctions.jl would indeed be quite nice. But in many ways, I see it as a higher priority that we at least have a package where serial special functions can live given that SpecialFunctions.jl appears to not be accepting any PRs anymore. There’s a bunch of work just sitting in the SpecialFunctions.jl PR queue just waiting to bit-rot.

Unfortunately, serial special functions aren’t a very sexy topic, they’re more of a utility that users expect to just exist.

3 Likes

Okay. Despite the name, SLEEFPirates is open for PRs that include implementations other than SLEEF ports, so perhaps a new library for SIMD implementations isn’t as needed.

EDIT: I’d also like to point out that @musm deserves all the credit for the the real work behind SLEEFPirates.

EDIT:
@viralbshah just showed me Libm.jl.

2 Likes

Note that Numerical Recipes has a very restrictive software license which means that no code from the book can be used. And basically you shouldn’t even look at the code. (Yes, that is ridiculous. But the code is pretty unreadable anyway.)

+1 for implementing special functions in Julia, but it’s (much) harder than it appears.

14 Likes

This seems to be a good, recent reference:

The Mathematical Function Computation Handbook, by Nelson Beebe

7 Likes

But the actual code for the library he describes is LGPL.

Implementing various special functions under an MIT license is a perfect idea for a GSOC (JSOC) project: it is well-defined, has immediate results, and valuable to the community.

I would recommend just starting from the literature, instead of looking at code with a license that is not compatible with a permissive one like MIT. Most libraries cite the relevant papers anyway.

Some knowledge of numerical analysis is helpful, but testing can be done by comparing with BigFloat calculations.

7 Likes

Love the idea of SIMD optimizations for special functions, but we can’t omit generic fallbacks. A big advantage of implementing these in Julia is automatic differentiation, and SIMD optimizations are a disadvantage for such purposes. Now, a SIMD-optimized dual numbers implementation…

4 Likes

SpecialFunctions is pretty much already covered by AD rules, so I don’t think it’ll do much here.

2 Likes

Good point. It’s worth keeping mind, though, for the day when someone implements pFq(a1,...,ap; b1,...,bp; x) (generalized hypergeometric function).

5 Likes

I looked at Gil, Segura, and Temme (2007) and it looks like another useful reference.

Some general thoughts about a native Julia package:

  1. I don’t think that fragmentation of “special” functions into smaller packages, depending on some natural constraints and details of implementation & testing, is necessarily a bad idea. Specifically, orthogonal families using recurrence relations and serial implementations of special functions like Bessel, Airy, and the error function just use different techniques and approaches.

    I would focus on (mostly) univariate (but possibly parametric) functions that use a clever approximation over \mathbb{R} to achieve a small relative error. In other words, the yucky formulas involving branches and tables of constants that come from continued fractions, Padé approximation, etc.

  2. Julia’s number types are generic, so it’s not clear what this package should target. I would suggest implementing methods for Float64, with the understanding that less and more precise types just fall back to this, with the implied extra cost or lack of possible precision. Interested contributors can extend these.

  3. The package should establish coding standard from the very beginning, including being explicit about errors we target (do we want custom variants, eg less precise functions?), testing against tables, BigFloat comparisons, random testing, etc; and importantly citation of sources and explicitly discouraging even looking at non-MIT-license compatible code. The primary sources should be papers.

  4. Ideally the package would work the following way: whenever someone needs a special function that is not available, they would just locate the relevant paper, code it, and make a PR, which is reviewed and merged in a reasonable timeframe. This way it would just channel contributions that may happen anyway. The benefit for contributors is that they get a code review (presumably they are implementing these functions because they are needed for their own work).

  5. Implementing functions within GSOC and similar should be considered — these are well-defined tasks and may serve as a great introduction to Julia but may require expertise that few students have.

  6. There should be 2–3 core contributors who have some numerical expertise and are willing to review and merge PRs.

14 Likes

One of my tasks is to study the fully generalized hypergeometric functions using geometry. I don’t claim to know every single thing necessary for this, but I am studying the general theory behind it and also pursuing this problem using some tricks. Fully generalized hypergeometric function is one of my long term goals.

3 Likes

I’m not sure why a package can’t be both a wrapper for openspecfun and also have native julia implementations. There are already a bunch of native implementations in there, and it seems simpler/less confusing to keep things together.

2 Likes

As far as my research on this goes, I will be creating my own repositories and not do GSOC or contribute to existing repositories. I work from first principles and do a completely original implementation from scratch, and will not want the hassle of dealing with a committee (as an unpaid volunteer, I want 100% full control of my repo). Also, my default license will probably be AGPL.

1 Like

Mmmm, if I may be so bold, let me say one thing… If some of the people around here were willing to act as a committee in any endeavour of mine and offer their expertise to me for free I would be on cloud nine, to be quite honest :smiley: . I used to think a bit like you in this respect, but then I discovered open source. I’m not interested in being in full control anymore, just want to share.

But back on topic, I agree a SpecialFunctions replacement is a key piece in the ecosystem, and deserves a good deal of attention to make it top notch and remain maintained. I’ve used it more that once in the past and found it very good. I didn’t even realize it was not pure Julia.

2 Likes

I’m not opposed to receiving guidance and free help from the community, I just dont like being required to answer to a committee when I’m a volunteer, as an obligation. Service obligations like that require payment.

I don’t think anyone expects that from you so no need to worry :). And if they do, what leverage do they have heh.

9 Likes

Apparently the wrapper for openspecfun (SpecialFunctions.jl) just wants to remain that.

Whether to separate wrapper and native libraries is a matter of style/preference. Personally I find it cleaner for many reasons, one of which is versioning: for a wrapper package you mostly follow changes in the library that is being wrapped, and aim for feature completeness, while a native library is more open.

Given how lightweight Julia packages are, I usually keep things modular unless there is a reason not to.

1 Like

Starting from first principles requires starting from scratch.