sici.c 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488
  1. /* sici.c
  2. *
  3. * Sine and cosine integrals
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double x, Ci, Si, sici();
  10. *
  11. * sici( x, &Si, &Ci );
  12. *
  13. *
  14. * DESCRIPTION:
  15. *
  16. * Evaluates the integrals
  17. *
  18. * x
  19. * -
  20. * | cos t - 1
  21. * Ci(x) = eul + ln x + | --------- dt,
  22. * | t
  23. * -
  24. * 0
  25. * x
  26. * -
  27. * | sin t
  28. * Si(x) = | ----- dt
  29. * | t
  30. * -
  31. * 0
  32. *
  33. * where eul = 0.57721566490153286061 is Euler's constant.
  34. * The integrals are approximated by rational functions.
  35. * For x > 8 auxiliary functions f(x) and g(x) are employed
  36. * such that
  37. *
  38. * Ci(x) = f(x) sin(x) - g(x) cos(x)
  39. * Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
  40. *
  41. *
  42. * ACCURACY:
  43. * Test interval = [0,50].
  44. * Absolute error, except relative when > 1:
  45. * arithmetic function # trials peak rms
  46. * IEEE Si 30000 4.4e-16 7.3e-17
  47. * IEEE Ci 30000 6.9e-16 5.1e-17
  48. * DEC Si 5000 4.4e-17 9.0e-18
  49. * DEC Ci 5300 7.9e-17 5.2e-18
  50. */
  51. /*
  52. Cephes Math Library Release 2.1: January, 1989
  53. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  54. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  55. */
  56. #include "mconf.h"
  57. #ifdef UNK
  58. static double SN[] = {
  59. -8.39167827910303881427E-11, 4.62591714427012837309E-8,
  60. -9.75759303843632795789E-6, 9.76945438170435310816E-4,
  61. -4.13470316229406538752E-2, 1.00000000000000000302E0,
  62. };
  63. static double SD[] = {
  64. 2.03269266195951942049E-12, 1.27997891179943299903E-9,
  65. 4.41827842801218905784E-7, 9.96412122043875552487E-5,
  66. 1.42085239326149893930E-2, 9.99999999999999996984E-1,
  67. };
  68. #endif
  69. #ifdef DEC
  70. static unsigned short SN[] = {
  71. 0127670, 0104362, 0167505, 0035161, 0032106, 0127177, 0032131, 0056461,
  72. 0134043, 0132213, 0000476, 0172351, 0035600, 0006331, 0064761, 0032665,
  73. 0137051, 0055601, 0044667, 0017645, 0040200, 0000000, 0000000, 0000000,
  74. };
  75. static unsigned short SD[] = {
  76. 0026417, 0004674, 0052064, 0001573, 0030657, 0165501, 0014666, 0131526,
  77. 0032755, 0032133, 0034147, 0024124, 0034720, 0173167, 0166624, 0154477,
  78. 0036550, 0145336, 0063534, 0063220, 0040200, 0000000, 0000000, 0000000,
  79. };
  80. #endif
  81. #ifdef IBMPC
  82. static unsigned short SN[] = {
  83. 0xa74e, 0x5de8, 0x111e, 0xbdd7, 0x2ba6, 0xe68b, 0xd5cf, 0x3e68,
  84. 0xde9d, 0x6027, 0x7691, 0xbee4, 0x26b7, 0x2d3e, 0x019b, 0x3f50,
  85. 0xe3f5, 0x2936, 0x2b70, 0xbfa5, 0x0000, 0x0000, 0x0000, 0x3ff0,
  86. };
  87. static unsigned short SD[] = {
  88. 0x806f, 0x8a86, 0xe137, 0x3d81, 0xd66b, 0x2336, 0xfd68, 0x3e15,
  89. 0xe50a, 0x670c, 0xa68b, 0x3e9d, 0x9b28, 0xfdb2, 0x1ece, 0x3f1a,
  90. 0x8cd2, 0xcceb, 0x195b, 0x3f8d, 0x0000, 0x0000, 0x0000, 0x3ff0,
  91. };
  92. #endif
  93. #ifdef MIEEE
  94. static unsigned short SN[] = {
  95. 0xbdd7, 0x111e, 0x5de8, 0xa74e, 0x3e68, 0xd5cf, 0xe68b, 0x2ba6,
  96. 0xbee4, 0x7691, 0x6027, 0xde9d, 0x3f50, 0x019b, 0x2d3e, 0x26b7,
  97. 0xbfa5, 0x2b70, 0x2936, 0xe3f5, 0x3ff0, 0x0000, 0x0000, 0x0000,
  98. };
  99. static unsigned short SD[] = {
  100. 0x3d81, 0xe137, 0x8a86, 0x806f, 0x3e15, 0xfd68, 0x2336, 0xd66b,
  101. 0x3e9d, 0xa68b, 0x670c, 0xe50a, 0x3f1a, 0x1ece, 0xfdb2, 0x9b28,
  102. 0x3f8d, 0x195b, 0xcceb, 0x8cd2, 0x3ff0, 0x0000, 0x0000, 0x0000,
  103. };
  104. #endif
  105. #ifdef UNK
  106. static double CN[] = {
  107. 2.02524002389102268789E-11, -1.35249504915790756375E-8,
  108. 3.59325051419993077021E-6, -4.74007206873407909465E-4,
  109. 2.89159652607555242092E-2, -1.00000000000000000080E0,
  110. };
  111. static double CD[] = {
  112. 4.07746040061880559506E-12, 3.06780997581887812692E-9,
  113. 1.23210355685883423679E-6, 3.17442024775032769882E-4,
  114. 5.10028056236446052392E-2, 4.00000000000000000080E0,
  115. };
  116. #endif
  117. #ifdef DEC
  118. static unsigned short CN[] = {
  119. 0027262, 0022131, 0160257, 0020166, 0131550, 0055534, 0077637, 0000557,
  120. 0033561, 0021622, 0161463, 0026575, 0135370, 0102053, 0116333, 0000466,
  121. 0036754, 0160454, 0122022, 0024622, 0140200, 0000000, 0000000, 0000000,
  122. };
  123. static unsigned short CD[] = {
  124. 0026617, 0073177, 0107543, 0104425, 0031122, 0150573, 0156453, 0041517,
  125. 0033245, 0057301, 0077706, 0110510, 0035246, 0067130, 0165424, 0044543,
  126. 0037120, 0164121, 0061206, 0053657, 0040600, 0000000, 0000000, 0000000,
  127. };
  128. #endif
  129. #ifdef IBMPC
  130. static unsigned short CN[] = {
  131. 0xe40f, 0x3c15, 0x448b, 0x3db6, 0xe02e, 0x8ff3, 0x0b6b, 0xbe4d,
  132. 0x65b0, 0x5c66, 0x2472, 0x3ece, 0x6027, 0x739b, 0x1085, 0xbf3f,
  133. 0x4532, 0x9482, 0x9c25, 0x3f9d, 0x0000, 0x0000, 0x0000, 0xbff0,
  134. };
  135. static unsigned short CD[] = {
  136. 0x7123, 0xf1ec, 0xeecf, 0x3d91, 0x686a, 0x7ba5, 0x5a2f, 0x3e2a,
  137. 0xd229, 0x2ff8, 0xabd8, 0x3eb4, 0x892c, 0x1d62, 0xcdcb, 0x3f34,
  138. 0xcaf6, 0x2c50, 0x1d0a, 0x3faa, 0x0000, 0x0000, 0x0000, 0x4010,
  139. };
  140. #endif
  141. #ifdef MIEEE
  142. static unsigned short CN[] = {
  143. 0x3db6, 0x448b, 0x3c15, 0xe40f, 0xbe4d, 0x0b6b, 0x8ff3, 0xe02e,
  144. 0x3ece, 0x2472, 0x5c66, 0x65b0, 0xbf3f, 0x1085, 0x739b, 0x6027,
  145. 0x3f9d, 0x9c25, 0x9482, 0x4532, 0xbff0, 0x0000, 0x0000, 0x0000,
  146. };
  147. static unsigned short CD[] = {
  148. 0x3d91, 0xeecf, 0xf1ec, 0x7123, 0x3e2a, 0x5a2f, 0x7ba5, 0x686a,
  149. 0x3eb4, 0xabd8, 0x2ff8, 0xd229, 0x3f34, 0xcdcb, 0x1d62, 0x892c,
  150. 0x3faa, 0x1d0a, 0x2c50, 0xcaf6, 0x4010, 0x0000, 0x0000, 0x0000,
  151. };
  152. #endif
  153. #ifdef UNK
  154. static double FN4[] = {
  155. 4.23612862892216586994E0, 5.45937717161812843388E0,
  156. 1.62083287701538329132E0, 1.67006611831323023771E-1,
  157. 6.81020132472518137426E-3, 1.08936580650328664411E-4,
  158. 5.48900223421373614008E-7,
  159. };
  160. static double FD4[] = {
  161. /* 1.00000000000000000000E0,*/
  162. 8.16496634205391016773E0, 7.30828822505564552187E0,
  163. 1.86792257950184183883E0, 1.78792052963149907262E-1,
  164. 7.01710668322789753610E-3, 1.10034357153915731354E-4,
  165. 5.48900252756255700982E-7,
  166. };
  167. #endif
  168. #ifdef DEC
  169. static unsigned short FN4[] = {
  170. 0040607, 0107135, 0120133, 0153471, 0040656, 0131467, 0140424,
  171. 0017567, 0040317, 0073563, 0121610, 0002511, 0037453, 0001710,
  172. 0000040, 0006334, 0036337, 0024033, 0176003, 0171425, 0034744,
  173. 0072341, 0121657, 0126035, 0033023, 0054042, 0154652, 0000451,
  174. };
  175. static unsigned short FD4[] = {
  176. /*0040200,0000000,0000000,0000000,*/
  177. 0041002, 0121663, 0137500, 0177450, 0040751, 0156577, 0042213,
  178. 0061552, 0040357, 0014026, 0045465, 0147265, 0037467, 0012503,
  179. 0110413, 0131772, 0036345, 0167701, 0155706, 0160551, 0034746,
  180. 0141076, 0162250, 0123547, 0033023, 0054043, 0056706, 0151050,
  181. };
  182. #endif
  183. #ifdef IBMPC
  184. static unsigned short FN4[] = {
  185. 0x7ae7, 0xb40b, 0xf1cb, 0x4010, 0x83ef, 0xf822, 0xd666,
  186. 0x4015, 0x00a9, 0x7471, 0xeeee, 0x3ff9, 0x019c, 0x0004,
  187. 0x6079, 0x3fc5, 0x7e63, 0x7f80, 0xe503, 0x3f7b, 0xf584,
  188. 0x3475, 0x8e9c, 0x3f1c, 0x4025, 0x5b35, 0x6b04, 0x3ea2,
  189. };
  190. static unsigned short FD4[] = {
  191. /*0x0000,0x0000,0x0000,0x3ff0,*/
  192. 0x1fe5, 0x77e8, 0x5476, 0x4020, 0x6c6d, 0xe891, 0x3baf,
  193. 0x401d, 0xb9d7, 0xc966, 0xe302, 0x3ffd, 0x767f, 0x7221,
  194. 0xe2a8, 0x3fc6, 0xdc2d, 0x3b78, 0xbdf8, 0x3f7c, 0x14ed,
  195. 0xdc95, 0xd847, 0x3f1c, 0xda45, 0x6bb8, 0x6b04, 0x3ea2,
  196. };
  197. #endif
  198. #ifdef MIEEE
  199. static unsigned short FN4[] = {
  200. 0x4010, 0xf1cb, 0xb40b, 0x7ae7, 0x4015, 0xd666, 0xf822,
  201. 0x83ef, 0x3ff9, 0xeeee, 0x7471, 0x00a9, 0x3fc5, 0x6079,
  202. 0x0004, 0x019c, 0x3f7b, 0xe503, 0x7f80, 0x7e63, 0x3f1c,
  203. 0x8e9c, 0x3475, 0xf584, 0x3ea2, 0x6b04, 0x5b35, 0x4025,
  204. };
  205. static unsigned short FD4[] = {
  206. /* 0x3ff0,0x0000,0x0000,0x0000,*/
  207. 0x4020, 0x5476, 0x77e8, 0x1fe5, 0x401d, 0x3baf, 0xe891,
  208. 0x6c6d, 0x3ffd, 0xe302, 0xc966, 0xb9d7, 0x3fc6, 0xe2a8,
  209. 0x7221, 0x767f, 0x3f7c, 0xbdf8, 0x3b78, 0xdc2d, 0x3f1c,
  210. 0xd847, 0xdc95, 0x14ed, 0x3ea2, 0x6b04, 0x6bb8, 0xda45,
  211. };
  212. #endif
  213. #ifdef UNK
  214. static double FN8[] = {
  215. 4.55880873470465315206E-1, 7.13715274100146711374E-1,
  216. 1.60300158222319456320E-1, 1.16064229408124407915E-2,
  217. 3.49556442447859055605E-4, 4.86215430826454749482E-6,
  218. 3.20092790091004902806E-8, 9.41779576128512936592E-11,
  219. 9.70507110881952024631E-14,
  220. };
  221. static double FD8[] = {
  222. /* 1.00000000000000000000E0,*/
  223. 9.17463611873684053703E-1, 1.78685545332074536321E-1,
  224. 1.22253594771971293032E-2, 3.58696481881851580297E-4,
  225. 4.92435064317881464393E-6, 3.21956939101046018377E-8,
  226. 9.43720590350276732376E-11, 9.70507110881952025725E-14,
  227. };
  228. #endif
  229. #ifdef DEC
  230. static unsigned short FN8[] = {
  231. 0037751, 0064467, 0142332, 0164573, 0040066, 0133013, 0050352, 0071102,
  232. 0037444, 0022671, 0102157, 0013535, 0036476, 0024335, 0136423, 0146444,
  233. 0035267, 0042253, 0164110, 0110460, 0033643, 0022626, 0062535, 0060320,
  234. 0032011, 0075223, 0010110, 0153413, 0027717, 0014572, 0011360, 0014034,
  235. 0025332, 0104755, 0004563, 0152354,
  236. };
  237. static unsigned short FD8[] = {
  238. /*0040200,0000000,0000000,0000000,*/
  239. 0040152, 0157345, 0030104, 0075616, 0037466, 0174527, 0172740, 0071060,
  240. 0036510, 0046337, 0144272, 0156552, 0035274, 0007555, 0042537, 0015572,
  241. 0033645, 0035731, 0112465, 0026474, 0032012, 0043612, 0030613, 0030123,
  242. 0027717, 0103277, 0004564, 0151000, 0025332, 0104755, 0004563, 0152354,
  243. };
  244. #endif
  245. #ifdef IBMPC
  246. static unsigned short FN8[] = {
  247. 0x5d2f, 0xf89b, 0x2d26, 0x3fdd, 0x4e48, 0x6a1d, 0xd6c1, 0x3fe6, 0xe2ec,
  248. 0x308d, 0x84b7, 0x3fc4, 0x79a4, 0xb7a2, 0xc51b, 0x3f87, 0x1226, 0x7d09,
  249. 0xe895, 0x3f36, 0xac1a, 0xccab, 0x64b2, 0x3ed4, 0x1ae1, 0x6209, 0x2f52,
  250. 0x3e61, 0x0304, 0x425e, 0xe32f, 0x3dd9, 0x7a9d, 0xa12e, 0x513d, 0x3d3b,
  251. };
  252. static unsigned short FD8[] = {
  253. /*0x0000,0x0000,0x0000,0x3ff0,*/
  254. 0x8f72, 0xa608, 0x5bdc, 0x3fed, 0x0e46, 0xfebc, 0xdf2a, 0x3fc6,
  255. 0x5bad, 0xf917, 0x099b, 0x3f89, 0xe36f, 0xa8ab, 0x81ed, 0x3f37,
  256. 0xa5a8, 0x32a6, 0xa77b, 0x3ed4, 0x660a, 0x4631, 0x48f1, 0x3e61,
  257. 0x9a40, 0xe12e, 0xf0d7, 0x3dd9, 0x7a9d, 0xa12e, 0x513d, 0x3d3b,
  258. };
  259. #endif
  260. #ifdef MIEEE
  261. static unsigned short FN8[] = {
  262. 0x3fdd, 0x2d26, 0xf89b, 0x5d2f, 0x3fe6, 0xd6c1, 0x6a1d, 0x4e48, 0x3fc4,
  263. 0x84b7, 0x308d, 0xe2ec, 0x3f87, 0xc51b, 0xb7a2, 0x79a4, 0x3f36, 0xe895,
  264. 0x7d09, 0x1226, 0x3ed4, 0x64b2, 0xccab, 0xac1a, 0x3e61, 0x2f52, 0x6209,
  265. 0x1ae1, 0x3dd9, 0xe32f, 0x425e, 0x0304, 0x3d3b, 0x513d, 0xa12e, 0x7a9d,
  266. };
  267. static unsigned short FD8[] = {
  268. /*0x3ff0,0x0000,0x0000,0x0000,*/
  269. 0x3fed, 0x5bdc, 0xa608, 0x8f72, 0x3fc6, 0xdf2a, 0xfebc, 0x0e46,
  270. 0x3f89, 0x099b, 0xf917, 0x5bad, 0x3f37, 0x81ed, 0xa8ab, 0xe36f,
  271. 0x3ed4, 0xa77b, 0x32a6, 0xa5a8, 0x3e61, 0x48f1, 0x4631, 0x660a,
  272. 0x3dd9, 0xf0d7, 0xe12e, 0x9a40, 0x3d3b, 0x513d, 0xa12e, 0x7a9d,
  273. };
  274. #endif
  275. #ifdef UNK
  276. static double GN4[] = {
  277. 8.71001698973114191777E-2, 6.11379109952219284151E-1,
  278. 3.97180296392337498885E-1, 7.48527737628469092119E-2,
  279. 5.38868681462177273157E-3, 1.61999794598934024525E-4,
  280. 1.97963874140963632189E-6, 7.82579040744090311069E-9,
  281. };
  282. static double GD4[] = {
  283. /* 1.00000000000000000000E0,*/
  284. 1.64402202413355338886E0, 6.66296701268987968381E-1,
  285. 9.88771761277688796203E-2, 6.22396345441768420760E-3,
  286. 1.73221081474177119497E-4, 2.02659182086343991969E-6,
  287. 7.82579218933534490868E-9,
  288. };
  289. #endif
  290. #ifdef DEC
  291. static unsigned short GN4[] = {
  292. 0037262, 0060622, 0164572, 0157515, 0040034, 0101527, 0061263, 0147204,
  293. 0037713, 0055467, 0037475, 0144512, 0037231, 0046151, 0035234, 0045261,
  294. 0036260, 0111624, 0150617, 0053536, 0035051, 0157175, 0016675, 0155456,
  295. 0033404, 0154757, 0041211, 0000055, 0031406, 0071060, 0130322, 0033322,
  296. };
  297. static unsigned short GD4[] = {
  298. /* 0040200,0000000,0000000,0000000,*/
  299. 0040322, 0067520, 0046707, 0053275, 0040052, 0111153, 0126542,
  300. 0005516, 0037312, 0100035, 0167121, 0014552, 0036313, 0171143,
  301. 0137176, 0014213, 0035065, 0121256, 0012033, 0150603, 0033410,
  302. 0000225, 0013121, 0071643, 0031406, 0071062, 0131152, 0150454,
  303. };
  304. #endif
  305. #ifdef IBMPC
  306. static unsigned short GN4[] = {
  307. 0x5bea, 0x5d2f, 0x4c32, 0x3fb6, 0x79d1, 0xec56, 0x906a, 0x3fe3,
  308. 0xb929, 0xe7e7, 0x6b66, 0x3fd9, 0x8956, 0x2753, 0x298d, 0x3fb3,
  309. 0xeaec, 0x9a31, 0x1272, 0x3f76, 0xbb66, 0xa3b7, 0x3bcf, 0x3f25,
  310. 0x2006, 0xe851, 0x9b3d, 0x3ec0, 0x46da, 0x161a, 0xce46, 0x3e40,
  311. };
  312. static unsigned short GD4[] = {
  313. /* 0x0000,0x0000,0x0000,0x3ff0,*/
  314. 0xead8, 0x09b8, 0x4dea, 0x3ffa, 0x416a, 0x75ac, 0x524d,
  315. 0x3fe5, 0x232d, 0xbdca, 0x5003, 0x3fb9, 0xc311, 0x77cf,
  316. 0x7e4c, 0x3f79, 0x7a30, 0xc283, 0xb455, 0x3f26, 0x2e74,
  317. 0xa2ca, 0x0012, 0x3ec1, 0x5a26, 0x564d, 0xce46, 0x3e40,
  318. };
  319. #endif
  320. #ifdef MIEEE
  321. static unsigned short GN4[] = {
  322. 0x3fb6, 0x4c32, 0x5d2f, 0x5bea, 0x3fe3, 0x906a, 0xec56, 0x79d1,
  323. 0x3fd9, 0x6b66, 0xe7e7, 0xb929, 0x3fb3, 0x298d, 0x2753, 0x8956,
  324. 0x3f76, 0x1272, 0x9a31, 0xeaec, 0x3f25, 0x3bcf, 0xa3b7, 0xbb66,
  325. 0x3ec0, 0x9b3d, 0xe851, 0x2006, 0x3e40, 0xce46, 0x161a, 0x46da,
  326. };
  327. static unsigned short GD4[] = {
  328. /*0x3ff0,0x0000,0x0000,0x0000,*/
  329. 0x3ffa, 0x4dea, 0x09b8, 0xead8, 0x3fe5, 0x524d, 0x75ac,
  330. 0x416a, 0x3fb9, 0x5003, 0xbdca, 0x232d, 0x3f79, 0x7e4c,
  331. 0x77cf, 0xc311, 0x3f26, 0xb455, 0xc283, 0x7a30, 0x3ec1,
  332. 0x0012, 0xa2ca, 0x2e74, 0x3e40, 0xce46, 0x564d, 0x5a26,
  333. };
  334. #endif
  335. #ifdef UNK
  336. static double GN8[] = {
  337. 6.97359953443276214934E-1, 3.30410979305632063225E-1,
  338. 3.84878767649974295920E-2, 1.71718239052347903558E-3,
  339. 3.48941165502279436777E-5, 3.47131167084116673800E-7,
  340. 1.70404452782044526189E-9, 3.85945925430276600453E-12,
  341. 3.14040098946363334640E-15,
  342. };
  343. static double GD8[] = {
  344. /* 1.00000000000000000000E0,*/
  345. 1.68548898811011640017E0, 4.87852258695304967486E-1,
  346. 4.67913194259625806320E-2, 1.90284426674399523638E-3,
  347. 3.68475504442561108162E-5, 3.57043223443740838771E-7,
  348. 1.72693748966316146736E-9, 3.87830166023954706752E-12,
  349. 3.14040098946363335242E-15,
  350. };
  351. #endif
  352. #ifdef DEC
  353. static unsigned short GN8[] = {
  354. 0040062, 0103056, 0110624, 0033123, 0037651, 0025640, 0136266, 0145647,
  355. 0037035, 0122566, 0137770, 0061777, 0035741, 0011424, 0065311, 0013370,
  356. 0034422, 0055505, 0134324, 0016755, 0032672, 0056530, 0022565, 0014747,
  357. 0030752, 0031674, 0114735, 0013162, 0026607, 0145353, 0022020, 0123625,
  358. 0024142, 0045054, 0060033, 0016505,
  359. };
  360. static unsigned short GD8[] = {
  361. /*0040200,0000000,0000000,0000000,*/
  362. 0040327, 0137032, 0064331, 0136425, 0037771, 0143705, 0070300, 0105711,
  363. 0037077, 0124101, 0025275, 0035356, 0035771, 0064333, 0145103, 0105357,
  364. 0034432, 0106301, 0105311, 0010713, 0032677, 0127645, 0120034, 0157551,
  365. 0030755, 0054466, 0010743, 0105566, 0026610, 0072242, 0142530, 0135744,
  366. 0024142, 0045054, 0060033, 0016505,
  367. };
  368. #endif
  369. #ifdef IBMPC
  370. static unsigned short GN8[] = {
  371. 0x86ca, 0xd232, 0x50c5, 0x3fe6, 0xd975, 0x1796, 0x2574, 0x3fd5, 0x0c80,
  372. 0xd7ff, 0xb4ae, 0x3fa3, 0x22df, 0x8d59, 0x2262, 0x3f5c, 0x83be, 0xb71a,
  373. 0x4b68, 0x3f02, 0xa33d, 0x04ae, 0x4bab, 0x3e97, 0xa2ce, 0x933b, 0x4677,
  374. 0x3e1d, 0x14f3, 0x6482, 0xf95d, 0x3d90, 0x63a9, 0x8c03, 0x4945, 0x3cec,
  375. };
  376. static unsigned short GD8[] = {
  377. /*0x0000,0x0000,0x0000,0x3ff0,*/
  378. 0x37a3, 0x4d1b, 0xf7c3, 0x3ffa, 0x1179, 0xae18, 0x38f8, 0x3fdf, 0xa75e,
  379. 0x2557, 0xf508, 0x3fa7, 0x715e, 0x7948, 0x2d1b, 0x3f5f, 0x2239, 0x3159,
  380. 0x5198, 0x3f03, 0x9bed, 0xb403, 0xf5f4, 0x3e97, 0x716f, 0xc23c, 0xab26,
  381. 0x3e1d, 0x177c, 0x58ab, 0x0e94, 0x3d91, 0x63a9, 0x8c03, 0x4945, 0x3cec,
  382. };
  383. #endif
  384. #ifdef MIEEE
  385. static unsigned short GN8[] = {
  386. 0x3fe6, 0x50c5, 0xd232, 0x86ca, 0x3fd5, 0x2574, 0x1796, 0xd975, 0x3fa3,
  387. 0xb4ae, 0xd7ff, 0x0c80, 0x3f5c, 0x2262, 0x8d59, 0x22df, 0x3f02, 0x4b68,
  388. 0xb71a, 0x83be, 0x3e97, 0x4bab, 0x04ae, 0xa33d, 0x3e1d, 0x4677, 0x933b,
  389. 0xa2ce, 0x3d90, 0xf95d, 0x6482, 0x14f3, 0x3cec, 0x4945, 0x8c03, 0x63a9,
  390. };
  391. static unsigned short GD8[] = {
  392. /*0x3ff0,0x0000,0x0000,0x0000,*/
  393. 0x3ffa, 0xf7c3, 0x4d1b, 0x37a3, 0x3fdf, 0x38f8, 0xae18, 0x1179, 0x3fa7,
  394. 0xf508, 0x2557, 0xa75e, 0x3f5f, 0x2d1b, 0x7948, 0x715e, 0x3f03, 0x5198,
  395. 0x3159, 0x2239, 0x3e97, 0xf5f4, 0xb403, 0x9bed, 0x3e1d, 0xab26, 0xc23c,
  396. 0x716f, 0x3d91, 0x0e94, 0x58ab, 0x177c, 0x3cec, 0x4945, 0x8c03, 0x63a9,
  397. };
  398. #endif
  399. #ifdef ANSIPROT
  400. extern double log(double);
  401. extern double sin(double);
  402. extern double cos(double);
  403. extern double polevl(double, void *, int);
  404. extern double p1evl(double, void *, int);
  405. #else
  406. double log(), sin(), cos(), polevl(), p1evl();
  407. #endif
  408. #define EUL 0.57721566490153286061
  409. extern double MAXNUM, PIO2, MACHEP;
  410. int sici(x, si, ci) double x;
  411. double *si, *ci;
  412. {
  413. double z, c, s, f, g;
  414. short sign;
  415. if (x < 0.0) {
  416. sign = -1;
  417. x = -x;
  418. } else
  419. sign = 0;
  420. if (x == 0.0) {
  421. *si = 0.0;
  422. *ci = -MAXNUM;
  423. return (0);
  424. }
  425. if (x > 1.0e9) {
  426. *si = PIO2 - cos(x) / x;
  427. *ci = sin(x) / x;
  428. return (0);
  429. }
  430. if (x > 4.0)
  431. goto asympt;
  432. z = x * x;
  433. s = x * polevl(z, SN, 5) / polevl(z, SD, 5);
  434. c = z * polevl(z, CN, 5) / polevl(z, CD, 5);
  435. if (sign)
  436. s = -s;
  437. *si = s;
  438. *ci = EUL + log(x) + c; /* real part if x < 0 */
  439. return (0);
  440. /* The auxiliary functions are:
  441. *
  442. *
  443. * *si = *si - PIO2;
  444. * c = cos(x);
  445. * s = sin(x);
  446. *
  447. * t = *ci * s - *si * c;
  448. * a = *ci * c + *si * s;
  449. *
  450. * *si = t;
  451. * *ci = -a;
  452. */
  453. asympt:
  454. s = sin(x);
  455. c = cos(x);
  456. z = 1.0 / (x * x);
  457. if (x < 8.0) {
  458. f = polevl(z, FN4, 6) / (x * p1evl(z, FD4, 7));
  459. g = z * polevl(z, GN4, 7) / p1evl(z, GD4, 7);
  460. } else {
  461. f = polevl(z, FN8, 8) / (x * p1evl(z, FD8, 8));
  462. g = z * polevl(z, GN8, 8) / p1evl(z, GD8, 9);
  463. }
  464. *si = PIO2 - f * c - g * s;
  465. if (sign)
  466. *si = -(*si);
  467. *ci = f * s - g * c;
  468. return (0);
  469. }