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