[RFC]: add `airy` to compute Airy functions

[rreusser]

Description

Airy functions are a common special function in scientific applications, particularly wave phenomena. This RFC proposes to add a single package to compute for real-valued inputs the Airy function of the first kind, Ai(x), the second kind, Bi(x), and respective derivatives, Ai'(x) and Bi'(x).

The four outputs share a lot of computation, so for similar reasons and in the style of @stdlib/math/base/special/ellipj, perhaps the structure that makes sense is a single package, airy with exports:

airy(x)           // [ Ai(x), Ai'(x), Bi(x), Bi'(x) ]
airy.ai( x )      // Ai(x)
airy.aip( x )     // Ai'(x)
airy.bi( x )      // Bi(x)
airy.bip( x )     // Bi'(x)
airy.assign( x, out, stride, offset )

Exports ai, aip, bi, and bip would be convenience functions that compute all four values and throw away three.

Related Issues

No response

Questions

No.

Other

See also:

scipy.special.airy
cephes airy

Checklist

• I have read and understood the Code of Conduct.
• Searched for existing issues and pull requests.
• The issue name begins with RFC:.

[kgryte]

@rreusser Yes, this makes sense to me and would be a welcome addition!

[kgryte]

@rreusser Are you thinking to follow the SciPy implementation or is there some other reference implementation (e.g., Boost) which serves as the gold standard?

[kgryte]

SciPy's implementation in xsf: https://github.com/scipy/xsf/blob/main/include/xsf/airy.h

[rreusser]

It looks to me like:

• cephes: relies heavily on rational polynomial approximations. Simple to port. (full disclosure: poked around and already ported it, though that shouldn't preclude making the best choice)
• xsf: relies on slightly more involved iteration, but looks similarly highly straightforward to port.
• boost: farms out to cylindrical bessel functions. Maybe alright to port and plus then you get some other functions, but more involved dependencies to implement and trace.

I don't have a strong opinion myself. My preference would be for cephes or xsf, then boost, mainly just since boost seems very, very well done, but in a way that requires a lot more untangling of dependencies and logic and probably a couple/few precursor packages. (I think it was for bessel where functions keep calling themselves with differently-typed arguments until templating sifts through it and selects the proper implementation. C++ certainly is certainly a higher-effort candidate for translation.)

Just one note, am I missing something that it looks to me like scipy depends on cephes and not xsf? See scipy.airy notes:

For real z in [-10, 10], the computation is carried out by calling the Cephes [1] airy routine, which uses power series summation for small z and rational minimax approximations for large z.

Outside this range, the AMOS [2] zairy and zbiry routines are employed.

Perhaps the notes are outdated?

[kgryte]

xsf is what SciPy is moving toward. They want a header-only implementation of scalar kernels to allow for GPU compilation via nvcc. Looking at xsf, it seems they do use both the Cephes and the Amos implementations based on value range. Both look straightforward to port.

[rreusser]

Ah, that makes sense. How about xsf then?

[kgryte]

I am fine with using the xsf implementation, which is ultimately derived (and occasionally improved) from Cephes.

[rreusser]

Well porting xsf took… ten minutes. 👍

[kgryte]

@rreusser Well, you are a pro.