1 /* airy.c
2 *
3 * Airy function
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, ai, aip, bi, bip;
10 * int airy();
11 *
12 * airy( 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 10000 1.6e-15 2.7e-16
38 * IEEE 0, 10 Ai 10000 2.3e-14* 1.8e-15*
39 * IEEE -10, 0 Ai' 10000 4.6e-15 7.6e-16
40 * IEEE 0, 10 Ai' 10000 1.8e-14* 1.5e-15*
41 * IEEE -10, 10 Bi 30000 4.2e-15 5.3e-16
42 * IEEE -10, 10 Bi' 30000 4.9e-15 7.3e-16
43 *
44 */
45 /* airy.c */
46
47 /*
48 * Cephes Math Library Release 2.8: June, 2000
49 * Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
50 */
51
52 #include "mconf.h"
53
54 static double c1 = 0.35502805388781723926;
55 static double c2 = 0.258819403792806798405;
56 static double sqrt3 = 1.732050807568877293527;
57 static double sqpii = 5.64189583547756286948E-1;
58
59 extern double MACHEP;
60
61 #ifdef UNK
62 #define MAXAIRY 25.77
63 #endif
64 #ifdef IBMPC
65 #define MAXAIRY 103.892
66 #endif
67 #ifdef MIEEE
68 #define MAXAIRY 103.892
69 #endif
70
71
72 static double AN[8] = {
73 3.46538101525629032477E-1,
74 1.20075952739645805542E1,
75 7.62796053615234516538E1,
76 1.68089224934630576269E2,
77 1.59756391350164413639E2,
78 7.05360906840444183113E1,
79 1.40264691163389668864E1,
80 9.99999999999999995305E-1,
81 };
82
83 static double AD[8] = {
84 5.67594532638770212846E-1,
85 1.47562562584847203173E1,
86 8.45138970141474626562E1,
87 1.77318088145400459522E2,
88 1.64234692871529701831E2,
89 7.14778400825575695274E1,
90 1.40959135607834029598E1,
91 1.00000000000000000470E0,
92 };
93
94 static double APN[8] = {
95 6.13759184814035759225E-1,
96 1.47454670787755323881E1,
97 8.20584123476060982430E1,
98 1.71184781360976385540E2,
99 1.59317847137141783523E2,
100 6.99778599330103016170E1,
101 1.39470856980481566958E1,
102 1.00000000000000000550E0,
103 };
104
105 static double APD[8] = {
106 3.34203677749736953049E-1,
107 1.11810297306158156705E1,
108 7.11727352147859965283E1,
109 1.58778084372838313640E2,
110 1.53206427475809220834E2,
111 6.86752304592780337944E1,
112 1.38498634758259442477E1,
113 9.99999999999999994502E-1,
114 };
115
116 static double BN16[5] = {
117 -2.53240795869364152689E-1,
118 5.75285167332467384228E-1,
119 -3.29907036873225371650E-1,
120 6.44404068948199951727E-2,
121 -3.82519546641336734394E-3,
122 };
123
124 static double BD16[5] = {
125 /* 1.00000000000000000000E0, */
126 -7.15685095054035237902E0,
127 1.06039580715664694291E1,
128 -5.23246636471251500874E0,
129 9.57395864378383833152E-1,
130 -5.50828147163549611107E-2,
131 };
132
133 static double BPPN[5] = {
134 4.65461162774651610328E-1,
135 -1.08992173800493920734E0,
136 6.38800117371827987759E-1,
137 -1.26844349553102907034E-1,
138 7.62487844342109852105E-3,
139 };
140
141 static double BPPD[5] = {
142 /* 1.00000000000000000000E0, */
143 -8.70622787633159124240E0,
144 1.38993162704553213172E1,
145 -7.14116144616431159572E0,
146 1.34008595960680518666E0,
147 -7.84273211323341930448E-2,
148 };
149
150 static double AFN[9] = {
151 -1.31696323418331795333E-1,
152 -6.26456544431912369773E-1,
153 -6.93158036036933542233E-1,
154 -2.79779981545119124951E-1,
155 -4.91900132609500318020E-2,
156 -4.06265923594885404393E-3,
157 -1.59276496239262096340E-4,
158 -2.77649108155232920844E-6,
159 -1.67787698489114633780E-8,
160 };
161
162 static double AFD[9] = {
163 /* 1.00000000000000000000E0, */
164 1.33560420706553243746E1,
165 3.26825032795224613948E1,
166 2.67367040941499554804E1,
167 9.18707402907259625840E0,
168 1.47529146771666414581E0,
169 1.15687173795188044134E-1,
170 4.40291641615211203805E-3,
171 7.54720348287414296618E-5,
172 4.51850092970580378464E-7,
173 };
174
175 static double AGN[11] = {
176 1.97339932091685679179E-2,
177 3.91103029615688277255E-1,
178 1.06579897599595591108E0,
179 9.39169229816650230044E-1,
180 3.51465656105547619242E-1,
181 6.33888919628925490927E-2,
182 5.85804113048388458567E-3,
183 2.82851600836737019778E-4,
184 6.98793669997260967291E-6,
185 8.11789239554389293311E-8,
186 3.41551784765923618484E-10,
187 };
188
189 static double AGD[10] = {
190 /* 1.00000000000000000000E0, */
191 9.30892908077441974853E0,
192 1.98352928718312140417E1,
193 1.55646628932864612953E1,
194 5.47686069422975497931E0,
195 9.54293611618961883998E-1,
196 8.64580826352392193095E-2,
197 4.12656523824222607191E-3,
198 1.01259085116509135510E-4,
199 1.17166733214413521882E-6,
200 4.91834570062930015649E-9,
201 };
202
203 static double APFN[9] = {
204 1.85365624022535566142E-1,
205 8.86712188052584095637E-1,
206 9.87391981747398547272E-1,
207 4.01241082318003734092E-1,
208 7.10304926289631174579E-2,
209 5.90618657995661810071E-3,
210 2.33051409401776799569E-4,
211 4.08718778289035454598E-6,
212 2.48379932900442457853E-8,
213 };
214
215 static double APFD[9] = {
216 /* 1.00000000000000000000E0, */
217 1.47345854687502542552E1,
218 3.75423933435489594466E1,
219 3.14657751203046424330E1,
220 1.09969125207298778536E1,
221 1.78885054766999417817E0,
222 1.41733275753662636873E-1,
223 5.44066067017226003627E-3,
224 9.39421290654511171663E-5,
225 5.65978713036027009243E-7,
226 };
227
228 static double APGN[11] = {
229 -3.55615429033082288335E-2,
230 -6.37311518129435504426E-1,
231 -1.70856738884312371053E0,
232 -1.50221872117316635393E0,
233 -5.63606665822102676611E-1,
234 -1.02101031120216891789E-1,
235 -9.48396695961445269093E-3,
236 -4.60325307486780994357E-4,
237 -1.14300836484517375919E-5,
238 -1.33415518685547420648E-7,
239 -5.63803833958893494476E-10,
240 };
241
242 static double APGD[11] = {
243 /* 1.00000000000000000000E0, */
244 9.85865801696130355144E0,
245 2.16401867356585941885E1,
246 1.73130776389749389525E1,
247 6.17872175280828766327E0,
248 1.08848694396321495475E0,
249 9.95005543440888479402E-2,
250 4.78468199683886610842E-3,
251 1.18159633322838625562E-4,
252 1.37480673554219441465E-6,
253 5.79912514929147598821E-9,
254 };
255
airy(x,ai,aip,bi,bip)256 int airy(x, ai, aip, bi, bip)
257 double x, *ai, *aip, *bi, *bip;
258 {
259 double z, zz, t, f, g, uf, ug, k, zeta, theta;
260 int domflg;
261
262 domflg = 0;
263 if (x > MAXAIRY) {
264 *ai = 0;
265 *aip = 0;
266 *bi = NPY_INFINITY;
267 *bip = NPY_INFINITY;
268 return (-1);
269 }
270
271 if (x < -2.09) {
272 domflg = 15;
273 t = sqrt(-x);
274 zeta = -2.0 * x * t / 3.0;
275 t = sqrt(t);
276 k = sqpii / t;
277 z = 1.0 / zeta;
278 zz = z * z;
279 uf = 1.0 + zz * polevl(zz, AFN, 8) / p1evl(zz, AFD, 9);
280 ug = z * polevl(zz, AGN, 10) / p1evl(zz, AGD, 10);
281 theta = zeta + 0.25 * NPY_PI;
282 f = sin(theta);
283 g = cos(theta);
284 *ai = k * (f * uf - g * ug);
285 *bi = k * (g * uf + f * ug);
286 uf = 1.0 + zz * polevl(zz, APFN, 8) / p1evl(zz, APFD, 9);
287 ug = z * polevl(zz, APGN, 10) / p1evl(zz, APGD, 10);
288 k = sqpii * t;
289 *aip = -k * (g * uf + f * ug);
290 *bip = k * (f * uf - g * ug);
291 return (0);
292 }
293
294 if (x >= 2.09) { /* cbrt(9) */
295 domflg = 5;
296 t = sqrt(x);
297 zeta = 2.0 * x * t / 3.0;
298 g = exp(zeta);
299 t = sqrt(t);
300 k = 2.0 * t * g;
301 z = 1.0 / zeta;
302 f = polevl(z, AN, 7) / polevl(z, AD, 7);
303 *ai = sqpii * f / k;
304 k = -0.5 * sqpii * t / g;
305 f = polevl(z, APN, 7) / polevl(z, APD, 7);
306 *aip = f * k;
307
308 if (x > 8.3203353) { /* zeta > 16 */
309 f = z * polevl(z, BN16, 4) / p1evl(z, BD16, 5);
310 k = sqpii * g;
311 *bi = k * (1.0 + f) / t;
312 f = z * polevl(z, BPPN, 4) / p1evl(z, BPPD, 5);
313 *bip = k * t * (1.0 + f);
314 return (0);
315 }
316 }
317
318 f = 1.0;
319 g = x;
320 t = 1.0;
321 uf = 1.0;
322 ug = x;
323 k = 1.0;
324 z = x * x * x;
325 while (t > MACHEP) {
326 uf *= z;
327 k += 1.0;
328 uf /= k;
329 ug *= z;
330 k += 1.0;
331 ug /= k;
332 uf /= k;
333 f += uf;
334 k += 1.0;
335 ug /= k;
336 g += ug;
337 t = fabs(uf / f);
338 }
339 uf = c1 * f;
340 ug = c2 * g;
341 if ((domflg & 1) == 0)
342 *ai = uf - ug;
343 if ((domflg & 2) == 0)
344 *bi = sqrt3 * (uf + ug);
345
346 /* the deriviative of ai */
347 k = 4.0;
348 uf = x * x / 2.0;
349 ug = z / 3.0;
350 f = uf;
351 g = 1.0 + ug;
352 uf /= 3.0;
353 t = 1.0;
354
355 while (t > MACHEP) {
356 uf *= z;
357 ug /= k;
358 k += 1.0;
359 ug *= z;
360 uf /= k;
361 f += uf;
362 k += 1.0;
363 ug /= k;
364 uf /= k;
365 g += ug;
366 k += 1.0;
367 t = fabs(ug / g);
368 }
369
370 uf = c1 * f;
371 ug = c2 * g;
372 if ((domflg & 4) == 0)
373 *aip = uf - ug;
374 if ((domflg & 8) == 0)
375 *bip = sqrt3 * (uf + ug);
376 return (0);
377 }
378