summaryrefslogtreecommitdiff
path: root/libm/float/airyf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/float/airyf.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/airyf.c')
-rw-r--r--libm/float/airyf.c377
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);
+}