1//  Ported from Cephes by
2//  Andreas Kloeckner (C) 2012
3//
4// Cephes Math Library Release 2.8:  June, 2000
5// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
6// What you see here may be used freely, but it comes with no support or
7// guarantee.
8
9#pragma once
10
11#include <pyopencl-eval-tbl.cl>
12
13__constant const double airy_maxairy = 103.892;
14
15__constant const double airy_sqrt3 = 1.732050807568877293527;
16__constant const double airy_sqpii = 5.64189583547756286948E-1;
17
18__constant const double airy_c1 = 0.35502805388781723926;
19__constant const double airy_c2 = 0.258819403792806798405;
20
21__constant const unsigned short AN[32] = {
220x3fd6,0x2dae,0x2537,0xb658,
230x4028,0x03e3,0x871a,0x9067,
240x4053,0x11e5,0x0de2,0xe1e3,
250x4065,0x02da,0xee40,0x073c,
260x4063,0xf834,0x5ba1,0xfddf,
270x4051,0xa24f,0x4f4c,0xea4f,
280x402c,0x0d8d,0x5c2a,0x0f4d,
290x3ff0,0x0000,0x0000,0x0000,
30};
31__constant const unsigned short AD[32] = {
320x3fe2,0x29bc,0x0262,0x4d31,
330x402d,0x8334,0x0533,0x2ca5,
340x4055,0x20e3,0xb04d,0x51a0,
350x4066,0x2a2d,0xc730,0xb7b0,
360x4064,0x8782,0x9a9f,0xfa61,
370x4051,0xde94,0xee91,0xd35f,
380x402c,0x311b,0x950d,0x9d81,
390x3ff0,0x0000,0x0000,0x0000,
40};
41
42__constant const unsigned short APN[32] = {
430x3fe3,0xa3ea,0x4d4c,0xab3e,
440x402d,0x7dad,0xdc67,0x2bcf,
450x4054,0x83bd,0x0724,0xa9a6,
460x4065,0x65e9,0xba99,0xc9ba,
470x4063,0xea2b,0xcdc2,0x64d7,
480x4051,0x7e95,0x41d4,0x1646,
490x402b,0xe4e8,0x6aa7,0x4099,
500x3ff0,0x0000,0x0000,0x0000,
51};
52__constant const unsigned short APD[32] = {
530x3fd5,0x6397,0xd288,0xd5b3,
540x4026,0x5caf,0xedc9,0x327e,
550x4051,0xcb0e,0x1800,0x97e6,
560x4063,0xd8e6,0x1132,0xdbd1,
570x4063,0x269b,0x0dcb,0x3316,
580x4051,0x2b36,0xf9d0,0xf72f,
590x402b,0xb321,0x4e35,0x7982,
600x3ff0,0x0000,0x0000,0x0000,
61};
62
63__constant const unsigned short BN16[20] = {
640xbfd0,0x3518,0xe211,0x6751,
650x3fe2,0x68bc,0x7072,0x2383,
660xbfd5,0x1d32,0x6785,0xcf29,
670x3fb0,0x7f2a,0xa027,0x78a8,
680xbf6f,0x5604,0x2dba,0xcd1b,
69};
70__constant const unsigned short BD16[20] = {
71/*0x3ff0,0x0000,0x0000,0x0000,*/
720xc01c,0xa09d,0x891b,0xab58,
730x4025,0x3539,0xfe0b,0x1101,
740xc014,0xee0b,0xa9a7,0x70e8,
750x3fee,0xa2fc,0xa6da,0x95ff,
760xbfac,0x33d0,0x8f8e,0x86c9,
77};
78
79__constant const unsigned short BPPN[20] = {
800x3fdd,0xca1d,0x9deb,0x377b,
810xbff1,0x7051,0xc6be,0xe420,
820x3fe4,0x710c,0xf199,0x5ff3,
830xbfc0,0x3c6f,0x8681,0xa8fa,
840x3f7f,0x3b43,0xb8ce,0xb896,
85};
86__constant const unsigned short BPPD[20] = {
87/*0x3ff0,0x0000,0x0000,0x0000,*/
880xc021,0x6996,0xb340,0xbc45,
890x402b,0xcc73,0x2ea4,0xbb8b,
900xc01c,0x908c,0xa04a,0xed59,
910x3ff5,0x70fd,0xf9a5,0x70a9,
920xbfb4,0x13d0,0x1b60,0x52e8,
93};
94
95__constant const unsigned short AFN[36] = {
960xbfc0,0xdb6c,0xd50a,0xe6fb,
970xbfe4,0x0bee,0x9856,0x6852,
980xbfe6,0x2e59,0xc2f7,0x9f7d,
990xbfd1,0xe7ea,0x4bb3,0xf40b,
1000xbfa9,0x2f6e,0xf47d,0xbd8a,
1010xbf70,0xa401,0xc8d9,0xe090,
1020xbf24,0xe06e,0xaf4b,0x009c,
1030xbec7,0x4a78,0x1d42,0x366d,
1040xbe52,0x041c,0xf68e,0xa2d2,
105};
106__constant const unsigned short AFD[36] = {
107/*0x3ff0,0x0000,0x0000,0x0000,*/
1080x402a,0xb64b,0x2572,0xedf2,
1090x4040,0x575c,0x4478,0x7b1a,
1100x403a,0xbc98,0xa3b7,0x3410,
1110x4022,0x5fc8,0x2ac9,0x9873,
1120x3ff7,0x9acb,0x39de,0x9319,
1130x3fbd,0x9dac,0xb404,0x5a2b,
1140x3f72,0x08ca,0xe03a,0xf617,
1150x3f13,0xc8d7,0xaf76,0xe73b,
1160x3e9e,0x52b9,0xb995,0x18a7,
117};
118
119__constant const unsigned short AGN[44] = {
1200x3f94,0x3525,0xddcf,0xbbde,
1210x3fd9,0x07d5,0x0064,0x37b7,
1220x3ff1,0x0d83,0x3a20,0x34eb,
1230x3fee,0x0dac,0xa0ef,0x1acb,
1240x3fd6,0x7e69,0xcea8,0xfe1d,
1250x3fb0,0x3a41,0x21e9,0x0978,
1260x3f77,0xfe99,0xf12f,0x5043,
1270x3f32,0x8976,0x600e,0x17a2,
1280x3edd,0x4f3d,0x69f8,0x574e,
1290x3e75,0xca92,0xbbad,0x11c8,
1300x3df7,0x78a4,0x7d97,0xee7a,
131};
132__constant const unsigned short AGD[40] = {
133/*0x3ff0,0x0000,0x0000,0x0000,*/
1340x4022,0x9e2b,0xf3d5,0x6b40,
1350x4033,0xd5d5,0xc0ef,0x18d4,
1360x402f,0x211b,0x7ea7,0xdc35,
1370x4015,0xe84e,0x2b79,0xdbce,
1380x3fee,0x8992,0xc195,0xece3,
1390x3fb6,0x221d,0xed64,0xa9ee,
1400x3f70,0xe704,0x6be3,0x93bb,
1410x3f1a,0x8b61,0xd603,0xa5a0,
1420x3eb3,0xa845,0xdb07,0x24e8,
1430x3e35,0x1fc7,0x3dd5,0x89d4,
144};
145
146__constant const unsigned short APFN[36] = {
1470x3fc7,0xba0f,0x8e7d,0x5db5,
1480x3fec,0x5ff2,0x3d14,0xd07e,
1490x3fef,0x98b7,0x11be,0x01af,
1500x3fd9,0xadef,0x1397,0x84a1,
1510x3fb2,0x2f0d,0xeadc,0x33d1,
1520x3f78,0x3115,0xe347,0xa140,
1530x3f2e,0x8be8,0x5d03,0x8059,
1540x3ed1,0x2495,0x9f80,0x12af,
1550x3e5a,0xab6a,0x654d,0x7d86,
156};
157__constant const unsigned short APFD[36] = {
158/*0x3ff0,0x0000,0x0000,0x0000,*/
1590x402d,0x781b,0x9628,0xcc60,
1600x4042,0xc56d,0x2524,0x0e31,
1610x403f,0x773d,0x09cc,0xffb8,
1620x4025,0xfe6b,0x5163,0x03f7,
1630x3ffc,0x9f21,0xc07a,0xc9fd,
1640x3fc2,0x2450,0xe40e,0xf796,
1650x3f76,0x48f2,0x3a5a,0x351a,
1660x3f18,0xa059,0x7cfb,0x63a1,
1670x3ea2,0xfdb8,0x5a24,0x1e2e,
168};
169
170__constant const unsigned short APGN[44] = {
1710xbfa2,0x351f,0x5f87,0xaf5b,
1720xbfe4,0x64db,0x1ff7,0x5c76,
1730xbffb,0x564a,0xc221,0x7e49,
1740xbff8,0x0916,0x7f6e,0x0b07,
1750xbfe2,0x0910,0xd8b0,0x6edb,
1760xbfba,0x234b,0x0d8c,0x9903,
1770xbf83,0x6c54,0x7f6c,0x50df,
1780xbf3e,0x2afa,0x2424,0x2ad0,
1790xbee7,0xf87a,0xbc17,0xf631,
1800xbe81,0xe81f,0x501e,0x6c10,
1810xbe03,0x5f45,0x5e46,0x870d,
182};
183__constant const unsigned short APGD[40] = {
184/*0x3ff0,0x0000,0x0000,0x0000,*/
1850x4023,0xb7a2,0x060a,0x9812,
1860x4035,0xa3e3,0x4724,0xfc96,
1870x4031,0x5025,0xdb2c,0x819a,
1880x4018,0xb702,0xd5cd,0x94e2,
1890x3ff1,0x6a71,0x4927,0x1eb1,
1900x3fb9,0x78de,0x4ad7,0x7bc5,
1910x3f73,0x991a,0x4b2b,0xc1d7,
1920x3f1e,0xf98f,0x0b16,0xbe1c,
1930x3eb7,0x10bf,0xfdde,0x4ef3,
1940x3e38,0xe834,0x9dc8,0x647e,
195};
196
197int airy( double x, double *ai, double *aip, double *bi, double *bip )
198{
199  typedef __constant const double *data_t;
200  double z, zz, t, f, g, uf, ug, k, zeta, theta;
201  int domflg;
202
203  domflg = 0;
204  if( x > airy_maxairy )
205  {
206    *ai = 0;
207    *aip = 0;
208    *bi = DBL_MAX;
209    *bip = DBL_MAX;
210    return(-1);
211  }
212
213  if( x < -2.09 )
214  {
215    domflg = 15;
216    t = sqrt(-x);
217    zeta = -2.0 * x * t / 3.0;
218    t = sqrt(t);
219    k = airy_sqpii / t;
220    z = 1.0/zeta;
221    zz = z * z;
222    uf = 1.0 + zz * cephes_polevl( zz, (data_t) AFN, 8 ) / cephes_p1evl( zz, (data_t) AFD, 9 );
223    ug = z * cephes_polevl( zz, (data_t) AGN, 10 ) / cephes_p1evl( zz, (data_t) AGD, 10 );
224    theta = zeta + 0.25 * M_PI;
225    f = sin( theta );
226    g = cos( theta );
227    *ai = k * (f * uf - g * ug);
228    *bi = k * (g * uf + f * ug);
229    uf = 1.0 + zz * cephes_polevl( zz, (data_t) APFN, 8 ) / cephes_p1evl( zz, (data_t) APFD, 9 );
230    ug = z * cephes_polevl( zz, (data_t) APGN, 10 ) / cephes_p1evl( zz, (data_t) APGD, 10 );
231    k = airy_sqpii * t;
232    *aip = -k * (g * uf + f * ug);
233    *bip = k * (f * uf - g * ug);
234    return(0);
235  }
236
237  if( x >= 2.09 )       /* cbrt(9) */
238  {
239    domflg = 5;
240    t = sqrt(x);
241    zeta = 2.0 * x * t / 3.0;
242    g = exp( zeta );
243    t = sqrt(t);
244    k = 2.0 * t * g;
245    z = 1.0/zeta;
246    f = cephes_polevl( z, (data_t) AN, 7 ) / cephes_polevl( z, (data_t) AD, 7 );
247    *ai = airy_sqpii * f / k;
248    k = -0.5 * airy_sqpii * t / g;
249    f = cephes_polevl( z, (data_t) APN, 7 ) / cephes_polevl( z, (data_t) APD, 7 );
250    *aip = f * k;
251
252    if( x > 8.3203353 ) /* zeta > 16 */
253    {
254      f = z * cephes_polevl( z, (data_t) BN16, 4 ) / cephes_p1evl( z, (data_t) BD16, 5 );
255      k = airy_sqpii * g;
256      *bi = k * (1.0 + f) / t;
257      f = z * cephes_polevl( z, (data_t) BPPN, 4 ) / cephes_p1evl( z, (data_t) BPPD, 5 );
258      *bip = k * t * (1.0 + f);
259      return(0);
260    }
261  }
262
263  f = 1.0;
264  g = x;
265  t = 1.0;
266  uf = 1.0;
267  ug = x;
268  k = 1.0;
269  z = x * x * x;
270  while( t > DBL_EPSILON )
271  {
272    uf *= z;
273    k += 1.0;
274    uf /=k;
275    ug *= z;
276    k += 1.0;
277    ug /=k;
278    uf /=k;
279    f += uf;
280    k += 1.0;
281    ug /=k;
282    g += ug;
283    t = fabs(uf/f);
284  }
285  uf = airy_c1 * f;
286  ug = airy_c2 * g;
287  if( (domflg & 1) == 0 )
288    *ai = uf - ug;
289  if( (domflg & 2) == 0 )
290    *bi = airy_sqrt3 * (uf + ug);
291
292  /* the deriviative of ai */
293  k = 4.0;
294  uf = x * x/2.0;
295  ug = z/3.0;
296  f = uf;
297  g = 1.0 + ug;
298  uf /= 3.0;
299  t = 1.0;
300
301  while( t > DBL_EPSILON )
302  {
303    uf *= z;
304    ug /=k;
305    k += 1.0;
306    ug *= z;
307    uf /=k;
308    f += uf;
309    k += 1.0;
310    ug /=k;
311    uf /=k;
312    g += ug;
313    k += 1.0;
314    t = fabs(ug/g);
315  }
316
317  uf = airy_c1 * f;
318  ug = airy_c2 * g;
319  if( (domflg & 4) == 0 )
320    *aip = uf - ug;
321  if( (domflg & 8) == 0 )
322    *bip = airy_sqrt3 * (uf + ug);
323  return(0);
324}
325