Module airy 

The Airy functions $\mathrm{Ai}$, $\mathrm{Bi}$ and their derivatives.

Airy is Bessel at order $\pm 1/3$ and $\pm 2/3$, with $\zeta = \tfrac{2}{3}\lvert x\rvert^{3/2}$. That is not a shortcut anyone invented here. It is how Boost.Math’s Airy functions are written, and the reason the fractional-order machinery in bessel_nu and bessel_ik_nu had to exist first.

\begin{aligned}
x < 0:\quad \mathrm{Ai} &= \tfrac{\sqrt{-x}}{3}\left(J_{1/3} + J_{-1/3}\right), &
            \mathrm{Bi} &= \sqrt{\tfrac{-x}{3}}\left(J_{-1/3} - J_{1/3}\right) \\
x > 0:\quad \mathrm{Ai} &= \tfrac{1}{\pi}\sqrt{\tfrac{x}{3}}\,K_{1/3}, &
            \mathrm{Bi} &= \sqrt{\tfrac{x}{3}}\left(I_{-1/3} + I_{1/3}\right)
\end{aligned}

with the derivatives the same shapes at order $2/3$. Note $\mathrm{Ai}$ on the positive axis goes through $K$ rather than the $I$ difference, and Boost’s comment says why: “the accuracy is horrible as we’re subtracting two very large values”. $\mathrm{Bi}$ uses the $I$ sum, which has no such problem.

Order 1/3 is where the rotation stops costing anything

Both branches need the negative order as well as the positive one, and a second Bessel evaluation is the expensive way to get it. One pass suffices: $J_{-\nu} = J_\nu\cos\nu\pi - Y_\nu\sin\nu\pi$ and $I_{-\nu} = I_\nu + \tfrac{2}{\pi}\sin(\nu\pi)K_\nu$, and both kernels return the pair.

At thirds those trigonometric factors are exact: $\cos(\pi/3) = 1/2$, $\cos(2\pi/3) = -1/2$, and $\sin(\pi/3) = \sin(2\pi/3) = \sqrt3/2$. So the rotation costs two multiplies and no transcendental at all. This is the whole of what a Thirds specialisation can buy (there is no cheaper algorithm for a third-order Bessel function in any library), and it is bought here rather than in the order dispatch.

Accuracy is set by zeta, which is why the scaled form exists

Measured against mpmath (notes/special/tools/model_airy.py), the error of the unscaled functions grows linearly in $\zeta$ on both sides of the origin, and for two different reasons that arrive at the same number:

x	mechanism	measured
-10	phase: $\mathrm{ulp}(\zeta)$ of argument error becomes phase error	21 eps
-300	same	3.7e3 eps
-1e4	same	3.2e5 eps
+10	the $e^{-\zeta}$ factor costs $\zeta\varepsilon/2$ relative	3.2 eps
+100	same	684 eps

On the negative axis that is irreducible without a two-word $\zeta$, and Boost has the same exposure. On the positive axis it is not: with SCALED the exponential is never formed, because bessel_ik_real already returns $e^{\zeta}K$ natively. The scaled pair holds 1-3 eps everywhere and, being free of the exponential, also has no range limit. $\zeta$ passes 710 at $x \approx 104$, where the unscaled $\mathrm{Ai}$ goes subnormal (zero by 108) and $\mathrm{Bi}$ overflows.

Boost ships no scaled Airy. SciPy does, as airye, and for exactly this reason.

Functions

airy_impl

$(\mathrm{Ai}(x),\; \mathrm{Ai}'(x),\; \mathrm{Bi}(x),\; \mathrm{Bi}'(x))$.