1 files changed, 377 insertions, 0 deletions
diff --git a/libm/float/airyf.c b/libm/float/airyf.c
new file mode 100644
index 000000000..a84a5c861
--- /dev/null
+++ b/libm/float/airyf.c
@@ -0,0 +1,377 @@
+/*							airy.c
+ *
+ *	Airy function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, ai, aip, bi, bip;
+ * int airyf();
+ *
+ * airyf( x, _&ai, _&aip, _&bi, _&bip );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Solution of the differential equation
+ *
+ *	y"(x) = xy.
+ *
+ * The function returns the two independent solutions Ai, Bi
+ * and their first derivatives Ai'(x), Bi'(x).
+ *
+ * Evaluation is by power series summation for small x,
+ * by rational minimax approximations for large x.
+ *
+ *
+ *
+ * ACCURACY:
+ * Error criterion is absolute when function <= 1, relative
+ * when function > 1, except * denotes relative error criterion.
+ * For large negative x, the absolute error increases as x^1.5.
+ * For large positive x, the relative error increases as x^1.5.
+ *
+ * Arithmetic  domain   function  # trials      peak         rms
+ * IEEE        -10, 0     Ai        50000       7.0e-7      1.2e-7
+ * IEEE          0, 10    Ai        50000       9.9e-6*     6.8e-7*
+ * IEEE        -10, 0     Ai'       50000       2.4e-6      3.5e-7
+ * IEEE          0, 10    Ai'       50000       8.7e-6*     6.2e-7*
+ * IEEE        -10, 10    Bi       100000       2.2e-6      2.6e-7
+ * IEEE        -10, 10    Bi'       50000       2.2e-6      3.5e-7
+ *
+ */
+/*							airy.c */
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+static float c1 = 0.35502805388781723926;
+static float c2 = 0.258819403792806798405;
+static float sqrt3 = 1.732050807568877293527;
+static float sqpii = 5.64189583547756286948E-1;
+extern float PIF;
+
+extern float MAXNUMF, MACHEPF;
+#define MAXAIRY 25.77
+
+/* Note, these expansions are for double precision accuracy;
+ * they have not yet been redesigned for single precision.
+ */
+static float AN[8] = {
+  3.46538101525629032477e-1,
+  1.20075952739645805542e1,
+  7.62796053615234516538e1,
+  1.68089224934630576269e2,
+  1.59756391350164413639e2,
+  7.05360906840444183113e1,
+  1.40264691163389668864e1,
+  9.99999999999999995305e-1,
+};
+static float AD[8] = {
+  5.67594532638770212846e-1,
+  1.47562562584847203173e1,
+  8.45138970141474626562e1,
+  1.77318088145400459522e2,
+  1.64234692871529701831e2,
+  7.14778400825575695274e1,
+  1.40959135607834029598e1,
+  1.00000000000000000470e0,
+};
+
+
+static float APN[8] = {
+  6.13759184814035759225e-1,
+  1.47454670787755323881e1,
+  8.20584123476060982430e1,
+  1.71184781360976385540e2,
+  1.59317847137141783523e2,
+  6.99778599330103016170e1,
+  1.39470856980481566958e1,
+  1.00000000000000000550e0,
+};
+static float APD[8] = {
+  3.34203677749736953049e-1,
+  1.11810297306158156705e1,
+  7.11727352147859965283e1,
+  1.58778084372838313640e2,
+  1.53206427475809220834e2,
+  6.86752304592780337944e1,
+  1.38498634758259442477e1,
+  9.99999999999999994502e-1,
+};
+
+static float BN16[5] = {
+-2.53240795869364152689e-1,
+ 5.75285167332467384228e-1,
+-3.29907036873225371650e-1,
+ 6.44404068948199951727e-2,
+-3.82519546641336734394e-3,
+};
+static float BD16[5] = {
+/* 1.00000000000000000000e0,*/
+-7.15685095054035237902e0,
+ 1.06039580715664694291e1,
+-5.23246636471251500874e0,
+ 9.57395864378383833152e-1,
+-5.50828147163549611107e-2,
+};
+
+static float BPPN[5] = {
+ 4.65461162774651610328e-1,
+-1.08992173800493920734e0,
+ 6.38800117371827987759e-1,
+-1.26844349553102907034e-1,
+ 7.62487844342109852105e-3,
+};
+static float BPPD[5] = {
+/* 1.00000000000000000000e0,*/
+-8.70622787633159124240e0,
+ 1.38993162704553213172e1,
+-7.14116144616431159572e0,
+ 1.34008595960680518666e0,
+-7.84273211323341930448e-2,
+};
+
+static float AFN[9] = {
+-1.31696323418331795333e-1,
+-6.26456544431912369773e-1,
+-6.93158036036933542233e-1,
+-2.79779981545119124951e-1,
+-4.91900132609500318020e-2,
+-4.06265923594885404393e-3,
+-1.59276496239262096340e-4,
+-2.77649108155232920844e-6,
+-1.67787698489114633780e-8,
+};
+static float AFD[9] = {
+/* 1.00000000000000000000e0,*/
+ 1.33560420706553243746e1,
+ 3.26825032795224613948e1,
+ 2.67367040941499554804e1,
+ 9.18707402907259625840e0,
+ 1.47529146771666414581e0,
+ 1.15687173795188044134e-1,
+ 4.40291641615211203805e-3,
+ 7.54720348287414296618e-5,
+ 4.51850092970580378464e-7,
+};
+
+static float AGN[11] = {
+  1.97339932091685679179e-2,
+  3.91103029615688277255e-1,
+  1.06579897599595591108e0,
+  9.39169229816650230044e-1,
+  3.51465656105547619242e-1,
+  6.33888919628925490927e-2,
+  5.85804113048388458567e-3,
+  2.82851600836737019778e-4,
+  6.98793669997260967291e-6,
+  8.11789239554389293311e-8,
+  3.41551784765923618484e-10,
+};
+static float AGD[10] = {
+/*  1.00000000000000000000e0,*/
+  9.30892908077441974853e0,
+  1.98352928718312140417e1,
+  1.55646628932864612953e1,
+  5.47686069422975497931e0,
+  9.54293611618961883998e-1,
+  8.64580826352392193095e-2,
+  4.12656523824222607191e-3,
+  1.01259085116509135510e-4,
+  1.17166733214413521882e-6,
+  4.91834570062930015649e-9,
+};
+
+static float APFN[9] = {
+  1.85365624022535566142e-1,
+  8.86712188052584095637e-1,
+  9.87391981747398547272e-1,
+  4.01241082318003734092e-1,
+  7.10304926289631174579e-2,
+  5.90618657995661810071e-3,
+  2.33051409401776799569e-4,
+  4.08718778289035454598e-6,
+  2.48379932900442457853e-8,
+};
+static float APFD[9] = {
+/*  1.00000000000000000000e0,*/
+  1.47345854687502542552e1,
+  3.75423933435489594466e1,
+  3.14657751203046424330e1,
+  1.09969125207298778536e1,
+  1.78885054766999417817e0,
+  1.41733275753662636873e-1,
+  5.44066067017226003627e-3,
+  9.39421290654511171663e-5,
+  5.65978713036027009243e-7,
+};
+
+static float APGN[11] = {
+-3.55615429033082288335e-2,
+-6.37311518129435504426e-1,
+-1.70856738884312371053e0,
+-1.50221872117316635393e0,
+-5.63606665822102676611e-1,
+-1.02101031120216891789e-1,
+-9.48396695961445269093e-3,
+-4.60325307486780994357e-4,
+-1.14300836484517375919e-5,
+-1.33415518685547420648e-7,
+-5.63803833958893494476e-10,
+};
+static float APGD[11] = {
+/*  1.00000000000000000000e0,*/
+  9.85865801696130355144e0,
+  2.16401867356585941885e1,
+  1.73130776389749389525e1,
+  6.17872175280828766327e0,
+  1.08848694396321495475e0,
+  9.95005543440888479402e-2,
+  4.78468199683886610842e-3,
+  1.18159633322838625562e-4,
+  1.37480673554219441465e-6,
+  5.79912514929147598821e-9,
+};
+
+#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
+
+float polevlf(float, float *, int);
+float p1evlf(float, float *, int);
+float sinf(float), cosf(float), expf(float), sqrtf(float);
+
+int airyf( float xx, float *ai, float *aip, float *bi, float *bip )
+{
+float x, z, zz, t, f, g, uf, ug, k, zeta, theta;
+int domflg;
+
+x = xx;
+domflg = 0;
+if( x > MAXAIRY )
+	{
+	*ai = 0;
+	*aip = 0;
+	*bi = MAXNUMF;
+	*bip = MAXNUMF;
+	return(-1);
+	}
+
+if( x < -2.09 )
+	{
+	domflg = 15;
+	t = sqrtf(-x);
+	zeta = -2.0 * x * t / 3.0;
+	t = sqrtf(t);
+	k = sqpii / t;
+	z = 1.0/zeta;
+	zz = z * z;
+	uf = 1.0 + zz * polevlf( zz, AFN, 8 ) / p1evlf( zz, AFD, 9 );
+	ug = z * polevlf( zz, AGN, 10 ) / p1evlf( zz, AGD, 10 );
+	theta = zeta + 0.25 * PIF;
+	f = sinf( theta );
+	g = cosf( theta );
+	*ai = k * (f * uf - g * ug);
+	*bi = k * (g * uf + f * ug);
+	uf = 1.0 + zz * polevlf( zz, APFN, 8 ) / p1evlf( zz, APFD, 9 );
+	ug = z * polevlf( zz, APGN, 10 ) / p1evlf( zz, APGD, 10 );
+	k = sqpii * t;
+	*aip = -k * (g * uf + f * ug);
+	*bip = k * (f * uf - g * ug);
+	return(0);
+	}
+
+if( x >= 2.09 )	/* cbrt(9) */
+	{
+	domflg = 5;
+	t = sqrtf(x);
+	zeta = 2.0 * x * t / 3.0;
+	g = expf( zeta );
+	t = sqrtf(t);
+	k = 2.0 * t * g;
+	z = 1.0/zeta;
+	f = polevlf( z, AN, 7 ) / polevlf( z, AD, 7 );
+	*ai = sqpii * f / k;
+	k = -0.5 * sqpii * t / g;
+	f = polevlf( z, APN, 7 ) / polevlf( z, APD, 7 );
+	*aip = f * k;
+
+	if( x > 8.3203353 )	/* zeta > 16 */
+		{
+		f = z * polevlf( z, BN16, 4 ) / p1evlf( z, BD16, 5 );
+		k = sqpii * g;
+		*bi = k * (1.0 + f) / t;
+		f = z * polevlf( z, BPPN, 4 ) / p1evlf( z, BPPD, 5 );
+		*bip = k * t * (1.0 + f);
+		return(0);
+		}
+	}
+
+f = 1.0;
+g = x;
+t = 1.0;
+uf = 1.0;
+ug = x;
+k = 1.0;
+z = x * x * x;
+while( t > MACHEPF )
+	{
+	uf *= z;
+	k += 1.0;
+	uf /=k;
+	ug *= z;
+	k += 1.0;
+	ug /=k;
+	uf /=k;
+	f += uf;
+	k += 1.0;
+	ug /=k;
+	g += ug;
+	t = fabsf(uf/f);
+	}
+uf = c1 * f;
+ug = c2 * g;
+if( (domflg & 1) == 0 )
+	*ai = uf - ug;
+if( (domflg & 2) == 0 )
+	*bi = sqrt3 * (uf + ug);
+
+/* the deriviative of ai */
+k = 4.0;
+uf = x * x/2.0;
+ug = z/3.0;
+f = uf;
+g = 1.0 + ug;
+uf /= 3.0;
+t = 1.0;
+
+while( t > MACHEPF )
+	{
+	uf *= z;
+	ug /=k;
+	k += 1.0;
+	ug *= z;
+	uf /=k;
+	f += uf;
+	k += 1.0;
+	ug /=k;
+	uf /=k;
+	g += ug;
+	k += 1.0;
+	t = fabsf(ug/g);
+	}
+
+uf = c1 * f;
+ug = c2 * g;
+if( (domflg & 4) == 0 )
+	*aip = uf - ug;
+if( (domflg & 8) == 0 )
+	*bip = sqrt3 * (uf + ug);
+return(0);
+}