lib/libm/ei.c

   1 /*                                                      ei.c
   2  *
   3  *      Exponential integral
   4  *
   5  *
   6  * SYNOPSIS:
   7  *
   8  * double x, y, ei();
   9  *
  10  * y = ei( x );
  11  *
  12  *
  13  *
  14  * DESCRIPTION:
  15  *
  16  *               x
  17  *                -     t
  18  *               | |   e
  19  *    Ei(x) =   -|-   ---  dt .
  20  *             | |     t
  21  *              -
  22  *             -inf
  23  *
  24  * Not defined for x <= 0.
  25  * See also expn.c.
  26  *
  27  *
  28  *
  29  * ACCURACY:
  30  *
  31  *                      Relative error:
  32  * arithmetic   domain     # trials      peak         rms
  33  *    IEEE       0,100       50000      8.6e-16     1.3e-16
  34  *
  35  */
  36
  37 /*
  38 Cephes Math Library Release 2.8:  May, 1999
  39 Copyright 1999 by Stephen L. Moshier
  40 */
  41
  42 #include "mconf.h"
  43 #ifdef ANSIPROT
  44 extern double log ( double );
  45 extern double exp ( double );
  46 extern double polevl ( double, void *, int );
  47 extern double p1evl ( double, void *, int );
  48 #else
  49 extern double log(), exp(), polevl(), p1evl();
  50 #endif
  51
  52 #define EUL 5.772156649015328606065e-1
  53
  54 /* 0 < x <= 2
  55    Ei(x) - EUL - ln(x) = x A(x)/B(x)
  56    Theoretical peak relative error 9.73e-18  */
  57 #if UNK
  58 static double A[6] = {
  59 -5.350447357812542947283E0,
  60  2.185049168816613393830E2,
  61 -4.176572384826693777058E3,
  62  5.541176756393557601232E4,
  63 -3.313381331178144034309E5,
  64  1.592627163384945414220E6,
  65 };
  66 static double B[6] = {
  67   /*  1.000000000000000000000E0, */
  68 -5.250547959112862969197E1,
  69  1.259616186786790571525E3,
  70 -1.756549581973534652631E4,
  71  1.493062117002725991967E5,
  72 -7.294949239640527645655E5,
  73  1.592627163384945429726E6,
  74 };
  75 #endif
  76 #if DEC
  77 static short A[24] = {
  78 0140653,0033335,0060230,0144217,
  79 0042132,0100502,0035625,0167413,
  80 0143202,0102224,0037176,0175403,
  81 0044130,0071704,0077421,0170343,
  82 0144641,0144504,0041200,0045154,
  83 0045302,0064631,0047234,0142052,
  84 };
  85 static short B[24] = {
  86   /* 0040200,0000000,0000000,0000000, */
  87 0141522,0002634,0070442,0142614,
  88 0042635,0071667,0146532,0027705,
  89 0143611,0035375,0156025,0114015,
  90 0044421,0147215,0106177,0046330,
  91 0145062,0014556,0144216,0103725,
  92 0045302,0064631,0047234,0142052,
  93 };
  94 #endif
  95 #if IBMPC
  96 static short A[24] = {
  97 0x1912,0xac13,0x66db,0xc015,
  98 0xbde1,0x4772,0x5028,0x406b,
  99 0xdf60,0x87cf,0x5092,0xc0b0,
 100 0x3e1c,0x8fe2,0x0e78,0x40eb,
 101 0x094e,0x8850,0x3928,0xc114,
 102 0x9885,0x29d3,0x4d33,0x4138,
 103 };
 104 static short B[24] = {
 105   /* 0x0000,0x0000,0x0000,0x3ff0, */
 106 0x58b1,0x8e24,0x40b3,0xc04a,
 107 0x45f9,0xf9ab,0xae76,0x4093,
 108 0xb302,0xbb82,0x275f,0xc0d1,
 109 0xe99b,0xb18f,0x39d1,0x4102,
 110 0xd0fb,0xd911,0x432d,0xc126,
 111 0x9885,0x29d3,0x4d33,0x4138,
 112 };
 113 #endif
 114 #if MIEEE
 115 static short A[24] = {
 116 0xc015,0x66db,0xac13,0x1912,
 117 0x406b,0x5028,0x4772,0xbde1,
 118 0xc0b0,0x5092,0x87cf,0xdf60,
 119 0x40eb,0x0e78,0x8fe2,0x3e1c,
 120 0xc114,0x3928,0x8850,0x094e,
 121 0x4138,0x4d33,0x29d3,0x9885,
 122 };
 123 static short B[24] = {
 124   /* 0x3ff0,0x0000,0x0000,0x0000, */
 125 0xc04a,0x40b3,0x8e24,0x58b1,
 126 0x4093,0xae76,0xf9ab,0x45f9,
 127 0xc0d1,0x275f,0xbb82,0xb302,
 128 0x4102,0x39d1,0xb18f,0xe99b,
 129 0xc126,0x432d,0xd911,0xd0fb,
 130 0x4138,0x4d33,0x29d3,0x9885,
 131 };
 132 #endif
 133
 134 #if 0
 135 /* 0 < x <= 4
 136    Ei(x) - EUL - ln(x) = x A(x)/B(x)
 137    Theoretical peak relative error 4.75e-17  */
 138 #if UNK
 139 static double A[7] = {
 140 -6.831869820732773831942E0,
 141  2.920190530726774500309E2,
 142 -1.195883839286649567993E4,
 143  1.761045255472548975666E5,
 144 -2.623034438354006526979E6,
 145  1.472430336917880803157E7,
 146 -8.205359388213261174960E7,
 147 };
 148 static double B[7] = {
 149   /* 1.000000000000000000000E0, */
 150 -7.731946237840033971071E1,
 151  2.751808700543578450827E3,
 152 -5.829268609072186897994E4,
 153  7.916610857961870631379E5,
 154 -6.873926904825733094076E6,
 155  3.523770183971164032710E7,
 156 -8.205359388213260785363E7,
 157 };
 158 #endif
 159 #if DEC
 160 static short A[28] = {
 161 0140732,0117255,0072522,0071743,
 162 0042222,0001160,0052302,0002334,
 163 0143472,0155532,0101650,0155462,
 164 0044453,0175041,0121220,0172022,
 165 0145440,0014351,0140337,0157550,
 166 0046140,0126317,0057202,0100233,
 167 0146634,0100473,0036072,0067054,
 168 };
 169 static short B[28] = {
 170   /* 0040200,0000000,0000000,0000000, */
 171 0141632,0121620,0111247,0010115,
 172 0043053,0176360,0067773,0027324,
 173 0144143,0132257,0121644,0036204,
 174 0045101,0043321,0057553,0151231,
 175 0145721,0143215,0147505,0050610,
 176 0046406,0065721,0072675,0152744,
 177 0146634,0100473,0036072,0067052,
 178 };
 179 #endif
 180 #if IBMPC
 181 static short A[28] = {
 182 0x4e7c,0xaeaa,0x53d5,0xc01b,
 183 0x409b,0x0a98,0x404e,0x4072,
 184 0x1b66,0x5075,0x5b6b,0xc0c7,
 185 0x1e82,0x3452,0x7f44,0x4105,
 186 0xfbed,0x381b,0x031d,0xc144,
 187 0x5013,0xebd0,0x1599,0x416c,
 188 0x4dc5,0x6787,0x9027,0xc193,
 189 };
 190 static short B[28] = {
 191   /* 0x0000,0x0000,0x0000,0x3ff0, */
 192 0xe20a,0x1254,0x5472,0xc053,
 193 0x65db,0x0dff,0x7f9e,0x40a5,
 194 0x8791,0xf474,0x7695,0xc0ec,
 195 0x7a53,0x2bed,0x28da,0x4128,
 196 0xaa31,0xb9e8,0x38d1,0xc15a,
 197 0xbabd,0x2eb7,0xcd7a,0x4180,
 198 0x4dc5,0x6787,0x9027,0xc193,
 199 };
 200 #endif
 201 #if MIEEE
 202 static short A[28] = {
 203 0xc01b,0x53d5,0xaeaa,0x4e7c,
 204 0x4072,0x404e,0x0a98,0x409b,
 205 0xc0c7,0x5b6b,0x5075,0x1b66,
 206 0x4105,0x7f44,0x3452,0x1e82,
 207 0xc144,0x031d,0x381b,0xfbed,
 208 0x416c,0x1599,0xebd0,0x5013,
 209 0xc193,0x9027,0x6787,0x4dc5,
 210 };
 211 static short B[28] = {
 212   /* 0x3ff0,0x0000,0x0000,0x0000, */
 213 0xc053,0x5472,0x1254,0xe20a,
 214 0x40a5,0x7f9e,0x0dff,0x65db,
 215 0xc0ec,0x7695,0xf474,0x8791,
 216 0x4128,0x28da,0x2bed,0x7a53,
 217 0xc15a,0x38d1,0xb9e8,0xaa31,
 218 0x4180,0xcd7a,0x2eb7,0xbabd,
 219 0xc193,0x9027,0x6787,0x4dc5,
 220 };
 221 #endif
 222 #endif /* 0 */
 223
 224 #if 0
 225 /* 0 < x <= 8
 226    Ei(x) - EUL - ln(x) = x A(x)/B(x)
 227    Theoretical peak relative error 2.14e-17  */
 228
 229 #if UNK
 230 static double A[9] = {
 231 -1.111230942210860450145E1,
 232  3.688203982071386319616E2,
 233 -4.924786153494029574350E4,
 234  1.050677503345557903241E6,
 235 -3.626713709916703688968E7,
 236  4.353499908839918635414E8,
 237 -6.454613717232006895409E9,
 238  3.408243056457762907071E10,
 239 -1.995466674647028468613E11,
 240 };
 241 static double B[9] = {
 242   /*  1.000000000000000000000E0, */
 243 -1.356757648138514017969E2,
 244  8.562181317107341736606E3,
 245 -3.298257180413775117555E5,
 246  8.543534058481435917210E6,
 247 -1.542380618535140055068E8,
 248  1.939251779195993632028E9,
 249 -1.636096210465615015435E10,
 250  8.396909743075306970605E10,
 251 -1.995466674647028425886E11,
 252 };
 253 #endif
 254 #if DEC
 255 static short A[36] = {
 256 0141061,0146004,0173357,0151553,
 257 0042270,0064402,0147366,0126701,
 258 0144100,0057734,0106615,0144356,
 259 0045200,0040654,0003332,0004456,
 260 0146412,0054440,0043130,0140263,
 261 0047317,0113517,0033422,0065123,
 262 0150300,0056313,0065235,0131147,
 263 0050775,0167423,0146222,0075760,
 264 0151471,0153642,0003442,0147667,
 265 };
 266 static short B[36] = {
 267   /* 0040200,0000000,0000000,0000000, */
 268 0142007,0126376,0166077,0043600,
 269 0043405,0144271,0125461,0014364,
 270 0144641,0006066,0175061,0164463,
 271 0046002,0056456,0007370,0121657,
 272 0147023,0013706,0156647,0177115,
 273 0047747,0026504,0103144,0054507,
 274 0150563,0146036,0007051,0177135,
 275 0051234,0063625,0173266,0003111,
 276 0151471,0153642,0003442,0147666,
 277 };
 278 #endif
 279 #if IBMPC
 280 static short A[36] = {
 281 0xfa6d,0x9edd,0x3980,0xc026,
 282 0xd5b8,0x59de,0x0d20,0x4077,
 283 0xb91e,0x91b1,0x0bfb,0xc0e8,
 284 0x4126,0x80db,0x0835,0x4130,
 285 0x1816,0x08cb,0x4b24,0xc181,
 286 0x4d4a,0xe6e2,0xf2e9,0x41b9,
 287 0xb64d,0x6d53,0x0b99,0xc1f8,
 288 0x4f7e,0x7992,0xbde2,0x421f,
 289 0x59f7,0x40e4,0x3af4,0xc247,
 290 };
 291 static short B[36] = {
 292   /* 0x0000,0x0000,0x0000,0x3ff0, */
 293 0xe8f0,0xdd87,0xf59f,0xc060,
 294 0x231e,0x3566,0xb917,0x40c0,
 295 0x3d26,0xdf46,0x2186,0xc114,
 296 0x1476,0xc1df,0x4ba5,0x4160,
 297 0xffca,0xdbb4,0x62f8,0xc1a2,
 298 0x8b29,0x90cc,0xe5a8,0x41dc,
 299 0x3fcc,0xc1c5,0x7983,0xc20e,
 300 0xc0c9,0xbed6,0x8cf2,0x4233,
 301 0x59f7,0x40e4,0x3af4,0xc247,
 302 };
 303 #endif
 304 #if MIEEE
 305 static short A[36] = {
 306 0xc026,0x3980,0x9edd,0xfa6d,
 307 0x4077,0x0d20,0x59de,0xd5b8,
 308 0xc0e8,0x0bfb,0x91b1,0xb91e,
 309 0x4130,0x0835,0x80db,0x4126,
 310 0xc181,0x4b24,0x08cb,0x1816,
 311 0x41b9,0xf2e9,0xe6e2,0x4d4a,
 312 0xc1f8,0x0b99,0x6d53,0xb64d,
 313 0x421f,0xbde2,0x7992,0x4f7e,
 314 0xc247,0x3af4,0x40e4,0x59f7,
 315 };
 316 static short B[36] = {
 317   /* 0x3ff0,0x0000,0x0000,0x0000, */
 318 0xc060,0xf59f,0xdd87,0xe8f0,
 319 0x40c0,0xb917,0x3566,0x231e,
 320 0xc114,0x2186,0xdf46,0x3d26,
 321 0x4160,0x4ba5,0xc1df,0x1476,
 322 0xc1a2,0x62f8,0xdbb4,0xffca,
 323 0x41dc,0xe5a8,0x90cc,0x8b29,
 324 0xc20e,0x7983,0xc1c5,0x3fcc,
 325 0x4233,0x8cf2,0xbed6,0xc0c9,
 326 0xc247,0x3af4,0x40e4,0x59f7,
 327 };
 328 #endif
 329 #endif /* 0 */
 330
 331 /* 8 <= x <= 20
 332    x exp(-x) Ei(x) - 1 = 1/x R(1/x)
 333    Theoretical peak absolute error = 1.07e-17  */
 334 #if UNK
 335 static double A2[10] = {
 336 -2.106934601691916512584E0,
 337  1.732733869664688041885E0,
 338 -2.423619178935841904839E-1,
 339  2.322724180937565842585E-2,
 340  2.372880440493179832059E-4,
 341 -8.343219561192552752335E-5,
 342  1.363408795605250394881E-5,
 343 -3.655412321999253963714E-7,
 344  1.464941733975961318456E-8,
 345  6.176407863710360207074E-10,
 346 };
 347 static double B2[9] = {
 348   /* 1.000000000000000000000E0, */
 349 -2.298062239901678075778E-1,
 350  1.105077041474037862347E-1,
 351 -1.566542966630792353556E-2,
 352  2.761106850817352773874E-3,
 353 -2.089148012284048449115E-4,
 354  1.708528938807675304186E-5,
 355 -4.459311796356686423199E-7,
 356  1.394634930353847498145E-8,
 357  6.150865933977338354138E-10,
 358 };
 359 #endif
 360 #if DEC
 361 static short A2[40] = {
 362 0140406,0154004,0035104,0173336,
 363 0040335,0145071,0031560,0150165,
 364 0137570,0026670,0176230,0055040,
 365 0036676,0043416,0077122,0054476,
 366 0035170,0150206,0034407,0175571,
 367 0134656,0174121,0123231,0021751,
 368 0034144,0136766,0036746,0121115,
 369 0132704,0037632,0135077,0107300,
 370 0031573,0126321,0117076,0004314,
 371 0030451,0143233,0041352,0172464,
 372 };
 373 static short B2[36] = {
 374   /* 0040200,0000000,0000000,0000000, */
 375 0137553,0051122,0120721,0170437,
 376 0037342,0050734,0175047,0032132,
 377 0136600,0052311,0101406,0147050,
 378 0036064,0171657,0120001,0071165,
 379 0135133,0010043,0151244,0066340,
 380 0034217,0051141,0026115,0043305,
 381 0132757,0064120,0106341,0051217,
 382 0031557,0114261,0060663,0135017,
 383 0030451,0011337,0001344,0175542,
 384 };
 385 #endif
 386 #if IBMPC
 387 static short A2[40] = {
 388 0x9edc,0x8748,0xdb00,0xc000,
 389 0x1a0f,0x266e,0xb947,0x3ffb,
 390 0x0b44,0x1f93,0x05b7,0xbfcf,
 391 0x4b28,0xcfca,0xc8e1,0x3f97,
 392 0xff6f,0xc720,0x1a10,0x3f2f,
 393 0x247d,0x34d3,0xdf0a,0xbf15,
 394 0xd44a,0xc7bc,0x97be,0x3eec,
 395 0xf1d8,0x5747,0x87f3,0xbe98,
 396 0xc119,0x33c7,0x759a,0x3e4f,
 397 0x5ea6,0x685d,0x38d3,0x3e05,
 398 };
 399 static short B2[36] = {
 400   /* 0x0000,0x0000,0x0000,0x3ff0, */
 401 0x3e24,0x543a,0x6a4a,0xbfcd,
 402 0xe68b,0x9f44,0x4a3b,0x3fbc,
 403 0xd9c5,0x3060,0x0a99,0xbf90,
 404 0x2e4f,0xf400,0x9e75,0x3f66,
 405 0x8d9c,0x7a54,0x6204,0xbf2b,
 406 0xa8d9,0x2589,0xea4c,0x3ef1,
 407 0x2a52,0x119c,0xed0a,0xbe9d,
 408 0x7742,0x2c36,0xf316,0x3e4d,
 409 0x9f6c,0xe05c,0x225b,0x3e05,
 410 };
 411 #endif
 412 #if MIEEE
 413 static short A2[40] = {
 414 0xc000,0xdb00,0x8748,0x9edc,
 415 0x3ffb,0xb947,0x266e,0x1a0f,
 416 0xbfcf,0x05b7,0x1f93,0x0b44,
 417 0x3f97,0xc8e1,0xcfca,0x4b28,
 418 0x3f2f,0x1a10,0xc720,0xff6f,
 419 0xbf15,0xdf0a,0x34d3,0x247d,
 420 0x3eec,0x97be,0xc7bc,0xd44a,
 421 0xbe98,0x87f3,0x5747,0xf1d8,
 422 0x3e4f,0x759a,0x33c7,0xc119,
 423 0x3e05,0x38d3,0x685d,0x5ea6,
 424 };
 425 static short B2[36] = {
 426   /* 0x3ff0,0x0000,0x0000,0x0000, */
 427 0xbfcd,0x6a4a,0x543a,0x3e24,
 428 0x3fbc,0x4a3b,0x9f44,0xe68b,
 429 0xbf90,0x0a99,0x3060,0xd9c5,
 430 0x3f66,0x9e75,0xf400,0x2e4f,
 431 0xbf2b,0x6204,0x7a54,0x8d9c,
 432 0x3ef1,0xea4c,0x2589,0xa8d9,
 433 0xbe9d,0xed0a,0x119c,0x2a52,
 434 0x3e4d,0xf316,0x2c36,0x7742,
 435 0x3e05,0x225b,0xe05c,0x9f6c,
 436 };
 437 #endif
 438
 439 /* x > 20
 440    x exp(-x) Ei(x) - 1  =  1/x A3(1/x)/B3(1/x)
 441    Theoretical absolute error = 6.15e-17  */
 442 #if UNK
 443 static double A3[9] = {
 444 -7.657847078286127362028E-1,
 445  6.886192415566705051750E-1,
 446 -2.132598113545206124553E-1,
 447  3.346107552384193813594E-2,
 448 -3.076541477344756050249E-3,
 449  1.747119316454907477380E-4,
 450 -6.103711682274170530369E-6,
 451  1.218032765428652199087E-7,
 452 -1.086076102793290233007E-9,
 453 };
 454 static double B3[9] = {
 455   /* 1.000000000000000000000E0, */
 456 -1.888802868662308731041E0,
 457  1.066691687211408896850E0,
 458 -2.751915982306380647738E-1,
 459  3.930852688233823569726E-2,
 460 -3.414684558602365085394E-3,
 461  1.866844370703555398195E-4,
 462 -6.345146083130515357861E-6,
 463  1.239754287483206878024E-7,
 464 -1.086076102793126632978E-9,
 465 };
 466 #endif
 467 #if DEC
 468 static short A3[36] = {
 469 0140104,0005167,0071746,0115510,
 470 0040060,0044531,0140741,0154556,
 471 0137532,0060307,0126506,0071123,
 472 0037011,0007173,0010405,0127224,
 473 0136111,0117715,0003654,0175577,
 474 0035067,0031340,0102657,0147714,
 475 0133714,0147173,0167473,0136640,
 476 0032402,0144407,0115547,0060114,
 477 0130625,0042347,0156431,0113425,
 478 };
 479 static short B3[36] = {
 480   /* 0040200,0000000,0000000,0000000, */
 481 0140361,0142112,0155277,0067714,
 482 0040210,0104532,0065676,0074326,
 483 0137614,0162751,0142421,0131033,
 484 0037041,0000772,0053236,0002632,
 485 0136137,0144346,0100536,0153136,
 486 0035103,0140270,0152211,0166215,
 487 0133724,0164143,0145763,0021153,
 488 0032405,0017033,0035333,0025736,
 489 0130625,0042347,0156431,0077134,
 490 };
 491 #endif
 492 #if IBMPC
 493 static short A3[36] = {
 494 0xd369,0xee7c,0x814e,0xbfe8,
 495 0x3b2e,0x383c,0x092b,0x3fe6,
 496 0xce4a,0xf5a8,0x4c18,0xbfcb,
 497 0xb5d2,0x6220,0x21cf,0x3fa1,
 498 0x9f70,0xa0f5,0x33f9,0xbf69,
 499 0xf9f9,0x10b5,0xe65c,0x3f26,
 500 0x77b4,0x7de7,0x99cf,0xbed9,
 501 0xec09,0xf36c,0x5920,0x3e80,
 502 0x32e3,0xfba3,0xa89c,0xbe12,
 503 };
 504 static short B3[36] = {
 505   /* 0x0000,0x0000,0x0000,0x3ff0, */
 506 0xedf9,0x5b57,0x3889,0xbffe,
 507 0xcf1b,0x4d77,0x112b,0x3ff1,
 508 0x3643,0x38a2,0x9cbd,0xbfd1,
 509 0xc0b3,0x4ad3,0x203f,0x3fa4,
 510 0xdacc,0xd02b,0xf91c,0xbf6b,
 511 0x3d92,0x1a91,0x7817,0x3f28,
 512 0x644d,0x797e,0x9d0c,0xbeda,
 513 0x657c,0x675b,0xa3c3,0x3e80,
 514 0x2fcb,0xfba3,0xa89c,0xbe12,
 515 };
 516 #endif
 517 #if MIEEE
 518 static short A3[36] = {
 519 0xbfe8,0x814e,0xee7c,0xd369,
 520 0x3fe6,0x092b,0x383c,0x3b2e,
 521 0xbfcb,0x4c18,0xf5a8,0xce4a,
 522 0x3fa1,0x21cf,0x6220,0xb5d2,
 523 0xbf69,0x33f9,0xa0f5,0x9f70,
 524 0x3f26,0xe65c,0x10b5,0xf9f9,
 525 0xbed9,0x99cf,0x7de7,0x77b4,
 526 0x3e80,0x5920,0xf36c,0xec09,
 527 0xbe12,0xa89c,0xfba3,0x32e3,
 528 };
 529 static short B3[36] = {
 530 /* 0x3ff0,0x0000,0x0000,0x0000, */
 531 0xbffe,0x3889,0x5b57,0xedf9,
 532 0x3ff1,0x112b,0x4d77,0xcf1b,
 533 0xbfd1,0x9cbd,0x38a2,0x3643,
 534 0x3fa4,0x203f,0x4ad3,0xc0b3,
 535 0xbf6b,0xf91c,0xd02b,0xdacc,
 536 0x3f28,0x7817,0x1a91,0x3d92,
 537 0xbeda,0x9d0c,0x797e,0x644d,
 538 0x3e80,0xa3c3,0x675b,0x657c,
 539 0xbe12,0xa89c,0xfba3,0x2fcb,
 540 };
 541 #endif
 542
 543 /* 16 <= x <= 32
 544    x exp(-x) Ei(x) - 1  =  1/x A4(1/x) / B4(1/x)
 545    Theoretical absolute error = 1.22e-17  */
 546 #if UNK
 547 static double A4[8] = {
 548 -2.458119367674020323359E-1,
 549 -1.483382253322077687183E-1,
 550  7.248291795735551591813E-2,
 551 -1.348315687380940523823E-2,
 552  1.342775069788636972294E-3,
 553 -7.942465637159712264564E-5,
 554  2.644179518984235952241E-6,
 555 -4.239473659313765177195E-8,
 556 };
 557 static double B4[8] = {
 558   /* 1.000000000000000000000E0, */
 559 -1.044225908443871106315E-1,
 560 -2.676453128101402655055E-1,
 561  9.695000254621984627876E-2,
 562 -1.601745692712991078208E-2,
 563  1.496414899205908021882E-3,
 564 -8.462452563778485013756E-5,
 565  2.728938403476726394024E-6,
 566 -4.239462431819542051337E-8,
 567 };
 568 #endif
 569 #if DEC
 570 static short A4[32] = {
 571 0137573,0133037,0152607,0113356,
 572 0137427,0162771,0145061,0126345,
 573 0037224,0070754,0110451,0174104,
 574 0136534,0164165,0072170,0063753,
 575 0035660,0000016,0002560,0147751,
 576 0134646,0110311,0123316,0047432,
 577 0033461,0071250,0101031,0075202,
 578 0132066,0012601,0077305,0170177,
 579 };
 580 static short B4[32] = {
 581   /* 0040200,0000000,0000000,0000000, */
 582 0137325,0155602,0162437,0030710,
 583 0137611,0004316,0071344,0176361,
 584 0037306,0106671,0011103,0155053,
 585 0136603,0033412,0132530,0175171,
 586 0035704,0021532,0015516,0166130,
 587 0134661,0074162,0036741,0073466,
 588 0033467,0021316,0003100,0171325,
 589 0132066,0012541,0162202,0150160,
 590 };
 591 #endif
 592 #if IBMPC
 593 static short A4[] = {
 594 0xf2de,0xfab0,0x76c3,0xbfcf,
 595 0x359d,0x3946,0xfcbf,0xbfc2,
 596 0x3f09,0x9225,0x8e3d,0x3fb2,
 597 0x0cfd,0xae8f,0x9d0e,0xbf8b,
 598 0x19fd,0xc0ae,0x0001,0x3f56,
 599 0xc9e3,0x34d9,0xd219,0xbf14,
 600 0x2f50,0x1043,0x2e55,0x3ec6,
 601 0xbe10,0x2fd8,0xc2b0,0xbe66,
 602 };
 603 static short B4[] = {
 604   /* 0x0000,0x0000,0x0000,0x3ff0, */
 605 0xe639,0x5ca3,0xbb70,0xbfba,
 606 0x9f9e,0xce5c,0x2119,0xbfd1,
 607 0x7b45,0x2248,0xd1b7,0x3fb8,
 608 0x1f4f,0x56ab,0x66e1,0xbf90,
 609 0xdd8b,0x4369,0x846b,0x3f58,
 610 0x2ee7,0x47bc,0x2f0e,0xbf16,
 611 0x1e5b,0xc0c8,0xe459,0x3ec6,
 612 0x5a0e,0x3c90,0xc2ac,0xbe66,
 613 };
 614 #endif
 615 #if MIEEE
 616 static short A4[32] = {
 617 0xbfcf,0x76c3,0xfab0,0xf2de,
 618 0xbfc2,0xfcbf,0x3946,0x359d,
 619 0x3fb2,0x8e3d,0x9225,0x3f09,
 620 0xbf8b,0x9d0e,0xae8f,0x0cfd,
 621 0x3f56,0x0001,0xc0ae,0x19fd,
 622 0xbf14,0xd219,0x34d9,0xc9e3,
 623 0x3ec6,0x2e55,0x1043,0x2f50,
 624 0xbe66,0xc2b0,0x2fd8,0xbe10,
 625 };
 626 static short B4[32] = {
 627   /* 0x3ff0,0x0000,0x0000,0x0000, */
 628 0xbfba,0xbb70,0x5ca3,0xe639,
 629 0xbfd1,0x2119,0xce5c,0x9f9e,
 630 0x3fb8,0xd1b7,0x2248,0x7b45,
 631 0xbf90,0x66e1,0x56ab,0x1f4f,
 632 0x3f58,0x846b,0x4369,0xdd8b,
 633 0xbf16,0x2f0e,0x47bc,0x2ee7,
 634 0x3ec6,0xe459,0xc0c8,0x1e5b,
 635 0xbe66,0xc2ac,0x3c90,0x5a0e,
 636 };
 637 #endif
 638
 639
 640 #if 0
 641 /* 20 <= x <= 40
 642    x exp(-x) Ei(x) - 1  =  1/x A4(1/x) / B4(1/x)
 643    Theoretical absolute error = 1.78e-17  */
 644 #if UNK
 645 static double A4[8] = {
 646  2.067245813525780707978E-1,
 647 -5.153749551345223645670E-1,
 648  1.928289589546695033096E-1,
 649 -3.124468842857260044075E-2,
 650  2.740283734277352539912E-3,
 651 -1.377775664366875175601E-4,
 652  3.803788980664744242323E-6,
 653 -4.611038277393688031154E-8,
 654 };
 655 static double B4[8] = {
 656   /*  1.000000000000000000000E0, */
 657 -8.544436025219516861531E-1,
 658  2.507436807692907385181E-1,
 659 -3.647688090228423114064E-2,
 660  3.008576950332041388892E-3,
 661 -1.452926405348421286334E-4,
 662  3.896007735260115431965E-6,
 663 -4.611037642697098234083E-8,
 664 };
 665 #endif
 666 #if DEC
 667 static short A4[32] = {
 668 0037523,0127633,0150301,0022031,
 669 0140003,0167634,0170572,0170420,
 670 0037505,0072364,0060672,0063220,
 671 0136777,0172334,0057456,0102640,
 672 0036063,0113125,0002476,0047251,
 673 0135020,0074142,0042600,0043630,
 674 0033577,0042230,0155372,0136105,
 675 0132106,0005346,0165333,0114541,
 676 };
 677 static short B4[28] = {
 678   /* 0040200,0000000,0000000,0000000, */
 679 0140132,0136320,0160433,0131535,
 680 0037600,0060571,0144452,0060214,
 681 0137025,0064310,0024220,0176472,
 682 0036105,0025613,0115762,0166605,
 683 0135030,0054662,0035454,0061763,
 684 0033602,0135163,0116430,0000066,
 685 0132106,0005345,0020602,0137133,
 686 };
 687 #endif
 688 #if IBMPC
 689 static short A4[32] = {
 690 0x2483,0x7a18,0x75f3,0x3fca,
 691 0x5e22,0x9e2f,0x7df3,0xbfe0,
 692 0x4cd2,0x8c37,0xae9e,0x3fc8,
 693 0xd0b4,0x8be5,0xfe9b,0xbf9f,
 694 0xc9d5,0xa0a7,0x72ca,0x3f66,
 695 0x08f3,0x48b0,0x0f0c,0xbf22,
 696 0x5789,0x1b5f,0xe893,0x3ecf,
 697 0x732c,0xdd5b,0xc15c,0xbe68,
 698 };
 699 static short B4[28] = {
 700   /* 0x0000,0x0000,0x0000,0x3ff0, */
 701 0x766c,0x1c23,0x579a,0xbfeb,
 702 0x4c11,0x3925,0x0c2f,0x3fd0,
 703 0x1fa7,0x0512,0xad19,0xbfa2,
 704 0x5db1,0x737e,0xa571,0x3f68,
 705 0x8c7e,0x4765,0x0b36,0xbf23,
 706 0x0007,0x73a3,0x574e,0x3ed0,
 707 0x57cb,0xa430,0xc15c,0xbe68,
 708 };
 709 #endif
 710 #if MIEEE
 711 static short A4[32] = {
 712 0x3fca,0x75f3,0x7a18,0x2483,
 713 0xbfe0,0x7df3,0x9e2f,0x5e22,
 714 0x3fc8,0xae9e,0x8c37,0x4cd2,
 715 0xbf9f,0xfe9b,0x8be5,0xd0b4,
 716 0x3f66,0x72ca,0xa0a7,0xc9d5,
 717 0xbf22,0x0f0c,0x48b0,0x08f3,
 718 0x3ecf,0xe893,0x1b5f,0x5789,
 719 0xbe68,0xc15c,0xdd5b,0x732c,
 720 };
 721 static short B4[28] = {
 722   /* 0x3ff0,0x0000,0x0000,0x0000, */
 723 0xbfeb,0x579a,0x1c23,0x766c,
 724 0x3fd0,0x0c2f,0x3925,0x4c11,
 725 0xbfa2,0xad19,0x0512,0x1fa7,
 726 0x3f68,0xa571,0x737e,0x5db1,
 727 0xbf23,0x0b36,0x4765,0x8c7e,
 728 0x3ed0,0x574e,0x73a3,0x0007,
 729 0xbe68,0xc15c,0xa430,0x57cb,
 730 };
 731 #endif
 732 #endif /* 0 */
 733
 734 /* 4 <= x <= 8
 735    x exp(-x) Ei(x) - 1  =  1/x A5(1/x) / B5(1/x)
 736    Theoretical absolute error = 2.20e-17  */
 737 #if UNK
 738 static double A5[8] = {
 739 -1.373215375871208729803E0,
 740 -7.084559133740838761406E-1,
 741  1.580806855547941010501E0,
 742 -2.601500427425622944234E-1,
 743  2.994674694113713763365E-2,
 744 -1.038086040188744005513E-3,
 745  4.371064420753005429514E-5,
 746  2.141783679522602903795E-6,
 747 };
 748 static double B5[8] = {
 749   /* 1.000000000000000000000E0, */
 750  8.585231423622028380768E-1,
 751  4.483285822873995129957E-1,
 752  7.687932158124475434091E-2,
 753  2.449868241021887685904E-2,
 754  8.832165941927796567926E-4,
 755  4.590952299511353531215E-4,
 756 -4.729848351866523044863E-6,
 757  2.665195537390710170105E-6,
 758 };
 759 #endif
 760 #if DEC
 761 static short A5[32] = {
 762 0140257,0142605,0076335,0113632,
 763 0140065,0056535,0161231,0074311,
 764 0040312,0053741,0004357,0076405,
 765 0137605,0031142,0165503,0136705,
 766 0036765,0051341,0053573,0007602,
 767 0135610,0010143,0027643,0110522,
 768 0034467,0052762,0062024,0120161,
 769 0033417,0135620,0036500,0062647,
 770 };
 771 static short B[32] = {
 772   /* 0040200,0000000,0000000,0000000, */
 773 0040133,0144054,0031516,0004100,
 774 0037745,0105522,0166622,0123146,
 775 0037235,0071347,0157560,0157464,
 776 0036710,0130565,0173747,0041670,
 777 0035547,0103651,0106243,0101240,
 778 0035360,0131267,0176263,0140257,
 779 0133636,0132426,0102537,0102531,
 780 0033462,0155665,0167503,0176350,
 781 };
 782 #endif
 783 #if IBMPC
 784 static short A5[32] = {
 785 0xb2f3,0xaf9b,0xf8b0,0xbff5,
 786 0x2f19,0xbc53,0xabab,0xbfe6,
 787 0xefa1,0x211d,0x4afc,0x3ff9,
 788 0x77b9,0x5d68,0xa64c,0xbfd0,
 789 0x61f0,0x2aef,0xaa5c,0x3f9e,
 790 0x722a,0x65f4,0x020c,0xbf51,
 791 0x940e,0x4c82,0xeabe,0x3f06,
 792 0x0cb5,0x07a8,0xf772,0x3ec1,
 793 };
 794 static short B5[32] = {
 795   /* 0x0000,0x0000,0x0000,0x3ff0, */
 796 0xc108,0x8669,0x7905,0x3feb,
 797 0x54cd,0x5db2,0xb16a,0x3fdc,
 798 0x1be7,0xfbee,0xae5c,0x3fb3,
 799 0xe877,0xbefc,0x162e,0x3f99,
 800 0x7054,0x3194,0xf0f5,0x3f4c,
 801 0x7816,0xff96,0x1656,0x3f3e,
 802 0xf0ab,0xd0ab,0xd6a2,0xbed3,
 803 0x7f9d,0xbde8,0x5b76,0x3ec6,
 804 };
 805 #endif
 806 #if MIEEE
 807 static short A5[32] = {
 808 0xbff5,0xf8b0,0xaf9b,0xb2f3,
 809 0xbfe6,0xabab,0xbc53,0x2f19,
 810 0x3ff9,0x4afc,0x211d,0xefa1,
 811 0xbfd0,0xa64c,0x5d68,0x77b9,
 812 0x3f9e,0xaa5c,0x2aef,0x61f0,
 813 0xbf51,0x020c,0x65f4,0x722a,
 814 0x3f06,0xeabe,0x4c82,0x940e,
 815 0x3ec1,0xf772,0x07a8,0x0cb5,
 816 };
 817 static short B5[32] = {
 818   /* 0x3ff0,0x0000,0x0000,0x0000, */
 819 0x3feb,0x7905,0x8669,0xc108,
 820 0x3fdc,0xb16a,0x5db2,0x54cd,
 821 0x3fb3,0xae5c,0xfbee,0x1be7,
 822 0x3f99,0x162e,0xbefc,0xe877,
 823 0x3f4c,0xf0f5,0x3194,0x7054,
 824 0x3f3e,0x1656,0xff96,0x7816,
 825 0xbed3,0xd6a2,0xd0ab,0xf0ab,
 826 0x3ec6,0x5b76,0xbde8,0x7f9d,
 827 };
 828 #endif
 829 /* 2 <= x <= 4
 830    x exp(-x) Ei(x) - 1  =  1/x A6(1/x) / B6(1/x)
 831    Theoretical absolute error = 4.89e-17  */
 832 #if UNK
 833 static double A6[8] = {
 834  1.981808503259689673238E-2,
 835 -1.271645625984917501326E0,
 836 -2.088160335681228318920E0,
 837  2.755544509187936721172E0,
 838 -4.409507048701600257171E-1,
 839  4.665623805935891391017E-2,
 840 -1.545042679673485262580E-3,
 841  7.059980605299617478514E-5,
 842 };
 843 static double B6[7] = {
 844   /* 1.000000000000000000000E0, */
 845  1.476498670914921440652E0,
 846  5.629177174822436244827E-1,
 847  1.699017897879307263248E-1,
 848  2.291647179034212017463E-2,
 849  4.450150439728752875043E-3,
 850  1.727439612206521482874E-4,
 851  3.953167195549672482304E-5,
 852 };
 853 #endif
 854 #if DEC
 855 static short A6[32] = {
 856 0036642,0054611,0061263,0000140,
 857 0140242,0142510,0125732,0072035,
 858 0140405,0122153,0037643,0104527,
 859 0040460,0055327,0055550,0116240,
 860 0137741,0142112,0070441,0103510,
 861 0037077,0015234,0104750,0146765,
 862 0135712,0101407,0107554,0020253,
 863 0034624,0007373,0072621,0063735,
 864 };
 865 static short B6[28] = {
 866   /* 0040200,0000000,0000000,0000000, */
 867 0040274,0176750,0110025,0061006,
 868 0040020,0015540,0021354,0155050,
 869 0037455,0175274,0015257,0021112,
 870 0036673,0135523,0016042,0117203,
 871 0036221,0151221,0046352,0144174,
 872 0035065,0021232,0117727,0152432,
 873 0034445,0147317,0037300,0067123,
 874 };
 875 #endif
 876 #if IBMPC
 877 static short A6[32] = {
 878 0x600c,0x2c56,0x4b31,0x3f94,
 879 0x4e84,0x157b,0x58a9,0xbff4,
 880 0x712b,0x67f4,0xb48d,0xc000,
 881 0x1394,0xeb6d,0x0b5a,0x4006,
 882 0x30e9,0x4e24,0x3889,0xbfdc,
 883 0x19bf,0x913d,0xe353,0x3fa7,
 884 0x8415,0xf1ed,0x5060,0xbf59,
 885 0x2cfc,0x6eb2,0x81df,0x3f12,
 886 };
 887 static short B6[28] = {
 888   /* 0x0000,0x0000,0x0000,0x3ff0, */
 889 0xac41,0x1202,0x9fbd,0x3ff7,
 890 0x9b45,0x045d,0x036c,0x3fe2,
 891 0xe449,0x8355,0xbf57,0x3fc5,
 892 0x53d0,0x6384,0x776a,0x3f97,
 893 0x590f,0x299d,0x3a52,0x3f72,
 894 0xfaa3,0x53fa,0xa453,0x3f26,
 895 0x0dca,0xe7d8,0xb9d9,0x3f04,
 896 };
 897 #endif
 898 #if MIEEE
 899 static short A6[32] = {
 900 0x3f94,0x4b31,0x2c56,0x600c,
 901 0xbff4,0x58a9,0x157b,0x4e84,
 902 0xc000,0xb48d,0x67f4,0x712b,
 903 0x4006,0x0b5a,0xeb6d,0x1394,
 904 0xbfdc,0x3889,0x4e24,0x30e9,
 905 0x3fa7,0xe353,0x913d,0x19bf,
 906 0xbf59,0x5060,0xf1ed,0x8415,
 907 0x3f12,0x81df,0x6eb2,0x2cfc,
 908 };
 909 static short B6[28] = {
 910   /* 0x3ff0,0x0000,0x0000,0x0000, */
 911 0x3ff7,0x9fbd,0x1202,0xac41,
 912 0x3fe2,0x036c,0x045d,0x9b45,
 913 0x3fc5,0xbf57,0x8355,0xe449,
 914 0x3f97,0x776a,0x6384,0x53d0,
 915 0x3f72,0x3a52,0x299d,0x590f,
 916 0x3f26,0xa453,0x53fa,0xfaa3,
 917 0x3f04,0xb9d9,0xe7d8,0x0dca,
 918 };
 919 #endif
 920 /* 32 <= x <= 64
 921    x exp(-x) Ei(x) - 1  =  1/x A7(1/x) / B7(1/x)
 922    Theoretical absolute error = 7.71e-18  */
 923 #if UNK
 924 static double A7[6] = {
 925  1.212561118105456670844E-1,
 926 -5.823133179043894485122E-1,
 927  2.348887314557016779211E-1,
 928 -3.040034318113248237280E-2,
 929  1.510082146865190661777E-3,
 930 -2.523137095499571377122E-5,
 931 };
 932 static double B7[5] = {
 933   /* 1.000000000000000000000E0, */
 934 -1.002252150365854016662E0,
 935  2.928709694872224144953E-1,
 936 -3.337004338674007801307E-2,
 937  1.560544881127388842819E-3,
 938 -2.523137093603234562648E-5,
 939 };
 940 #endif
 941 #if DEC
 942 static short A7[24] = {
 943 0037370,0052437,0152524,0150125,
 944 0140025,0011174,0050154,0131330,
 945 0037560,0103253,0167464,0062245,
 946 0136771,0005043,0174001,0023345,
 947 0035705,0166762,0157300,0016451,
 948 0134323,0123764,0157767,0134477,
 949 };
 950 static short B7[20] = {
 951   /* 0040200,0000000,0000000,0000000, */
 952 0140200,0044714,0064025,0060324,
 953 0037625,0171457,0003712,0073131,
 954 0137010,0127406,0150061,0141746,
 955 0035714,0105462,0072356,0103712,
 956 0134323,0123764,0156514,0077414,
 957 };
 958 #endif
 959 #if IBMPC
 960 static short A7[24] = {
 961 0x9a0b,0xfaaa,0x0aa3,0x3fbf,
 962 0x965b,0x8a0d,0xa24f,0xbfe2,
 963 0x8c95,0x7de6,0x10d5,0x3fce,
 964 0x24dd,0x7f00,0x2144,0xbf9f,
 965 0x03a5,0x5bd8,0xbdbe,0x3f58,
 966 0xf728,0x9bfe,0x74fe,0xbefa,
 967 };
 968 static short B7[20] = {
 969   /* 0x0000,0x0000,0x0000,0x3ff0, */
 970 0xac1a,0x8d02,0x0939,0xbff0,
 971 0x4ecb,0xe0f9,0xbe65,0x3fd2,
 972 0x387d,0xda06,0x15e0,0xbfa1,
 973 0xd0f9,0x4e9d,0x9166,0x3f59,
 974 0x8fe2,0x9ba9,0x74fe,0xbefa,
 975 };
 976 #endif
 977 #if MIEEE
 978 static short A7[24] = {
 979 0x3fbf,0x0aa3,0xfaaa,0x9a0b,
 980 0xbfe2,0xa24f,0x8a0d,0x965b,
 981 0x3fce,0x10d5,0x7de6,0x8c95,
 982 0xbf9f,0x2144,0x7f00,0x24dd,
 983 0x3f58,0xbdbe,0x5bd8,0x03a5,
 984 0xbefa,0x74fe,0x9bfe,0xf728,
 985 };
 986 static short B7[20] = {
 987   /* 0x3ff0,0x0000,0x0000,0x0000, */
 988 0xbff0,0x0939,0x8d02,0xac1a,
 989 0x3fd2,0xbe65,0xe0f9,0x4ecb,
 990 0xbfa1,0x15e0,0xda06,0x387d,
 991 0x3f59,0x9166,0x4e9d,0xd0f9,
 992 0xbefa,0x74fe,0x9ba9,0x8fe2,
 993 };
 994 #endif
 995
 996 double ei (x)
 997 double x;
 998 {
 999   double f, w;
1000
1001   if (x <= 0.0)
1002     {
1003       mtherr("ei", DOMAIN);
1004       return 0.0;
1005     }
1006   else if (x < 2.0)
1007     {
1008   /* Power series.
1009                             inf    n
1010                              -    x
1011      Ei(x) = EUL + ln x  +   >   ----
1012                              -   n n!
1013                             n=1
1014   */
1015       f = polevl(x,A,5) / p1evl(x,B,6);
1016       /*      f = polevl(x,A,6) / p1evl(x,B,7); */
1017       /*      f = polevl(x,A,8) / p1evl(x,B,9); */
1018       return (EUL + log(x) + x * f);
1019     }
1020   else if (x < 4.0)
1021     {
1022   /* Asymptotic expansion.
1023                             1       2       6
1024     x exp(-x) Ei(x) =  1 + ---  +  ---  +  ---- + ...
1025                             x        2       3
1026                                     x       x
1027   */
1028       w = 1.0/x;
1029       f = polevl(w,A6,7) / p1evl(w,B6,7);
1030       return (exp(x) * w * (1.0 + w * f));
1031     }
1032   else if (x < 8.0)
1033     {
1034       w = 1.0/x;
1035       f = polevl(w,A5,7) / p1evl(w,B5,8);
1036       return (exp(x) * w * (1.0 + w * f));
1037     }
1038   else if (x < 16.0)
1039     {
1040       w = 1.0/x;
1041       f = polevl(w,A2,9) / p1evl(w,B2,9);
1042       return (exp(x) * w * (1.0 + w * f));
1043     }
1044   else if (x < 32.0)
1045     {
1046       w = 1.0/x;
1047       f = polevl(w,A4,7) / p1evl(w,B4,8);
1048       return (exp(x) * w * (1.0 + w * f));
1049     }
1050   else if (x < 64.0)
1051     {
1052       w = 1.0/x;
1053       f = polevl(w,A7,5) / p1evl(w,B7,5);
1054       return (exp(x) * w * (1.0 + w * f));
1055     }
1056   else
1057     {
1058       w = 1.0/x;
1059       f = polevl(w,A3,8) / p1evl(w,B3,9);
1060       return (exp(x) * w * (1.0 + w * f));
1061     }
1062 }