From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/float/airyf.c | 377 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 377 insertions(+) create mode 100644 libm/float/airyf.c (limited to 'libm/float/airyf.c') 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 + +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); +} -- cgit v1.2.3