libm/double/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 \f
  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 <math.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 #ifdef ANSIPROT
 579 extern double log ( double );
 580 extern double sin ( double );
 581 extern double cos ( double );
 582 extern double polevl ( double, void *, int );
 583 extern double p1evl ( double, void *, int );
 584 #else
 585 double log(), sin(), cos(), polevl(), p1evl();
 586 #endif
 587 #define EUL 0.57721566490153286061
 588 extern double MAXNUM, PIO2, MACHEP;
 589
 590
 591 int sici( x, si, ci )
 592 double x;
 593 double *si, *ci;
 594 {
 595 double z, c, s, f, g;
 596 short sign;
 597
 598 if( x < 0.0 )
 599         {
 600         sign = -1;
 601         x = -x;
 602         }
 603 else
 604         sign = 0;
 605
 606
 607 if( x == 0.0 )
 608         {
 609         *si = 0.0;
 610         *ci = -MAXNUM;
 611         return( 0 );
 612         }
 613
 614
 615 if( x > 1.0e9 )
 616         {
 617         *si = PIO2 - cos(x)/x;
 618         *ci = sin(x)/x;
 619         return( 0 );
 620         }
 621
 622
 623
 624 if( x > 4.0 )
 625         goto asympt;
 626
 627 z = x * x;
 628 s = x * polevl( z, SN, 5 ) / polevl( z, SD, 5 );
 629 c = z * polevl( z, CN, 5 ) / polevl( z, CD, 5 );
 630
 631 if( sign )
 632         s = -s;
 633 *si = s;
 634 *ci = EUL + log(x) + c; /* real part if x < 0 */
 635 return(0);
 636
 637
 638
 639 /* The auxiliary functions are:
 640  *
 641  *
 642  * *si = *si - PIO2;
 643  * c = cos(x);
 644  * s = sin(x);
 645  *
 646  * t = *ci * s - *si * c;
 647  * a = *ci * c + *si * s;
 648  *
 649  * *si = t;
 650  * *ci = -a;
 651  */
 652
 653
 654 asympt:
 655
 656 s = sin(x);
 657 c = cos(x);
 658 z = 1.0/(x*x);
 659 if( x < 8.0 )
 660         {
 661         f = polevl( z, FN4, 6 ) / (x * p1evl( z, FD4, 7 ));
 662         g = z * polevl( z, GN4, 7 ) / p1evl( z, GD4, 7 );
 663         }
 664 else
 665         {
 666         f = polevl( z, FN8, 8 ) / (x * p1evl( z, FD8, 8 ));
 667         g = z * polevl( z, GN8, 8 ) / p1evl( z, GD8, 9 );
 668         }
 669 *si = PIO2 - f * c - g * s;
 670 if( sign )
 671         *si = -( *si );
 672 *ci = f * s - g * c;
 673
 674 return(0);
 675 }