airy.c 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965
  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. * DEC -10, 0 Ai 5000 1.7e-16 2.8e-17
  44. * DEC 0, 10 Ai 5000 2.1e-15* 1.7e-16*
  45. * DEC -10, 0 Ai' 5000 4.7e-16 7.8e-17
  46. * DEC 0, 10 Ai' 12000 1.8e-15* 1.5e-16*
  47. * DEC -10, 10 Bi 10000 5.5e-16 6.8e-17
  48. * DEC -10, 10 Bi' 7000 5.3e-16 8.7e-17
  49. *
  50. */
  51. /* airy.c */
  52. /*
  53. Cephes Math Library Release 2.8: June, 2000
  54. Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
  55. */
  56. #include <math.h>
  57. static double c1 = 0.35502805388781723926;
  58. static double c2 = 0.258819403792806798405;
  59. static double sqrt3 = 1.732050807568877293527;
  60. static double sqpii = 5.64189583547756286948E-1;
  61. extern double PI;
  62. extern double MAXNUM, MACHEP;
  63. #ifdef UNK
  64. #define MAXAIRY 25.77
  65. #endif
  66. #ifdef DEC
  67. #define MAXAIRY 25.77
  68. #endif
  69. #ifdef IBMPC
  70. #define MAXAIRY 103.892
  71. #endif
  72. #ifdef MIEEE
  73. #define MAXAIRY 103.892
  74. #endif
  75. #ifdef UNK
  76. static double AN[8] = {
  77. 3.46538101525629032477E-1,
  78. 1.20075952739645805542E1,
  79. 7.62796053615234516538E1,
  80. 1.68089224934630576269E2,
  81. 1.59756391350164413639E2,
  82. 7.05360906840444183113E1,
  83. 1.40264691163389668864E1,
  84. 9.99999999999999995305E-1,
  85. };
  86. static double AD[8] = {
  87. 5.67594532638770212846E-1,
  88. 1.47562562584847203173E1,
  89. 8.45138970141474626562E1,
  90. 1.77318088145400459522E2,
  91. 1.64234692871529701831E2,
  92. 7.14778400825575695274E1,
  93. 1.40959135607834029598E1,
  94. 1.00000000000000000470E0,
  95. };
  96. #endif
  97. #ifdef DEC
  98. static unsigned short AN[32] = {
  99. 0037661,0066561,0024675,0131301,
  100. 0041100,0017434,0034324,0101466,
  101. 0041630,0107450,0067427,0007430,
  102. 0042050,0013327,0071000,0034737,
  103. 0042037,0140642,0156417,0167366,
  104. 0041615,0011172,0075147,0051165,
  105. 0041140,0066152,0160520,0075146,
  106. 0040200,0000000,0000000,0000000,
  107. };
  108. static unsigned short AD[32] = {
  109. 0040021,0046740,0011422,0064606,
  110. 0041154,0014640,0024631,0062450,
  111. 0041651,0003435,0101152,0106401,
  112. 0042061,0050556,0034605,0136602,
  113. 0042044,0036024,0152377,0151414,
  114. 0041616,0172247,0072216,0115374,
  115. 0041141,0104334,0124154,0166007,
  116. 0040200,0000000,0000000,0000000,
  117. };
  118. #endif
  119. #ifdef IBMPC
  120. static unsigned short AN[32] = {
  121. 0xb658,0x2537,0x2dae,0x3fd6,
  122. 0x9067,0x871a,0x03e3,0x4028,
  123. 0xe1e3,0x0de2,0x11e5,0x4053,
  124. 0x073c,0xee40,0x02da,0x4065,
  125. 0xfddf,0x5ba1,0xf834,0x4063,
  126. 0xea4f,0x4f4c,0xa24f,0x4051,
  127. 0x0f4d,0x5c2a,0x0d8d,0x402c,
  128. 0x0000,0x0000,0x0000,0x3ff0,
  129. };
  130. static unsigned short AD[32] = {
  131. 0x4d31,0x0262,0x29bc,0x3fe2,
  132. 0x2ca5,0x0533,0x8334,0x402d,
  133. 0x51a0,0xb04d,0x20e3,0x4055,
  134. 0xb7b0,0xc730,0x2a2d,0x4066,
  135. 0xfa61,0x9a9f,0x8782,0x4064,
  136. 0xd35f,0xee91,0xde94,0x4051,
  137. 0x9d81,0x950d,0x311b,0x402c,
  138. 0x0000,0x0000,0x0000,0x3ff0,
  139. };
  140. #endif
  141. #ifdef MIEEE
  142. static unsigned short AN[32] = {
  143. 0x3fd6,0x2dae,0x2537,0xb658,
  144. 0x4028,0x03e3,0x871a,0x9067,
  145. 0x4053,0x11e5,0x0de2,0xe1e3,
  146. 0x4065,0x02da,0xee40,0x073c,
  147. 0x4063,0xf834,0x5ba1,0xfddf,
  148. 0x4051,0xa24f,0x4f4c,0xea4f,
  149. 0x402c,0x0d8d,0x5c2a,0x0f4d,
  150. 0x3ff0,0x0000,0x0000,0x0000,
  151. };
  152. static unsigned short AD[32] = {
  153. 0x3fe2,0x29bc,0x0262,0x4d31,
  154. 0x402d,0x8334,0x0533,0x2ca5,
  155. 0x4055,0x20e3,0xb04d,0x51a0,
  156. 0x4066,0x2a2d,0xc730,0xb7b0,
  157. 0x4064,0x8782,0x9a9f,0xfa61,
  158. 0x4051,0xde94,0xee91,0xd35f,
  159. 0x402c,0x311b,0x950d,0x9d81,
  160. 0x3ff0,0x0000,0x0000,0x0000,
  161. };
  162. #endif
  163. #ifdef UNK
  164. static double APN[8] = {
  165. 6.13759184814035759225E-1,
  166. 1.47454670787755323881E1,
  167. 8.20584123476060982430E1,
  168. 1.71184781360976385540E2,
  169. 1.59317847137141783523E2,
  170. 6.99778599330103016170E1,
  171. 1.39470856980481566958E1,
  172. 1.00000000000000000550E0,
  173. };
  174. static double APD[8] = {
  175. 3.34203677749736953049E-1,
  176. 1.11810297306158156705E1,
  177. 7.11727352147859965283E1,
  178. 1.58778084372838313640E2,
  179. 1.53206427475809220834E2,
  180. 6.86752304592780337944E1,
  181. 1.38498634758259442477E1,
  182. 9.99999999999999994502E-1,
  183. };
  184. #endif
  185. #ifdef DEC
  186. static unsigned short APN[32] = {
  187. 0040035,0017522,0065145,0054755,
  188. 0041153,0166556,0161471,0057174,
  189. 0041644,0016750,0034445,0046462,
  190. 0042053,0027515,0152316,0046717,
  191. 0042037,0050536,0067023,0023264,
  192. 0041613,0172252,0007240,0131055,
  193. 0041137,0023503,0052472,0002305,
  194. 0040200,0000000,0000000,0000000,
  195. };
  196. static unsigned short APD[32] = {
  197. 0037653,0016276,0112106,0126625,
  198. 0041062,0162577,0067111,0111761,
  199. 0041616,0054160,0140004,0137455,
  200. 0042036,0143460,0104626,0157206,
  201. 0042031,0032330,0067131,0114260,
  202. 0041611,0054667,0147207,0134564,
  203. 0041135,0114412,0070653,0146015,
  204. 0040200,0000000,0000000,0000000,
  205. };
  206. #endif
  207. #ifdef IBMPC
  208. static unsigned short APN[32] = {
  209. 0xab3e,0x4d4c,0xa3ea,0x3fe3,
  210. 0x2bcf,0xdc67,0x7dad,0x402d,
  211. 0xa9a6,0x0724,0x83bd,0x4054,
  212. 0xc9ba,0xba99,0x65e9,0x4065,
  213. 0x64d7,0xcdc2,0xea2b,0x4063,
  214. 0x1646,0x41d4,0x7e95,0x4051,
  215. 0x4099,0x6aa7,0xe4e8,0x402b,
  216. 0x0000,0x0000,0x0000,0x3ff0,
  217. };
  218. static unsigned short APD[32] = {
  219. 0xd5b3,0xd288,0x6397,0x3fd5,
  220. 0x327e,0xedc9,0x5caf,0x4026,
  221. 0x97e6,0x1800,0xcb0e,0x4051,
  222. 0xdbd1,0x1132,0xd8e6,0x4063,
  223. 0x3316,0x0dcb,0x269b,0x4063,
  224. 0xf72f,0xf9d0,0x2b36,0x4051,
  225. 0x7982,0x4e35,0xb321,0x402b,
  226. 0x0000,0x0000,0x0000,0x3ff0,
  227. };
  228. #endif
  229. #ifdef MIEEE
  230. static unsigned short APN[32] = {
  231. 0x3fe3,0xa3ea,0x4d4c,0xab3e,
  232. 0x402d,0x7dad,0xdc67,0x2bcf,
  233. 0x4054,0x83bd,0x0724,0xa9a6,
  234. 0x4065,0x65e9,0xba99,0xc9ba,
  235. 0x4063,0xea2b,0xcdc2,0x64d7,
  236. 0x4051,0x7e95,0x41d4,0x1646,
  237. 0x402b,0xe4e8,0x6aa7,0x4099,
  238. 0x3ff0,0x0000,0x0000,0x0000,
  239. };
  240. static unsigned short APD[32] = {
  241. 0x3fd5,0x6397,0xd288,0xd5b3,
  242. 0x4026,0x5caf,0xedc9,0x327e,
  243. 0x4051,0xcb0e,0x1800,0x97e6,
  244. 0x4063,0xd8e6,0x1132,0xdbd1,
  245. 0x4063,0x269b,0x0dcb,0x3316,
  246. 0x4051,0x2b36,0xf9d0,0xf72f,
  247. 0x402b,0xb321,0x4e35,0x7982,
  248. 0x3ff0,0x0000,0x0000,0x0000,
  249. };
  250. #endif
  251. #ifdef UNK
  252. static double BN16[5] = {
  253. -2.53240795869364152689E-1,
  254. 5.75285167332467384228E-1,
  255. -3.29907036873225371650E-1,
  256. 6.44404068948199951727E-2,
  257. -3.82519546641336734394E-3,
  258. };
  259. static double BD16[5] = {
  260. /* 1.00000000000000000000E0,*/
  261. -7.15685095054035237902E0,
  262. 1.06039580715664694291E1,
  263. -5.23246636471251500874E0,
  264. 9.57395864378383833152E-1,
  265. -5.50828147163549611107E-2,
  266. };
  267. #endif
  268. #ifdef DEC
  269. static unsigned short BN16[20] = {
  270. 0137601,0124307,0010213,0035210,
  271. 0040023,0042743,0101621,0016031,
  272. 0137650,0164623,0036056,0074511,
  273. 0037203,0174525,0000473,0142474,
  274. 0136172,0130041,0066726,0064324,
  275. };
  276. static unsigned short BD16[20] = {
  277. /*0040200,0000000,0000000,0000000,*/
  278. 0140745,0002354,0044335,0055276,
  279. 0041051,0124717,0170130,0104013,
  280. 0140647,0070135,0046473,0103501,
  281. 0040165,0013745,0033324,0127766,
  282. 0137141,0117204,0076164,0033107,
  283. };
  284. #endif
  285. #ifdef IBMPC
  286. static unsigned short BN16[20] = {
  287. 0x6751,0xe211,0x3518,0xbfd0,
  288. 0x2383,0x7072,0x68bc,0x3fe2,
  289. 0xcf29,0x6785,0x1d32,0xbfd5,
  290. 0x78a8,0xa027,0x7f2a,0x3fb0,
  291. 0xcd1b,0x2dba,0x5604,0xbf6f,
  292. };
  293. static unsigned short BD16[20] = {
  294. /*0x0000,0x0000,0x0000,0x3ff0,*/
  295. 0xab58,0x891b,0xa09d,0xc01c,
  296. 0x1101,0xfe0b,0x3539,0x4025,
  297. 0x70e8,0xa9a7,0xee0b,0xc014,
  298. 0x95ff,0xa6da,0xa2fc,0x3fee,
  299. 0x86c9,0x8f8e,0x33d0,0xbfac,
  300. };
  301. #endif
  302. #ifdef MIEEE
  303. static unsigned short BN16[20] = {
  304. 0xbfd0,0x3518,0xe211,0x6751,
  305. 0x3fe2,0x68bc,0x7072,0x2383,
  306. 0xbfd5,0x1d32,0x6785,0xcf29,
  307. 0x3fb0,0x7f2a,0xa027,0x78a8,
  308. 0xbf6f,0x5604,0x2dba,0xcd1b,
  309. };
  310. static unsigned short BD16[20] = {
  311. /*0x3ff0,0x0000,0x0000,0x0000,*/
  312. 0xc01c,0xa09d,0x891b,0xab58,
  313. 0x4025,0x3539,0xfe0b,0x1101,
  314. 0xc014,0xee0b,0xa9a7,0x70e8,
  315. 0x3fee,0xa2fc,0xa6da,0x95ff,
  316. 0xbfac,0x33d0,0x8f8e,0x86c9,
  317. };
  318. #endif
  319. #ifdef UNK
  320. static double BPPN[5] = {
  321. 4.65461162774651610328E-1,
  322. -1.08992173800493920734E0,
  323. 6.38800117371827987759E-1,
  324. -1.26844349553102907034E-1,
  325. 7.62487844342109852105E-3,
  326. };
  327. static double BPPD[5] = {
  328. /* 1.00000000000000000000E0,*/
  329. -8.70622787633159124240E0,
  330. 1.38993162704553213172E1,
  331. -7.14116144616431159572E0,
  332. 1.34008595960680518666E0,
  333. -7.84273211323341930448E-2,
  334. };
  335. #endif
  336. #ifdef DEC
  337. static unsigned short BPPN[20] = {
  338. 0037756,0050354,0167531,0135731,
  339. 0140213,0101216,0032767,0020375,
  340. 0040043,0104147,0106312,0177632,
  341. 0137401,0161574,0032015,0043714,
  342. 0036371,0155035,0143165,0142262,
  343. };
  344. static unsigned short BPPD[20] = {
  345. /*0040200,0000000,0000000,0000000,*/
  346. 0141013,0046265,0115005,0161053,
  347. 0041136,0061631,0072445,0156131,
  348. 0140744,0102145,0001127,0065304,
  349. 0040253,0103757,0146453,0102513,
  350. 0137240,0117200,0155402,0113500,
  351. };
  352. #endif
  353. #ifdef IBMPC
  354. static unsigned short BPPN[20] = {
  355. 0x377b,0x9deb,0xca1d,0x3fdd,
  356. 0xe420,0xc6be,0x7051,0xbff1,
  357. 0x5ff3,0xf199,0x710c,0x3fe4,
  358. 0xa8fa,0x8681,0x3c6f,0xbfc0,
  359. 0xb896,0xb8ce,0x3b43,0x3f7f,
  360. };
  361. static unsigned short BPPD[20] = {
  362. /*0x0000,0x0000,0x0000,0x3ff0,*/
  363. 0xbc45,0xb340,0x6996,0xc021,
  364. 0xbb8b,0x2ea4,0xcc73,0x402b,
  365. 0xed59,0xa04a,0x908c,0xc01c,
  366. 0x70a9,0xf9a5,0x70fd,0x3ff5,
  367. 0x52e8,0x1b60,0x13d0,0xbfb4,
  368. };
  369. #endif
  370. #ifdef MIEEE
  371. static unsigned short BPPN[20] = {
  372. 0x3fdd,0xca1d,0x9deb,0x377b,
  373. 0xbff1,0x7051,0xc6be,0xe420,
  374. 0x3fe4,0x710c,0xf199,0x5ff3,
  375. 0xbfc0,0x3c6f,0x8681,0xa8fa,
  376. 0x3f7f,0x3b43,0xb8ce,0xb896,
  377. };
  378. static unsigned short BPPD[20] = {
  379. /*0x3ff0,0x0000,0x0000,0x0000,*/
  380. 0xc021,0x6996,0xb340,0xbc45,
  381. 0x402b,0xcc73,0x2ea4,0xbb8b,
  382. 0xc01c,0x908c,0xa04a,0xed59,
  383. 0x3ff5,0x70fd,0xf9a5,0x70a9,
  384. 0xbfb4,0x13d0,0x1b60,0x52e8,
  385. };
  386. #endif
  387. #ifdef UNK
  388. static double AFN[9] = {
  389. -1.31696323418331795333E-1,
  390. -6.26456544431912369773E-1,
  391. -6.93158036036933542233E-1,
  392. -2.79779981545119124951E-1,
  393. -4.91900132609500318020E-2,
  394. -4.06265923594885404393E-3,
  395. -1.59276496239262096340E-4,
  396. -2.77649108155232920844E-6,
  397. -1.67787698489114633780E-8,
  398. };
  399. static double AFD[9] = {
  400. /* 1.00000000000000000000E0,*/
  401. 1.33560420706553243746E1,
  402. 3.26825032795224613948E1,
  403. 2.67367040941499554804E1,
  404. 9.18707402907259625840E0,
  405. 1.47529146771666414581E0,
  406. 1.15687173795188044134E-1,
  407. 4.40291641615211203805E-3,
  408. 7.54720348287414296618E-5,
  409. 4.51850092970580378464E-7,
  410. };
  411. #endif
  412. #ifdef DEC
  413. static unsigned short AFN[36] = {
  414. 0137406,0155546,0124127,0033732,
  415. 0140040,0057564,0141263,0041222,
  416. 0140061,0071316,0013674,0175754,
  417. 0137617,0037522,0056637,0120130,
  418. 0137111,0075567,0121755,0166122,
  419. 0136205,0020016,0043317,0002201,
  420. 0135047,0001565,0075130,0002334,
  421. 0133472,0051700,0165021,0131551,
  422. 0131620,0020347,0132165,0013215,
  423. };
  424. static unsigned short AFD[36] = {
  425. /*0040200,0000000,0000000,0000000,*/
  426. 0041125,0131131,0025627,0067623,
  427. 0041402,0135342,0021703,0154315,
  428. 0041325,0162305,0016671,0120175,
  429. 0041022,0177101,0053114,0141632,
  430. 0040274,0153131,0147364,0114306,
  431. 0037354,0166545,0120042,0150530,
  432. 0036220,0043127,0000727,0130273,
  433. 0034636,0043275,0075667,0034733,
  434. 0032762,0112715,0146250,0142474,
  435. };
  436. #endif
  437. #ifdef IBMPC
  438. static unsigned short AFN[36] = {
  439. 0xe6fb,0xd50a,0xdb6c,0xbfc0,
  440. 0x6852,0x9856,0x0bee,0xbfe4,
  441. 0x9f7d,0xc2f7,0x2e59,0xbfe6,
  442. 0xf40b,0x4bb3,0xe7ea,0xbfd1,
  443. 0xbd8a,0xf47d,0x2f6e,0xbfa9,
  444. 0xe090,0xc8d9,0xa401,0xbf70,
  445. 0x009c,0xaf4b,0xe06e,0xbf24,
  446. 0x366d,0x1d42,0x4a78,0xbec7,
  447. 0xa2d2,0xf68e,0x041c,0xbe52,
  448. };
  449. static unsigned short AFD[36] = {
  450. /*0x0000,0x0000,0x0000,0x3ff0,*/
  451. 0xedf2,0x2572,0xb64b,0x402a,
  452. 0x7b1a,0x4478,0x575c,0x4040,
  453. 0x3410,0xa3b7,0xbc98,0x403a,
  454. 0x9873,0x2ac9,0x5fc8,0x4022,
  455. 0x9319,0x39de,0x9acb,0x3ff7,
  456. 0x5a2b,0xb404,0x9dac,0x3fbd,
  457. 0xf617,0xe03a,0x08ca,0x3f72,
  458. 0xe73b,0xaf76,0xc8d7,0x3f13,
  459. 0x18a7,0xb995,0x52b9,0x3e9e,
  460. };
  461. #endif
  462. #ifdef MIEEE
  463. static unsigned short AFN[36] = {
  464. 0xbfc0,0xdb6c,0xd50a,0xe6fb,
  465. 0xbfe4,0x0bee,0x9856,0x6852,
  466. 0xbfe6,0x2e59,0xc2f7,0x9f7d,
  467. 0xbfd1,0xe7ea,0x4bb3,0xf40b,
  468. 0xbfa9,0x2f6e,0xf47d,0xbd8a,
  469. 0xbf70,0xa401,0xc8d9,0xe090,
  470. 0xbf24,0xe06e,0xaf4b,0x009c,
  471. 0xbec7,0x4a78,0x1d42,0x366d,
  472. 0xbe52,0x041c,0xf68e,0xa2d2,
  473. };
  474. static unsigned short AFD[36] = {
  475. /*0x3ff0,0x0000,0x0000,0x0000,*/
  476. 0x402a,0xb64b,0x2572,0xedf2,
  477. 0x4040,0x575c,0x4478,0x7b1a,
  478. 0x403a,0xbc98,0xa3b7,0x3410,
  479. 0x4022,0x5fc8,0x2ac9,0x9873,
  480. 0x3ff7,0x9acb,0x39de,0x9319,
  481. 0x3fbd,0x9dac,0xb404,0x5a2b,
  482. 0x3f72,0x08ca,0xe03a,0xf617,
  483. 0x3f13,0xc8d7,0xaf76,0xe73b,
  484. 0x3e9e,0x52b9,0xb995,0x18a7,
  485. };
  486. #endif
  487. #ifdef UNK
  488. static double AGN[11] = {
  489. 1.97339932091685679179E-2,
  490. 3.91103029615688277255E-1,
  491. 1.06579897599595591108E0,
  492. 9.39169229816650230044E-1,
  493. 3.51465656105547619242E-1,
  494. 6.33888919628925490927E-2,
  495. 5.85804113048388458567E-3,
  496. 2.82851600836737019778E-4,
  497. 6.98793669997260967291E-6,
  498. 8.11789239554389293311E-8,
  499. 3.41551784765923618484E-10,
  500. };
  501. static double AGD[10] = {
  502. /* 1.00000000000000000000E0,*/
  503. 9.30892908077441974853E0,
  504. 1.98352928718312140417E1,
  505. 1.55646628932864612953E1,
  506. 5.47686069422975497931E0,
  507. 9.54293611618961883998E-1,
  508. 8.64580826352392193095E-2,
  509. 4.12656523824222607191E-3,
  510. 1.01259085116509135510E-4,
  511. 1.17166733214413521882E-6,
  512. 4.91834570062930015649E-9,
  513. };
  514. #endif
  515. #ifdef DEC
  516. static unsigned short AGN[44] = {
  517. 0036641,0124456,0167175,0157354,
  518. 0037710,0037250,0001441,0136671,
  519. 0040210,0066031,0150401,0123532,
  520. 0040160,0066545,0003570,0153133,
  521. 0037663,0171516,0072507,0170345,
  522. 0037201,0151011,0007510,0045702,
  523. 0036277,0172317,0104572,0101030,
  524. 0035224,0045663,0000160,0136422,
  525. 0033752,0074753,0047702,0135160,
  526. 0032256,0052225,0156550,0107103,
  527. 0030273,0142443,0166277,0071720,
  528. };
  529. static unsigned short AGD[40] = {
  530. /*0040200,0000000,0000000,0000000,*/
  531. 0041024,0170537,0117253,0055003,
  532. 0041236,0127256,0003570,0143240,
  533. 0041171,0004333,0172476,0160645,
  534. 0040657,0041161,0055716,0157161,
  535. 0040164,0046226,0006257,0063431,
  536. 0037261,0010357,0065445,0047563,
  537. 0036207,0034043,0057434,0116732,
  538. 0034724,0055416,0130035,0026377,
  539. 0033235,0041056,0154071,0023502,
  540. 0031250,0177071,0167254,0047242,
  541. };
  542. #endif
  543. #ifdef IBMPC
  544. static unsigned short AGN[44] = {
  545. 0xbbde,0xddcf,0x3525,0x3f94,
  546. 0x37b7,0x0064,0x07d5,0x3fd9,
  547. 0x34eb,0x3a20,0x0d83,0x3ff1,
  548. 0x1acb,0xa0ef,0x0dac,0x3fee,
  549. 0xfe1d,0xcea8,0x7e69,0x3fd6,
  550. 0x0978,0x21e9,0x3a41,0x3fb0,
  551. 0x5043,0xf12f,0xfe99,0x3f77,
  552. 0x17a2,0x600e,0x8976,0x3f32,
  553. 0x574e,0x69f8,0x4f3d,0x3edd,
  554. 0x11c8,0xbbad,0xca92,0x3e75,
  555. 0xee7a,0x7d97,0x78a4,0x3df7,
  556. };
  557. static unsigned short AGD[40] = {
  558. /*0x0000,0x0000,0x0000,0x3ff0,*/
  559. 0x6b40,0xf3d5,0x9e2b,0x4022,
  560. 0x18d4,0xc0ef,0xd5d5,0x4033,
  561. 0xdc35,0x7ea7,0x211b,0x402f,
  562. 0xdbce,0x2b79,0xe84e,0x4015,
  563. 0xece3,0xc195,0x8992,0x3fee,
  564. 0xa9ee,0xed64,0x221d,0x3fb6,
  565. 0x93bb,0x6be3,0xe704,0x3f70,
  566. 0xa5a0,0xd603,0x8b61,0x3f1a,
  567. 0x24e8,0xdb07,0xa845,0x3eb3,
  568. 0x89d4,0x3dd5,0x1fc7,0x3e35,
  569. };
  570. #endif
  571. #ifdef MIEEE
  572. static unsigned short AGN[44] = {
  573. 0x3f94,0x3525,0xddcf,0xbbde,
  574. 0x3fd9,0x07d5,0x0064,0x37b7,
  575. 0x3ff1,0x0d83,0x3a20,0x34eb,
  576. 0x3fee,0x0dac,0xa0ef,0x1acb,
  577. 0x3fd6,0x7e69,0xcea8,0xfe1d,
  578. 0x3fb0,0x3a41,0x21e9,0x0978,
  579. 0x3f77,0xfe99,0xf12f,0x5043,
  580. 0x3f32,0x8976,0x600e,0x17a2,
  581. 0x3edd,0x4f3d,0x69f8,0x574e,
  582. 0x3e75,0xca92,0xbbad,0x11c8,
  583. 0x3df7,0x78a4,0x7d97,0xee7a,
  584. };
  585. static unsigned short AGD[40] = {
  586. /*0x3ff0,0x0000,0x0000,0x0000,*/
  587. 0x4022,0x9e2b,0xf3d5,0x6b40,
  588. 0x4033,0xd5d5,0xc0ef,0x18d4,
  589. 0x402f,0x211b,0x7ea7,0xdc35,
  590. 0x4015,0xe84e,0x2b79,0xdbce,
  591. 0x3fee,0x8992,0xc195,0xece3,
  592. 0x3fb6,0x221d,0xed64,0xa9ee,
  593. 0x3f70,0xe704,0x6be3,0x93bb,
  594. 0x3f1a,0x8b61,0xd603,0xa5a0,
  595. 0x3eb3,0xa845,0xdb07,0x24e8,
  596. 0x3e35,0x1fc7,0x3dd5,0x89d4,
  597. };
  598. #endif
  599. #ifdef UNK
  600. static double APFN[9] = {
  601. 1.85365624022535566142E-1,
  602. 8.86712188052584095637E-1,
  603. 9.87391981747398547272E-1,
  604. 4.01241082318003734092E-1,
  605. 7.10304926289631174579E-2,
  606. 5.90618657995661810071E-3,
  607. 2.33051409401776799569E-4,
  608. 4.08718778289035454598E-6,
  609. 2.48379932900442457853E-8,
  610. };
  611. static double APFD[9] = {
  612. /* 1.00000000000000000000E0,*/
  613. 1.47345854687502542552E1,
  614. 3.75423933435489594466E1,
  615. 3.14657751203046424330E1,
  616. 1.09969125207298778536E1,
  617. 1.78885054766999417817E0,
  618. 1.41733275753662636873E-1,
  619. 5.44066067017226003627E-3,
  620. 9.39421290654511171663E-5,
  621. 5.65978713036027009243E-7,
  622. };
  623. #endif
  624. #ifdef DEC
  625. static unsigned short APFN[36] = {
  626. 0037475,0150174,0071752,0166651,
  627. 0040142,0177621,0164246,0101757,
  628. 0040174,0142670,0106760,0006573,
  629. 0037715,0067570,0116274,0022404,
  630. 0037221,0074157,0053341,0117207,
  631. 0036301,0104257,0015075,0004777,
  632. 0035164,0057502,0164034,0001313,
  633. 0033611,0022254,0176000,0112565,
  634. 0031725,0055523,0025153,0166057,
  635. };
  636. static unsigned short APFD[36] = {
  637. /*0040200,0000000,0000000,0000000,*/
  638. 0041153,0140334,0130506,0061402,
  639. 0041426,0025551,0024440,0070611,
  640. 0041373,0134750,0047147,0176702,
  641. 0041057,0171532,0105430,0017674,
  642. 0040344,0174416,0001726,0047754,
  643. 0037421,0021207,0020167,0136264,
  644. 0036262,0043621,0151321,0124324,
  645. 0034705,0001313,0163733,0016407,
  646. 0033027,0166702,0150440,0170561,
  647. };
  648. #endif
  649. #ifdef IBMPC
  650. static unsigned short APFN[36] = {
  651. 0x5db5,0x8e7d,0xba0f,0x3fc7,
  652. 0xd07e,0x3d14,0x5ff2,0x3fec,
  653. 0x01af,0x11be,0x98b7,0x3fef,
  654. 0x84a1,0x1397,0xadef,0x3fd9,
  655. 0x33d1,0xeadc,0x2f0d,0x3fb2,
  656. 0xa140,0xe347,0x3115,0x3f78,
  657. 0x8059,0x5d03,0x8be8,0x3f2e,
  658. 0x12af,0x9f80,0x2495,0x3ed1,
  659. 0x7d86,0x654d,0xab6a,0x3e5a,
  660. };
  661. static unsigned short APFD[36] = {
  662. /*0x0000,0x0000,0x0000,0x3ff0,*/
  663. 0xcc60,0x9628,0x781b,0x402d,
  664. 0x0e31,0x2524,0xc56d,0x4042,
  665. 0xffb8,0x09cc,0x773d,0x403f,
  666. 0x03f7,0x5163,0xfe6b,0x4025,
  667. 0xc9fd,0xc07a,0x9f21,0x3ffc,
  668. 0xf796,0xe40e,0x2450,0x3fc2,
  669. 0x351a,0x3a5a,0x48f2,0x3f76,
  670. 0x63a1,0x7cfb,0xa059,0x3f18,
  671. 0x1e2e,0x5a24,0xfdb8,0x3ea2,
  672. };
  673. #endif
  674. #ifdef MIEEE
  675. static unsigned short APFN[36] = {
  676. 0x3fc7,0xba0f,0x8e7d,0x5db5,
  677. 0x3fec,0x5ff2,0x3d14,0xd07e,
  678. 0x3fef,0x98b7,0x11be,0x01af,
  679. 0x3fd9,0xadef,0x1397,0x84a1,
  680. 0x3fb2,0x2f0d,0xeadc,0x33d1,
  681. 0x3f78,0x3115,0xe347,0xa140,
  682. 0x3f2e,0x8be8,0x5d03,0x8059,
  683. 0x3ed1,0x2495,0x9f80,0x12af,
  684. 0x3e5a,0xab6a,0x654d,0x7d86,
  685. };
  686. static unsigned short APFD[36] = {
  687. /*0x3ff0,0x0000,0x0000,0x0000,*/
  688. 0x402d,0x781b,0x9628,0xcc60,
  689. 0x4042,0xc56d,0x2524,0x0e31,
  690. 0x403f,0x773d,0x09cc,0xffb8,
  691. 0x4025,0xfe6b,0x5163,0x03f7,
  692. 0x3ffc,0x9f21,0xc07a,0xc9fd,
  693. 0x3fc2,0x2450,0xe40e,0xf796,
  694. 0x3f76,0x48f2,0x3a5a,0x351a,
  695. 0x3f18,0xa059,0x7cfb,0x63a1,
  696. 0x3ea2,0xfdb8,0x5a24,0x1e2e,
  697. };
  698. #endif
  699. #ifdef UNK
  700. static double APGN[11] = {
  701. -3.55615429033082288335E-2,
  702. -6.37311518129435504426E-1,
  703. -1.70856738884312371053E0,
  704. -1.50221872117316635393E0,
  705. -5.63606665822102676611E-1,
  706. -1.02101031120216891789E-1,
  707. -9.48396695961445269093E-3,
  708. -4.60325307486780994357E-4,
  709. -1.14300836484517375919E-5,
  710. -1.33415518685547420648E-7,
  711. -5.63803833958893494476E-10,
  712. };
  713. static double APGD[11] = {
  714. /* 1.00000000000000000000E0,*/
  715. 9.85865801696130355144E0,
  716. 2.16401867356585941885E1,
  717. 1.73130776389749389525E1,
  718. 6.17872175280828766327E0,
  719. 1.08848694396321495475E0,
  720. 9.95005543440888479402E-2,
  721. 4.78468199683886610842E-3,
  722. 1.18159633322838625562E-4,
  723. 1.37480673554219441465E-6,
  724. 5.79912514929147598821E-9,
  725. };
  726. #endif
  727. #ifdef DEC
  728. static unsigned short APGN[44] = {
  729. 0137021,0124372,0176075,0075331,
  730. 0140043,0023330,0177672,0161655,
  731. 0140332,0131126,0010413,0171112,
  732. 0140300,0044263,0175560,0054070,
  733. 0140020,0044206,0142603,0073324,
  734. 0137321,0015130,0066144,0144033,
  735. 0136433,0061243,0175542,0103373,
  736. 0135361,0053721,0020441,0053203,
  737. 0134077,0141725,0160277,0130612,
  738. 0132417,0040372,0100363,0060200,
  739. 0130432,0175052,0171064,0034147,
  740. };
  741. static unsigned short APGD[40] = {
  742. /*0040200,0000000,0000000,0000000,*/
  743. 0041035,0136420,0030124,0140220,
  744. 0041255,0017432,0034447,0162256,
  745. 0041212,0100456,0154544,0006321,
  746. 0040705,0134026,0127154,0123414,
  747. 0040213,0051612,0044470,0172607,
  748. 0037313,0143362,0053273,0157051,
  749. 0036234,0144322,0054536,0007264,
  750. 0034767,0146170,0054265,0170342,
  751. 0033270,0102777,0167362,0073631,
  752. 0031307,0040644,0167103,0021763,
  753. };
  754. #endif
  755. #ifdef IBMPC
  756. static unsigned short APGN[44] = {
  757. 0xaf5b,0x5f87,0x351f,0xbfa2,
  758. 0x5c76,0x1ff7,0x64db,0xbfe4,
  759. 0x7e49,0xc221,0x564a,0xbffb,
  760. 0x0b07,0x7f6e,0x0916,0xbff8,
  761. 0x6edb,0xd8b0,0x0910,0xbfe2,
  762. 0x9903,0x0d8c,0x234b,0xbfba,
  763. 0x50df,0x7f6c,0x6c54,0xbf83,
  764. 0x2ad0,0x2424,0x2afa,0xbf3e,
  765. 0xf631,0xbc17,0xf87a,0xbee7,
  766. 0x6c10,0x501e,0xe81f,0xbe81,
  767. 0x870d,0x5e46,0x5f45,0xbe03,
  768. };
  769. static unsigned short APGD[40] = {
  770. /*0x0000,0x0000,0x0000,0x3ff0,*/
  771. 0x9812,0x060a,0xb7a2,0x4023,
  772. 0xfc96,0x4724,0xa3e3,0x4035,
  773. 0x819a,0xdb2c,0x5025,0x4031,
  774. 0x94e2,0xd5cd,0xb702,0x4018,
  775. 0x1eb1,0x4927,0x6a71,0x3ff1,
  776. 0x7bc5,0x4ad7,0x78de,0x3fb9,
  777. 0xc1d7,0x4b2b,0x991a,0x3f73,
  778. 0xbe1c,0x0b16,0xf98f,0x3f1e,
  779. 0x4ef3,0xfdde,0x10bf,0x3eb7,
  780. 0x647e,0x9dc8,0xe834,0x3e38,
  781. };
  782. #endif
  783. #ifdef MIEEE
  784. static unsigned short APGN[44] = {
  785. 0xbfa2,0x351f,0x5f87,0xaf5b,
  786. 0xbfe4,0x64db,0x1ff7,0x5c76,
  787. 0xbffb,0x564a,0xc221,0x7e49,
  788. 0xbff8,0x0916,0x7f6e,0x0b07,
  789. 0xbfe2,0x0910,0xd8b0,0x6edb,
  790. 0xbfba,0x234b,0x0d8c,0x9903,
  791. 0xbf83,0x6c54,0x7f6c,0x50df,
  792. 0xbf3e,0x2afa,0x2424,0x2ad0,
  793. 0xbee7,0xf87a,0xbc17,0xf631,
  794. 0xbe81,0xe81f,0x501e,0x6c10,
  795. 0xbe03,0x5f45,0x5e46,0x870d,
  796. };
  797. static unsigned short APGD[40] = {
  798. /*0x3ff0,0x0000,0x0000,0x0000,*/
  799. 0x4023,0xb7a2,0x060a,0x9812,
  800. 0x4035,0xa3e3,0x4724,0xfc96,
  801. 0x4031,0x5025,0xdb2c,0x819a,
  802. 0x4018,0xb702,0xd5cd,0x94e2,
  803. 0x3ff1,0x6a71,0x4927,0x1eb1,
  804. 0x3fb9,0x78de,0x4ad7,0x7bc5,
  805. 0x3f73,0x991a,0x4b2b,0xc1d7,
  806. 0x3f1e,0xf98f,0x0b16,0xbe1c,
  807. 0x3eb7,0x10bf,0xfdde,0x4ef3,
  808. 0x3e38,0xe834,0x9dc8,0x647e,
  809. };
  810. #endif
  811. #ifdef ANSIPROT
  812. extern double fabs ( double );
  813. extern double exp ( double );
  814. extern double sqrt ( double );
  815. extern double polevl ( double, void *, int );
  816. extern double p1evl ( double, void *, int );
  817. extern double sin ( double );
  818. extern double cos ( double );
  819. #else
  820. double fabs(), exp(), sqrt();
  821. double polevl(), p1evl(), sin(), cos();
  822. #endif
  823. int airy( x, ai, aip, bi, bip )
  824. double x, *ai, *aip, *bi, *bip;
  825. {
  826. double z, zz, t, f, g, uf, ug, k, zeta, theta;
  827. int domflg;
  828. domflg = 0;
  829. if( x > MAXAIRY )
  830. {
  831. *ai = 0;
  832. *aip = 0;
  833. *bi = MAXNUM;
  834. *bip = MAXNUM;
  835. return(-1);
  836. }
  837. if( x < -2.09 )
  838. {
  839. domflg = 15;
  840. t = sqrt(-x);
  841. zeta = -2.0 * x * t / 3.0;
  842. t = sqrt(t);
  843. k = sqpii / t;
  844. z = 1.0/zeta;
  845. zz = z * z;
  846. uf = 1.0 + zz * polevl( zz, AFN, 8 ) / p1evl( zz, AFD, 9 );
  847. ug = z * polevl( zz, AGN, 10 ) / p1evl( zz, AGD, 10 );
  848. theta = zeta + 0.25 * PI;
  849. f = sin( theta );
  850. g = cos( theta );
  851. *ai = k * (f * uf - g * ug);
  852. *bi = k * (g * uf + f * ug);
  853. uf = 1.0 + zz * polevl( zz, APFN, 8 ) / p1evl( zz, APFD, 9 );
  854. ug = z * polevl( zz, APGN, 10 ) / p1evl( zz, APGD, 10 );
  855. k = sqpii * t;
  856. *aip = -k * (g * uf + f * ug);
  857. *bip = k * (f * uf - g * ug);
  858. return(0);
  859. }
  860. if( x >= 2.09 ) /* cbrt(9) */
  861. {
  862. domflg = 5;
  863. t = sqrt(x);
  864. zeta = 2.0 * x * t / 3.0;
  865. g = exp( zeta );
  866. t = sqrt(t);
  867. k = 2.0 * t * g;
  868. z = 1.0/zeta;
  869. f = polevl( z, AN, 7 ) / polevl( z, AD, 7 );
  870. *ai = sqpii * f / k;
  871. k = -0.5 * sqpii * t / g;
  872. f = polevl( z, APN, 7 ) / polevl( z, APD, 7 );
  873. *aip = f * k;
  874. if( x > 8.3203353 ) /* zeta > 16 */
  875. {
  876. f = z * polevl( z, BN16, 4 ) / p1evl( z, BD16, 5 );
  877. k = sqpii * g;
  878. *bi = k * (1.0 + f) / t;
  879. f = z * polevl( z, BPPN, 4 ) / p1evl( z, BPPD, 5 );
  880. *bip = k * t * (1.0 + f);
  881. return(0);
  882. }
  883. }
  884. f = 1.0;
  885. g = x;
  886. t = 1.0;
  887. uf = 1.0;
  888. ug = x;
  889. k = 1.0;
  890. z = x * x * x;
  891. while( t > MACHEP )
  892. {
  893. uf *= z;
  894. k += 1.0;
  895. uf /=k;
  896. ug *= z;
  897. k += 1.0;
  898. ug /=k;
  899. uf /=k;
  900. f += uf;
  901. k += 1.0;
  902. ug /=k;
  903. g += ug;
  904. t = fabs(uf/f);
  905. }
  906. uf = c1 * f;
  907. ug = c2 * g;
  908. if( (domflg & 1) == 0 )
  909. *ai = uf - ug;
  910. if( (domflg & 2) == 0 )
  911. *bi = sqrt3 * (uf + ug);
  912. /* the deriviative of ai */
  913. k = 4.0;
  914. uf = x * x/2.0;
  915. ug = z/3.0;
  916. f = uf;
  917. g = 1.0 + ug;
  918. uf /= 3.0;
  919. t = 1.0;
  920. while( t > MACHEP )
  921. {
  922. uf *= z;
  923. ug /=k;
  924. k += 1.0;
  925. ug *= z;
  926. uf /=k;
  927. f += uf;
  928. k += 1.0;
  929. ug /=k;
  930. uf /=k;
  931. g += ug;
  932. k += 1.0;
  933. t = fabs(ug/g);
  934. }
  935. uf = c1 * f;
  936. ug = c2 * g;
  937. if( (domflg & 4) == 0 )
  938. *aip = uf - ug;
  939. if( (domflg & 8) == 0 )
  940. *bip = sqrt3 * (uf + ug);
  941. return(0);
  942. }