lib/libm/airyf.c

   1 /*                                                      airy.c
   2  *
   3  *      Airy function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, ai, aip, bi, bip;
  10  * int airyf();
  11  *
  12  * airyf( x, _&ai, _&aip, _&bi, _&bip );
  13  *
  14  *
  15  *
  16  * DESCRIPTION:
  17  *
  18  * Solution of the differential equation
  19  *
  20  *      y"(x) = xy.
  21  *
  22  * The function returns the two independent solutions Ai, Bi
  23  * and their first derivatives Ai'(x), Bi'(x).
  24  *
  25  * Evaluation is by power series summation for small x,
  26  * by rational minimax approximations for large x.
  27  *
  28  *
  29  *
  30  * ACCURACY:
  31  * Error criterion is absolute when function <= 1, relative
  32  * when function > 1, except * denotes relative error criterion.
  33  * For large negative x, the absolute error increases as x^1.5.
  34  * For large positive x, the relative error increases as x^1.5.
  35  *
  36  * Arithmetic  domain   function  # trials      peak         rms
  37  * IEEE        -10, 0     Ai        50000       7.0e-7      1.2e-7
  38  * IEEE          0, 10    Ai        50000       9.9e-6*     6.8e-7*
  39  * IEEE        -10, 0     Ai'       50000       2.4e-6      3.5e-7
  40  * IEEE          0, 10    Ai'       50000       8.7e-6*     6.2e-7*
  41  * IEEE        -10, 10    Bi       100000       2.2e-6      2.6e-7
  42  * IEEE        -10, 10    Bi'       50000       2.2e-6      3.5e-7
  43  *
  44  */
  45 \f/*                                                     airy.c */
  46
  47 /*
  48 Cephes Math Library Release 2.2: June, 1992
  49 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  50 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  51 */
  52
  53 #include "mconf.h"
  54
  55 static float c1 = 0.35502805388781723926;
  56 static float c2 = 0.258819403792806798405;
  57 static float sqrt3 = 1.732050807568877293527;
  58 static float sqpii = 5.64189583547756286948E-1;
  59 extern float PIF;
  60
  61 extern float MAXNUMF, MACHEPF;
  62 #define MAXAIRY 25.77
  63
  64 /* Note, these expansions are for double precision accuracy;
  65  * they have not yet been redesigned for single precision.
  66  */
  67 static float AN[8] = {
  68   3.46538101525629032477e-1,
  69   1.20075952739645805542e1,
  70   7.62796053615234516538e1,
  71   1.68089224934630576269e2,
  72   1.59756391350164413639e2,
  73   7.05360906840444183113e1,
  74   1.40264691163389668864e1,
  75   9.99999999999999995305e-1,
  76 };
  77 static float AD[8] = {
  78   5.67594532638770212846e-1,
  79   1.47562562584847203173e1,
  80   8.45138970141474626562e1,
  81   1.77318088145400459522e2,
  82   1.64234692871529701831e2,
  83   7.14778400825575695274e1,
  84   1.40959135607834029598e1,
  85   1.00000000000000000470e0,
  86 };
  87
  88
  89 static float APN[8] = {
  90   6.13759184814035759225e-1,
  91   1.47454670787755323881e1,
  92   8.20584123476060982430e1,
  93   1.71184781360976385540e2,
  94   1.59317847137141783523e2,
  95   6.99778599330103016170e1,
  96   1.39470856980481566958e1,
  97   1.00000000000000000550e0,
  98 };
  99 static float APD[8] = {
 100   3.34203677749736953049e-1,
 101   1.11810297306158156705e1,
 102   7.11727352147859965283e1,
 103   1.58778084372838313640e2,
 104   1.53206427475809220834e2,
 105   6.86752304592780337944e1,
 106   1.38498634758259442477e1,
 107   9.99999999999999994502e-1,
 108 };
 109
 110 static float BN16[5] = {
 111 -2.53240795869364152689e-1,
 112  5.75285167332467384228e-1,
 113 -3.29907036873225371650e-1,
 114  6.44404068948199951727e-2,
 115 -3.82519546641336734394e-3,
 116 };
 117 static float BD16[5] = {
 118 /* 1.00000000000000000000e0,*/
 119 -7.15685095054035237902e0,
 120  1.06039580715664694291e1,
 121 -5.23246636471251500874e0,
 122  9.57395864378383833152e-1,
 123 -5.50828147163549611107e-2,
 124 };
 125
 126 static float BPPN[5] = {
 127  4.65461162774651610328e-1,
 128 -1.08992173800493920734e0,
 129  6.38800117371827987759e-1,
 130 -1.26844349553102907034e-1,
 131  7.62487844342109852105e-3,
 132 };
 133 static float BPPD[5] = {
 134 /* 1.00000000000000000000e0,*/
 135 -8.70622787633159124240e0,
 136  1.38993162704553213172e1,
 137 -7.14116144616431159572e0,
 138  1.34008595960680518666e0,
 139 -7.84273211323341930448e-2,
 140 };
 141
 142 static float AFN[9] = {
 143 -1.31696323418331795333e-1,
 144 -6.26456544431912369773e-1,
 145 -6.93158036036933542233e-1,
 146 -2.79779981545119124951e-1,
 147 -4.91900132609500318020e-2,
 148 -4.06265923594885404393e-3,
 149 -1.59276496239262096340e-4,
 150 -2.77649108155232920844e-6,
 151 -1.67787698489114633780e-8,
 152 };
 153 static float AFD[9] = {
 154 /* 1.00000000000000000000e0,*/
 155  1.33560420706553243746e1,
 156  3.26825032795224613948e1,
 157  2.67367040941499554804e1,
 158  9.18707402907259625840e0,
 159  1.47529146771666414581e0,
 160  1.15687173795188044134e-1,
 161  4.40291641615211203805e-3,
 162  7.54720348287414296618e-5,
 163  4.51850092970580378464e-7,
 164 };
 165
 166 static float AGN[11] = {
 167   1.97339932091685679179e-2,
 168   3.91103029615688277255e-1,
 169   1.06579897599595591108e0,
 170   9.39169229816650230044e-1,
 171   3.51465656105547619242e-1,
 172   6.33888919628925490927e-2,
 173   5.85804113048388458567e-3,
 174   2.82851600836737019778e-4,
 175   6.98793669997260967291e-6,
 176   8.11789239554389293311e-8,
 177   3.41551784765923618484e-10,
 178 };
 179 static float AGD[10] = {
 180 /*  1.00000000000000000000e0,*/
 181   9.30892908077441974853e0,
 182   1.98352928718312140417e1,
 183   1.55646628932864612953e1,
 184   5.47686069422975497931e0,
 185   9.54293611618961883998e-1,
 186   8.64580826352392193095e-2,
 187   4.12656523824222607191e-3,
 188   1.01259085116509135510e-4,
 189   1.17166733214413521882e-6,
 190   4.91834570062930015649e-9,
 191 };
 192
 193 static float APFN[9] = {
 194   1.85365624022535566142e-1,
 195   8.86712188052584095637e-1,
 196   9.87391981747398547272e-1,
 197   4.01241082318003734092e-1,
 198   7.10304926289631174579e-2,
 199   5.90618657995661810071e-3,
 200   2.33051409401776799569e-4,
 201   4.08718778289035454598e-6,
 202   2.48379932900442457853e-8,
 203 };
 204 static float APFD[9] = {
 205 /*  1.00000000000000000000e0,*/
 206   1.47345854687502542552e1,
 207   3.75423933435489594466e1,
 208   3.14657751203046424330e1,
 209   1.09969125207298778536e1,
 210   1.78885054766999417817e0,
 211   1.41733275753662636873e-1,
 212   5.44066067017226003627e-3,
 213   9.39421290654511171663e-5,
 214   5.65978713036027009243e-7,
 215 };
 216
 217 static float APGN[11] = {
 218 -3.55615429033082288335e-2,
 219 -6.37311518129435504426e-1,
 220 -1.70856738884312371053e0,
 221 -1.50221872117316635393e0,
 222 -5.63606665822102676611e-1,
 223 -1.02101031120216891789e-1,
 224 -9.48396695961445269093e-3,
 225 -4.60325307486780994357e-4,
 226 -1.14300836484517375919e-5,
 227 -1.33415518685547420648e-7,
 228 -5.63803833958893494476e-10,
 229 };
 230 static float APGD[11] = {
 231 /*  1.00000000000000000000e0,*/
 232   9.85865801696130355144e0,
 233   2.16401867356585941885e1,
 234   1.73130776389749389525e1,
 235   6.17872175280828766327e0,
 236   1.08848694396321495475e0,
 237   9.95005543440888479402e-2,
 238   4.78468199683886610842e-3,
 239   1.18159633322838625562e-4,
 240   1.37480673554219441465e-6,
 241   5.79912514929147598821e-9,
 242 };
 243
 244 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
 245
 246 #ifdef ANSIC
 247 float polevlf(float, float *, int);
 248 float p1evlf(float, float *, int);
 249 float sinf(float), cosf(float), expf(float), sqrtf(float);
 250
 251 int airyf( float xx, float *ai, float *aip, float *bi, float *bip )
 252
 253 #else
 254 float polevlf(), p1evlf(), sinf(), cosf(), expf(), sqrtf();
 255
 256 int airyf( xx, ai, aip, bi, bip )
 257 double xx;
 258 float *ai, *aip, *bi, *bip;
 259 #endif
 260 {
 261 float x, z, zz, t, f, g, uf, ug, k, zeta, theta;
 262 int domflg;
 263
 264 x = xx;
 265 domflg = 0;
 266 if( x > MAXAIRY )
 267         {
 268         *ai = 0;
 269         *aip = 0;
 270         *bi = MAXNUMF;
 271         *bip = MAXNUMF;
 272         return(-1);
 273         }
 274
 275 if( x < -2.09 )
 276         {
 277         domflg = 15;
 278         t = sqrtf(-x);
 279         zeta = -2.0 * x * t / 3.0;
 280         t = sqrtf(t);
 281         k = sqpii / t;
 282         z = 1.0/zeta;
 283         zz = z * z;
 284         uf = 1.0 + zz * polevlf( zz, AFN, 8 ) / p1evlf( zz, AFD, 9 );
 285         ug = z * polevlf( zz, AGN, 10 ) / p1evlf( zz, AGD, 10 );
 286         theta = zeta + 0.25 * PIF;
 287         f = sinf( theta );
 288         g = cosf( theta );
 289         *ai = k * (f * uf - g * ug);
 290         *bi = k * (g * uf + f * ug);
 291         uf = 1.0 + zz * polevlf( zz, APFN, 8 ) / p1evlf( zz, APFD, 9 );
 292         ug = z * polevlf( zz, APGN, 10 ) / p1evlf( zz, APGD, 10 );
 293         k = sqpii * t;
 294         *aip = -k * (g * uf + f * ug);
 295         *bip = k * (f * uf - g * ug);
 296         return(0);
 297         }
 298
 299 if( x >= 2.09 ) /* cbrt(9) */
 300         {
 301         domflg = 5;
 302         t = sqrtf(x);
 303         zeta = 2.0 * x * t / 3.0;
 304         g = expf( zeta );
 305         t = sqrtf(t);
 306         k = 2.0 * t * g;
 307         z = 1.0/zeta;
 308         f = polevlf( z, AN, 7 ) / polevlf( z, AD, 7 );
 309         *ai = sqpii * f / k;
 310         k = -0.5 * sqpii * t / g;
 311         f = polevlf( z, APN, 7 ) / polevlf( z, APD, 7 );
 312         *aip = f * k;
 313
 314         if( x > 8.3203353 )     /* zeta > 16 */
 315                 {
 316                 f = z * polevlf( z, BN16, 4 ) / p1evlf( z, BD16, 5 );
 317                 k = sqpii * g;
 318                 *bi = k * (1.0 + f) / t;
 319                 f = z * polevlf( z, BPPN, 4 ) / p1evlf( z, BPPD, 5 );
 320                 *bip = k * t * (1.0 + f);
 321                 return(0);
 322                 }
 323         }
 324
 325 f = 1.0;
 326 g = x;
 327 t = 1.0;
 328 uf = 1.0;
 329 ug = x;
 330 k = 1.0;
 331 z = x * x * x;
 332 while( t > MACHEPF )
 333         {
 334         uf *= z;
 335         k += 1.0;
 336         uf /=k;
 337         ug *= z;
 338         k += 1.0;
 339         ug /=k;
 340         uf /=k;
 341         f += uf;
 342         k += 1.0;
 343         ug /=k;
 344         g += ug;
 345         t = fabsf(uf/f);
 346         }
 347 uf = c1 * f;
 348 ug = c2 * g;
 349 if( (domflg & 1) == 0 )
 350         *ai = uf - ug;
 351 if( (domflg & 2) == 0 )
 352         *bi = sqrt3 * (uf + ug);
 353
 354 /* the deriviative of ai */
 355 k = 4.0;
 356 uf = x * x/2.0;
 357 ug = z/3.0;
 358 f = uf;
 359 g = 1.0 + ug;
 360 uf /= 3.0;
 361 t = 1.0;
 362
 363 while( t > MACHEPF )
 364         {
 365         uf *= z;
 366         ug /=k;
 367         k += 1.0;
 368         ug *= z;
 369         uf /=k;
 370         f += uf;
 371         k += 1.0;
 372         ug /=k;
 373         uf /=k;
 374         g += ug;
 375         k += 1.0;
 376         t = fabsf(ug/g);
 377         }
 378
 379 uf = c1 * f;
 380 ug = c2 * g;
 381 if( (domflg & 4) == 0 )
 382         *aip = uf - ug;
 383 if( (domflg & 8) == 0 )
 384         *bip = sqrt3 * (uf + ug);
 385 return(0);
 386 }