← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/special/cephes/sici.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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