airy.c 26 KB

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