Martín et al. (1992) provides two simple formulas for approximating $\mathrm{Ai}\left(x\right)$ to graphical accuracy, one for $$, the other for $$.
Moshier (1989, §6.14) provides minimax rational approximations for calculating $\mathrm{Ai}\left(x\right)$, ${\mathrm{Ai}}^{\prime}\left(x\right)$, $\mathrm{Bi}\left(x\right)$, ${\mathrm{Bi}}^{\prime}\left(x\right)$. They are in terms of the variable $\zeta $, where $\zeta =\frac{2}{3}{x}^{3/2}$ when $x$ is positive, $\zeta =\frac{2}{3}{(-x)}^{3/2}$ when $x$ is negative, and $\zeta =0$ when $x=0$. The approximations apply when $$, that is, when $$ or $$. The precision in the coefficients is 21S.
These expansions are for real arguments $x$ and are supplied in sets of four for each function, corresponding to intervals $$, $a\le x\le 0$, $0\le x\le b$, $$. The constants $a$ and $b$ are chosen numerically, with a view to equalizing the effort required for summing the series.
Németh (1992, Chapter 8) covers $\mathrm{Ai}\left(x\right)$, ${\mathrm{Ai}}^{\prime}\left(x\right)$, $\mathrm{Bi}\left(x\right)$, ${\mathrm{Bi}}^{\prime}\left(x\right)$, and integrals ${\int}_{0}^{x}\mathrm{Ai}\left(t\right)dt$, ${\int}_{0}^{x}\mathrm{Bi}\left(t\right)dt$, ${\int}_{0}^{x}{\int}_{0}^{v}\mathrm{Ai}\left(t\right)dtdv$, ${\int}_{0}^{x}{\int}_{0}^{v}\mathrm{Bi}\left(t\right)dtdv$ (see also (9.10.20) and (9.10.21)). The Chebyshev coefficients are given to 15D. Chebyshev coefficients are also given for expansions of the second and higher (real) zeros of $\mathrm{Ai}\left(x\right)$, ${\mathrm{Ai}}^{\prime}\left(x\right)$, $\mathrm{Bi}\left(x\right)$, ${\mathrm{Bi}}^{\prime}\left(x\right)$, again to 15D.
Razaz and Schonfelder (1980) covers $\mathrm{Ai}\left(x\right)$, ${\mathrm{Ai}}^{\prime}\left(x\right)$, $\mathrm{Bi}\left(x\right)$, ${\mathrm{Bi}}^{\prime}\left(x\right)$. The Chebyshev coefficients are given to 30D.
Corless et al. (1992) describe a method of approximation based on subdividing $\mathrm{\u2102}$ into a triangular mesh, with values of $\mathrm{Ai}\left(z\right)$, ${\mathrm{Ai}}^{\prime}\left(z\right)$ stored at the nodes. $\mathrm{Ai}\left(z\right)$ and ${\mathrm{Ai}}^{\prime}\left(z\right)$ are then computed from Taylor-series expansions centered at one of the nearest nodes. The Taylor coefficients are generated by recursion, starting from the stored values of $\mathrm{Ai}\left(z\right)$, ${\mathrm{Ai}}^{\prime}\left(z\right)$ at the node. Similarly for $\mathrm{Bi}\left(z\right)$, ${\mathrm{Bi}}^{\prime}\left(z\right)$.
MacLeod (1994) supplies Chebyshev-series expansions to cover $\mathrm{Gi}\left(x\right)$ for $$ and $\mathrm{Hi}\left(x\right)$ for $$. The Chebyshev coefficients are given to 20D.