Quelle  ei.c

/* ei.c
 *
 * Exponential integral
 *
 *
 * SYNOPSIS:
 *
 * double x, y, ei();
 *
 * y = ei( x );
 *
 *
 *
 * DESCRIPTION:
 *
 *               x
 *                -     t
 *               | |   e
 *    Ei(x) =   -|-   ---  dt .
 *             | |     t
 *              -
 *             -inf
 *
 * Not defined for x <= 0.
 * See also expn.c.
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE       0,100       50000      8.6e-16     1.3e-16
 *
 */

/*
Cephes Math Library Release 2.8:  May, 1999
Copyright 1999 by Stephen L. Moshier
*/

#include "mconf.h"
#ifdef ANSIPROT
extern double log(double);
extern double exp(double);
extern double polevl(double, void *, int);
extern double p1evl(double, void *, int);
#else
extern double log(), exp(), polevl(), p1evl();
#endif

#define EUL 5.772156649015328606065e-1

/* 0 < x <= 2
   Ei(x) - EUL - ln(x) = x A(x)/B(x)
   Theoretical peak relative error 9.73e-18  */
#if UNK
static double A[6] = {
    -5.350447357812542947283E0, 2.185049168816613393830E2,
    -4.176572384826693777058E3, 5.541176756393557601232E4,
    -3.313381331178144034309E5, 1.592627163384945414220E6,
};
static double B[6] = {
    /*  1.000000000000000000000E0, */
    -5.250547959112862969197E1, 1.259616186786790571525E3,
    -1.756549581973534652631E4, 1.493062117002725991967E5,
    -7.294949239640527645655E5, 1.592627163384945429726E6,
};
#endif
#if DEC
static short A[24] = {
    0140653, 0033335, 0060230, 0144217, 0042132, 0100502, 0035625, 0167413,
    0143202, 0102224, 0037176, 0175403, 0044130, 0071704, 0077421, 0170343,
    0144641, 0144504, 0041200, 0045154, 0045302, 0064631, 0047234, 0142052,
};
static short B[24] = {
    /* 0040200,0000000,0000000,0000000, */
    0141522, 0002634, 0070442, 0142614, 0042635, 0071667, 0146532, 0027705,
    0143611, 0035375, 0156025, 0114015, 0044421, 0147215, 0106177, 0046330,
    0145062, 0014556, 0144216, 0103725, 0045302, 0064631, 0047234, 0142052,
};
#endif
#if IBMPC
static short A[24] = {
    0x1912, 0xac13, 0x66db, 0xc015, 0xbde1, 0x4772, 0x5028, 0x406b,
    0xdf60, 0x87cf, 0x5092, 0xc0b0, 0x3e1c, 0x8fe2, 0x0e78, 0x40eb,
    0x094e, 0x8850, 0x3928, 0xc114, 0x9885, 0x29d3, 0x4d33, 0x4138,
};
static short B[24] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0x58b1, 0x8e24, 0x40b3, 0xc04a, 0x45f9, 0xf9ab, 0xae76, 0x4093,
    0xb302, 0xbb82, 0x275f, 0xc0d1, 0xe99b, 0xb18f, 0x39d1, 0x4102,
    0xd0fb, 0xd911, 0x432d, 0xc126, 0x9885, 0x29d3, 0x4d33, 0x4138,
};
#endif
#if MIEEE
static short A[24] = {
    0xc015, 0x66db, 0xac13, 0x1912, 0x406b, 0x5028, 0x4772, 0xbde1,
    0xc0b0, 0x5092, 0x87cf, 0xdf60, 0x40eb, 0x0e78, 0x8fe2, 0x3e1c,
    0xc114, 0x3928, 0x8850, 0x094e, 0x4138, 0x4d33, 0x29d3, 0x9885,
};
static short B[24] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0xc04a, 0x40b3, 0x8e24, 0x58b1, 0x4093, 0xae76, 0xf9ab, 0x45f9,
    0xc0d1, 0x275f, 0xbb82, 0xb302, 0x4102, 0x39d1, 0xb18f, 0xe99b,
    0xc126, 0x432d, 0xd911, 0xd0fb, 0x4138, 0x4d33, 0x29d3, 0x9885,
};
#endif

#if 0
/* 0 < x <= 4
   Ei(x) - EUL - ln(x) = x A(x)/B(x)
   Theoretical peak relative error 4.75e-17  */
#if UNK
static double A[7] = {
-6.831869820732773831942E0,
 2.920190530726774500309E2,
-1.195883839286649567993E4,
 1.761045255472548975666E5,
-2.623034438354006526979E6,
 1.472430336917880803157E7,
-8.205359388213261174960E7,
};
static double B[7] = {
  /* 1.000000000000000000000E0, */
-7.731946237840033971071E1,
 2.751808700543578450827E3,
-5.829268609072186897994E4,
 7.916610857961870631379E5,
-6.873926904825733094076E6,
 3.523770183971164032710E7,
-8.205359388213260785363E7,
};
#endif
#if DEC
static short A[28] = {
0140732,0117255,0072522,0071743,
0042222,0001160,0052302,0002334,
0143472,0155532,0101650,0155462,
0044453,0175041,0121220,0172022,
0145440,0014351,0140337,0157550,
0046140,0126317,0057202,0100233,
0146634,0100473,0036072,0067054,
};
static short B[28] = {
  /* 0040200,0000000,0000000,0000000, */
0141632,0121620,0111247,0010115,
0043053,0176360,0067773,0027324,
0144143,0132257,0121644,0036204,
0045101,0043321,0057553,0151231,
0145721,0143215,0147505,0050610,
0046406,0065721,0072675,0152744,
0146634,0100473,0036072,0067052,
};
#endif
#if IBMPC
static short A[28] = {
0x4e7c,0xaeaa,0x53d5,0xc01b,
0x409b,0x0a98,0x404e,0x4072,
0x1b66,0x5075,0x5b6b,0xc0c7,
0x1e82,0x3452,0x7f44,0x4105,
0xfbed,0x381b,0x031d,0xc144,
0x5013,0xebd0,0x1599,0x416c,
0x4dc5,0x6787,0x9027,0xc193,
};
static short B[28] = {
  /* 0x0000,0x0000,0x0000,0x3ff0, */
0xe20a,0x1254,0x5472,0xc053,
0x65db,0x0dff,0x7f9e,0x40a5,
0x8791,0xf474,0x7695,0xc0ec,
0x7a53,0x2bed,0x28da,0x4128,
0xaa31,0xb9e8,0x38d1,0xc15a,
0xbabd,0x2eb7,0xcd7a,0x4180,
0x4dc5,0x6787,0x9027,0xc193,
};
#endif
#if MIEEE
static short A[28] = {
0xc01b,0x53d5,0xaeaa,0x4e7c,
0x4072,0x404e,0x0a98,0x409b,
0xc0c7,0x5b6b,0x5075,0x1b66,
0x4105,0x7f44,0x3452,0x1e82,
0xc144,0x031d,0x381b,0xfbed,
0x416c,0x1599,0xebd0,0x5013,
0xc193,0x9027,0x6787,0x4dc5,
};
static short B[28] = {
  /* 0x3ff0,0x0000,0x0000,0x0000, */
0xc053,0x5472,0x1254,0xe20a,
0x40a5,0x7f9e,0x0dff,0x65db,
0xc0ec,0x7695,0xf474,0x8791,
0x4128,0x28da,0x2bed,0x7a53,
0xc15a,0x38d1,0xb9e8,0xaa31,
0x4180,0xcd7a,0x2eb7,0xbabd,
0xc193,0x9027,0x6787,0x4dc5,
};
#endif
#endif /* 0 */

#if 0
/* 0 < x <= 8
   Ei(x) - EUL - ln(x) = x A(x)/B(x)
   Theoretical peak relative error 2.14e-17  */

#if UNK
static double A[9] = {
-1.111230942210860450145E1,
 3.688203982071386319616E2,
-4.924786153494029574350E4,
 1.050677503345557903241E6,
-3.626713709916703688968E7,
 4.353499908839918635414E8,
-6.454613717232006895409E9,
 3.408243056457762907071E10,
-1.995466674647028468613E11,
};
static double B[9] = {
  /*  1.000000000000000000000E0, */
-1.356757648138514017969E2,
 8.562181317107341736606E3,
-3.298257180413775117555E5,
 8.543534058481435917210E6,
-1.542380618535140055068E8,
 1.939251779195993632028E9,
-1.636096210465615015435E10,
 8.396909743075306970605E10,
-1.995466674647028425886E11,
};
#endif
#if DEC
static short A[36] = {
0141061,0146004,0173357,0151553,
0042270,0064402,0147366,0126701,
0144100,0057734,0106615,0144356,
0045200,0040654,0003332,0004456,
0146412,0054440,0043130,0140263,
0047317,0113517,0033422,0065123,
0150300,0056313,0065235,0131147,
0050775,0167423,0146222,0075760,
0151471,0153642,0003442,0147667,
};
static short B[36] = {
  /* 0040200,0000000,0000000,0000000, */
0142007,0126376,0166077,0043600,
0043405,0144271,0125461,0014364,
0144641,0006066,0175061,0164463,
0046002,0056456,0007370,0121657,
0147023,0013706,0156647,0177115,
0047747,0026504,0103144,0054507,
0150563,0146036,0007051,0177135,
0051234,0063625,0173266,0003111,
0151471,0153642,0003442,0147666,
};
#endif
#if IBMPC
static short A[36] = {
0xfa6d,0x9edd,0x3980,0xc026,
0xd5b8,0x59de,0x0d20,0x4077,
0xb91e,0x91b1,0x0bfb,0xc0e8,
0x4126,0x80db,0x0835,0x4130,
0x1816,0x08cb,0x4b24,0xc181,
0x4d4a,0xe6e2,0xf2e9,0x41b9,
0xb64d,0x6d53,0x0b99,0xc1f8,
0x4f7e,0x7992,0xbde2,0x421f,
0x59f7,0x40e4,0x3af4,0xc247,
};
static short B[36] = {
  /* 0x0000,0x0000,0x0000,0x3ff0, */
0xe8f0,0xdd87,0xf59f,0xc060,
0x231e,0x3566,0xb917,0x40c0,
0x3d26,0xdf46,0x2186,0xc114,
0x1476,0xc1df,0x4ba5,0x4160,
0xffca,0xdbb4,0x62f8,0xc1a2,
0x8b29,0x90cc,0xe5a8,0x41dc,
0x3fcc,0xc1c5,0x7983,0xc20e,
0xc0c9,0xbed6,0x8cf2,0x4233,
0x59f7,0x40e4,0x3af4,0xc247,
};
#endif
#if MIEEE
static short A[36] = {
0xc026,0x3980,0x9edd,0xfa6d,
0x4077,0x0d20,0x59de,0xd5b8,
0xc0e8,0x0bfb,0x91b1,0xb91e,
0x4130,0x0835,0x80db,0x4126,
0xc181,0x4b24,0x08cb,0x1816,
0x41b9,0xf2e9,0xe6e2,0x4d4a,
0xc1f8,0x0b99,0x6d53,0xb64d,
0x421f,0xbde2,0x7992,0x4f7e,
0xc247,0x3af4,0x40e4,0x59f7,
};
static short B[36] = {
  /* 0x3ff0,0x0000,0x0000,0x0000, */
0xc060,0xf59f,0xdd87,0xe8f0,
0x40c0,0xb917,0x3566,0x231e,
0xc114,0x2186,0xdf46,0x3d26,
0x4160,0x4ba5,0xc1df,0x1476,
0xc1a2,0x62f8,0xdbb4,0xffca,
0x41dc,0xe5a8,0x90cc,0x8b29,
0xc20e,0x7983,0xc1c5,0x3fcc,
0x4233,0x8cf2,0xbed6,0xc0c9,
0xc247,0x3af4,0x40e4,0x59f7,
};
#endif
#endif /* 0 */

/* 8 <= x <= 20
   x exp(-x) Ei(x) - 1 = 1/x R(1/x)
   Theoretical peak absolute error = 1.07e-17  */
#if UNK
static double A2[10] = {
    -2.106934601691916512584E0,  1.732733869664688041885E0,
    -2.423619178935841904839E-1, 2.322724180937565842585E-2,
    2.372880440493179832059E-4,  -8.343219561192552752335E-5,
    1.363408795605250394881E-5,  -3.655412321999253963714E-7,
    1.464941733975961318456E-8,  6.176407863710360207074E-10,
};
static double B2[9] = {
    /* 1.000000000000000000000E0, */
    -2.298062239901678075778E-1, 1.105077041474037862347E-1,
    -1.566542966630792353556E-2, 2.761106850817352773874E-3,
    -2.089148012284048449115E-4, 1.708528938807675304186E-5,
    -4.459311796356686423199E-7, 1.394634930353847498145E-8,
    6.150865933977338354138E-10,
};
#endif
#if DEC
static short A2[40] = {
    0140406, 0154004, 0035104, 0173336, 0040335, 0145071, 0031560, 0150165,
    0137570, 0026670, 0176230, 0055040, 0036676, 0043416, 0077122, 0054476,
    0035170, 0150206, 0034407, 0175571, 0134656, 0174121, 0123231, 0021751,
    0034144, 0136766, 0036746, 0121115, 0132704, 0037632, 0135077, 0107300,
    0031573, 0126321, 0117076, 0004314, 0030451, 0143233, 0041352, 0172464,
};
static short B2[36] = {
    /* 0040200,0000000,0000000,0000000, */
    0137553, 0051122, 0120721, 0170437, 0037342, 0050734, 0175047, 0032132,
    0136600, 0052311, 0101406, 0147050, 0036064, 0171657, 0120001, 0071165,
    0135133, 0010043, 0151244, 0066340, 0034217, 0051141, 0026115, 0043305,
    0132757, 0064120, 0106341, 0051217, 0031557, 0114261, 0060663, 0135017,
    0030451, 0011337, 0001344, 0175542,
};
#endif
#if IBMPC
static short A2[40] = {
    0x9edc, 0x8748, 0xdb00, 0xc000, 0x1a0f, 0x266e, 0xb947, 0x3ffb,
    0x0b44, 0x1f93, 0x05b7, 0xbfcf, 0x4b28, 0xcfca, 0xc8e1, 0x3f97,
    0xff6f, 0xc720, 0x1a10, 0x3f2f, 0x247d, 0x34d3, 0xdf0a, 0xbf15,
    0xd44a, 0xc7bc, 0x97be, 0x3eec, 0xf1d8, 0x5747, 0x87f3, 0xbe98,
    0xc119, 0x33c7, 0x759a, 0x3e4f, 0x5ea6, 0x685d, 0x38d3, 0x3e05,
};
static short B2[36] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0x3e24, 0x543a, 0x6a4a, 0xbfcd, 0xe68b, 0x9f44, 0x4a3b, 0x3fbc, 0xd9c5,
    0x3060, 0x0a99, 0xbf90, 0x2e4f, 0xf400, 0x9e75, 0x3f66, 0x8d9c, 0x7a54,
    0x6204, 0xbf2b, 0xa8d9, 0x2589, 0xea4c, 0x3ef1, 0x2a52, 0x119c, 0xed0a,
    0xbe9d, 0x7742, 0x2c36, 0xf316, 0x3e4d, 0x9f6c, 0xe05c, 0x225b, 0x3e05,
};
#endif
#if MIEEE
static short A2[40] = {
    0xc000, 0xdb00, 0x8748, 0x9edc, 0x3ffb, 0xb947, 0x266e, 0x1a0f,
    0xbfcf, 0x05b7, 0x1f93, 0x0b44, 0x3f97, 0xc8e1, 0xcfca, 0x4b28,
    0x3f2f, 0x1a10, 0xc720, 0xff6f, 0xbf15, 0xdf0a, 0x34d3, 0x247d,
    0x3eec, 0x97be, 0xc7bc, 0xd44a, 0xbe98, 0x87f3, 0x5747, 0xf1d8,
    0x3e4f, 0x759a, 0x33c7, 0xc119, 0x3e05, 0x38d3, 0x685d, 0x5ea6,
};
static short B2[36] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0xbfcd, 0x6a4a, 0x543a, 0x3e24, 0x3fbc, 0x4a3b, 0x9f44, 0xe68b, 0xbf90,
    0x0a99, 0x3060, 0xd9c5, 0x3f66, 0x9e75, 0xf400, 0x2e4f, 0xbf2b, 0x6204,
    0x7a54, 0x8d9c, 0x3ef1, 0xea4c, 0x2589, 0xa8d9, 0xbe9d, 0xed0a, 0x119c,
    0x2a52, 0x3e4d, 0xf316, 0x2c36, 0x7742, 0x3e05, 0x225b, 0xe05c, 0x9f6c,
};
#endif

/* x > 20
   x exp(-x) Ei(x) - 1  =  1/x A3(1/x)/B3(1/x)
   Theoretical absolute error = 6.15e-17  */
#if UNK
static double A3[9] = {
    -7.657847078286127362028E-1, 6.886192415566705051750E-1,
    -2.132598113545206124553E-1, 3.346107552384193813594E-2,
    -3.076541477344756050249E-3, 1.747119316454907477380E-4,
    -6.103711682274170530369E-6, 1.218032765428652199087E-7,
    -1.086076102793290233007E-9,
};
static double B3[9] = {
    /* 1.000000000000000000000E0, */
    -1.888802868662308731041E0,  1.066691687211408896850E0,
    -2.751915982306380647738E-1, 3.930852688233823569726E-2,
    -3.414684558602365085394E-3, 1.866844370703555398195E-4,
    -6.345146083130515357861E-6, 1.239754287483206878024E-7,
    -1.086076102793126632978E-9,
};
#endif
#if DEC
static short A3[36] = {
    0140104, 0005167, 0071746, 0115510, 0040060, 0044531, 0140741, 0154556,
    0137532, 0060307, 0126506, 0071123, 0037011, 0007173, 0010405, 0127224,
    0136111, 0117715, 0003654, 0175577, 0035067, 0031340, 0102657, 0147714,
    0133714, 0147173, 0167473, 0136640, 0032402, 0144407, 0115547, 0060114,
    0130625, 0042347, 0156431, 0113425,
};
static short B3[36] = {
    /* 0040200,0000000,0000000,0000000, */
    0140361, 0142112, 0155277, 0067714, 0040210, 0104532, 0065676, 0074326,
    0137614, 0162751, 0142421, 0131033, 0037041, 0000772, 0053236, 0002632,
    0136137, 0144346, 0100536, 0153136, 0035103, 0140270, 0152211, 0166215,
    0133724, 0164143, 0145763, 0021153, 0032405, 0017033, 0035333, 0025736,
    0130625, 0042347, 0156431, 0077134,
};
#endif
#if IBMPC
static short A3[36] = {
    0xd369, 0xee7c, 0x814e, 0xbfe8, 0x3b2e, 0x383c, 0x092b, 0x3fe6, 0xce4a,
    0xf5a8, 0x4c18, 0xbfcb, 0xb5d2, 0x6220, 0x21cf, 0x3fa1, 0x9f70, 0xa0f5,
    0x33f9, 0xbf69, 0xf9f9, 0x10b5, 0xe65c, 0x3f26, 0x77b4, 0x7de7, 0x99cf,
    0xbed9, 0xec09, 0xf36c, 0x5920, 0x3e80, 0x32e3, 0xfba3, 0xa89c, 0xbe12,
};
static short B3[36] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0xedf9, 0x5b57, 0x3889, 0xbffe, 0xcf1b, 0x4d77, 0x112b, 0x3ff1, 0x3643,
    0x38a2, 0x9cbd, 0xbfd1, 0xc0b3, 0x4ad3, 0x203f, 0x3fa4, 0xdacc, 0xd02b,
    0xf91c, 0xbf6b, 0x3d92, 0x1a91, 0x7817, 0x3f28, 0x644d, 0x797e, 0x9d0c,
    0xbeda, 0x657c, 0x675b, 0xa3c3, 0x3e80, 0x2fcb, 0xfba3, 0xa89c, 0xbe12,
};
#endif
#if MIEEE
static short A3[36] = {
    0xbfe8, 0x814e, 0xee7c, 0xd369, 0x3fe6, 0x092b, 0x383c, 0x3b2e, 0xbfcb,
    0x4c18, 0xf5a8, 0xce4a, 0x3fa1, 0x21cf, 0x6220, 0xb5d2, 0xbf69, 0x33f9,
    0xa0f5, 0x9f70, 0x3f26, 0xe65c, 0x10b5, 0xf9f9, 0xbed9, 0x99cf, 0x7de7,
    0x77b4, 0x3e80, 0x5920, 0xf36c, 0xec09, 0xbe12, 0xa89c, 0xfba3, 0x32e3,
};
static short B3[36] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0xbffe, 0x3889, 0x5b57, 0xedf9, 0x3ff1, 0x112b, 0x4d77, 0xcf1b, 0xbfd1,
    0x9cbd, 0x38a2, 0x3643, 0x3fa4, 0x203f, 0x4ad3, 0xc0b3, 0xbf6b, 0xf91c,
    0xd02b, 0xdacc, 0x3f28, 0x7817, 0x1a91, 0x3d92, 0xbeda, 0x9d0c, 0x797e,
    0x644d, 0x3e80, 0xa3c3, 0x675b, 0x657c, 0xbe12, 0xa89c, 0xfba3, 0x2fcb,
};
#endif

/* 16 <= x <= 32
   x exp(-x) Ei(x) - 1  =  1/x A4(1/x) / B4(1/x)
   Theoretical absolute error = 1.22e-17  */
#if UNK
static double A4[8] = {
    -2.458119367674020323359E-1, -1.483382253322077687183E-1,
    7.248291795735551591813E-2,  -1.348315687380940523823E-2,
    1.342775069788636972294E-3,  -7.942465637159712264564E-5,
    2.644179518984235952241E-6,  -4.239473659313765177195E-8,
};
static double B4[8] = {
    /* 1.000000000000000000000E0, */
    -1.044225908443871106315E-1, -2.676453128101402655055E-1,
    9.695000254621984627876E-2,  -1.601745692712991078208E-2,
    1.496414899205908021882E-3,  -8.462452563778485013756E-5,
    2.728938403476726394024E-6,  -4.239462431819542051337E-8,
};
#endif
#if DEC
static short A4[32] = {
    0137573, 0133037, 0152607, 0113356, 0137427, 0162771, 0145061, 0126345,
    0037224, 0070754, 0110451, 0174104, 0136534, 0164165, 0072170, 0063753,
    0035660, 0000016, 0002560, 0147751, 0134646, 0110311, 0123316, 0047432,
    0033461, 0071250, 0101031, 0075202, 0132066, 0012601, 0077305, 0170177,
};
static short B4[32] = {
    /* 0040200,0000000,0000000,0000000, */
    0137325, 0155602, 0162437, 0030710, 0137611, 0004316, 0071344, 0176361,
    0037306, 0106671, 0011103, 0155053, 0136603, 0033412, 0132530, 0175171,
    0035704, 0021532, 0015516, 0166130, 0134661, 0074162, 0036741, 0073466,
    0033467, 0021316, 0003100, 0171325, 0132066, 0012541, 0162202, 0150160,
};
#endif
#if IBMPC
static short A4[] = {
    0xf2de, 0xfab0, 0x76c3, 0xbfcf, 0x359d, 0x3946, 0xfcbf, 0xbfc2,
    0x3f09, 0x9225, 0x8e3d, 0x3fb2, 0x0cfd, 0xae8f, 0x9d0e, 0xbf8b,
    0x19fd, 0xc0ae, 0x0001, 0x3f56, 0xc9e3, 0x34d9, 0xd219, 0xbf14,
    0x2f50, 0x1043, 0x2e55, 0x3ec6, 0xbe10, 0x2fd8, 0xc2b0, 0xbe66,
};
static short B4[] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0xe639, 0x5ca3, 0xbb70, 0xbfba, 0x9f9e, 0xce5c, 0x2119, 0xbfd1,
    0x7b45, 0x2248, 0xd1b7, 0x3fb8, 0x1f4f, 0x56ab, 0x66e1, 0xbf90,
    0xdd8b, 0x4369, 0x846b, 0x3f58, 0x2ee7, 0x47bc, 0x2f0e, 0xbf16,
    0x1e5b, 0xc0c8, 0xe459, 0x3ec6, 0x5a0e, 0x3c90, 0xc2ac, 0xbe66,
};
#endif
#if MIEEE
static short A4[32] = {
    0xbfcf, 0x76c3, 0xfab0, 0xf2de, 0xbfc2, 0xfcbf, 0x3946, 0x359d,
    0x3fb2, 0x8e3d, 0x9225, 0x3f09, 0xbf8b, 0x9d0e, 0xae8f, 0x0cfd,
    0x3f56, 0x0001, 0xc0ae, 0x19fd, 0xbf14, 0xd219, 0x34d9, 0xc9e3,
    0x3ec6, 0x2e55, 0x1043, 0x2f50, 0xbe66, 0xc2b0, 0x2fd8, 0xbe10,
};
static short B4[32] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0xbfba, 0xbb70, 0x5ca3, 0xe639, 0xbfd1, 0x2119, 0xce5c, 0x9f9e,
    0x3fb8, 0xd1b7, 0x2248, 0x7b45, 0xbf90, 0x66e1, 0x56ab, 0x1f4f,
    0x3f58, 0x846b, 0x4369, 0xdd8b, 0xbf16, 0x2f0e, 0x47bc, 0x2ee7,
    0x3ec6, 0xe459, 0xc0c8, 0x1e5b, 0xbe66, 0xc2ac, 0x3c90, 0x5a0e,
};
#endif

#if 0
/* 20 <= x <= 40
   x exp(-x) Ei(x) - 1  =  1/x A4(1/x) / B4(1/x)
   Theoretical absolute error = 1.78e-17  */
#if UNK
static double A4[8] = {
 2.067245813525780707978E-1,
-5.153749551345223645670E-1,
 1.928289589546695033096E-1,
-3.124468842857260044075E-2,
 2.740283734277352539912E-3,
-1.377775664366875175601E-4,
 3.803788980664744242323E-6,
-4.611038277393688031154E-8,
};
static double B4[8] = {
  /*  1.000000000000000000000E0, */
-8.544436025219516861531E-1,
 2.507436807692907385181E-1,
-3.647688090228423114064E-2,
 3.008576950332041388892E-3,
-1.452926405348421286334E-4,
 3.896007735260115431965E-6,
-4.611037642697098234083E-8,
};
#endif
#if DEC
static short A4[32] = {
0037523,0127633,0150301,0022031,
0140003,0167634,0170572,0170420,
0037505,0072364,0060672,0063220,
0136777,0172334,0057456,0102640,
0036063,0113125,0002476,0047251,
0135020,0074142,0042600,0043630,
0033577,0042230,0155372,0136105,
0132106,0005346,0165333,0114541,
};
static short B4[28] = {
  /* 0040200,0000000,0000000,0000000, */
0140132,0136320,0160433,0131535,
0037600,0060571,0144452,0060214,
0137025,0064310,0024220,0176472,
0036105,0025613,0115762,0166605,
0135030,0054662,0035454,0061763,
0033602,0135163,0116430,0000066,
0132106,0005345,0020602,0137133,
};
#endif
#if IBMPC
static short A4[32] = {
0x2483,0x7a18,0x75f3,0x3fca,
0x5e22,0x9e2f,0x7df3,0xbfe0,
0x4cd2,0x8c37,0xae9e,0x3fc8,
0xd0b4,0x8be5,0xfe9b,0xbf9f,
0xc9d5,0xa0a7,0x72ca,0x3f66,
0x08f3,0x48b0,0x0f0c,0xbf22,
0x5789,0x1b5f,0xe893,0x3ecf,
0x732c,0xdd5b,0xc15c,0xbe68,
};
static short B4[28] = {
  /* 0x0000,0x0000,0x0000,0x3ff0, */
0x766c,0x1c23,0x579a,0xbfeb,
0x4c11,0x3925,0x0c2f,0x3fd0,
0x1fa7,0x0512,0xad19,0xbfa2,
0x5db1,0x737e,0xa571,0x3f68,
0x8c7e,0x4765,0x0b36,0xbf23,
0x0007,0x73a3,0x574e,0x3ed0,
0x57cb,0xa430,0xc15c,0xbe68,
};
#endif
#if MIEEE
static short A4[32] = {
0x3fca,0x75f3,0x7a18,0x2483,
0xbfe0,0x7df3,0x9e2f,0x5e22,
0x3fc8,0xae9e,0x8c37,0x4cd2,
0xbf9f,0xfe9b,0x8be5,0xd0b4,
0x3f66,0x72ca,0xa0a7,0xc9d5,
0xbf22,0x0f0c,0x48b0,0x08f3,
0x3ecf,0xe893,0x1b5f,0x5789,
0xbe68,0xc15c,0xdd5b,0x732c,
};
static short B4[28] = {
  /* 0x3ff0,0x0000,0x0000,0x0000, */
0xbfeb,0x579a,0x1c23,0x766c,
0x3fd0,0x0c2f,0x3925,0x4c11,
0xbfa2,0xad19,0x0512,0x1fa7,
0x3f68,0xa571,0x737e,0x5db1,
0xbf23,0x0b36,0x4765,0x8c7e,
0x3ed0,0x574e,0x73a3,0x0007,
0xbe68,0xc15c,0xa430,0x57cb,
};
#endif
#endif /* 0 */

/* 4 <= x <= 8
   x exp(-x) Ei(x) - 1  =  1/x A5(1/x) / B5(1/x)
   Theoretical absolute error = 2.20e-17  */
#if UNK
static double A5[8] = {
    -1.373215375871208729803E0, -7.084559133740838761406E-1,
    1.580806855547941010501E0,  -2.601500427425622944234E-1,
    2.994674694113713763365E-2, -1.038086040188744005513E-3,
    4.371064420753005429514E-5, 2.141783679522602903795E-6,
};
static double B5[8] = {
    /* 1.000000000000000000000E0, */
    8.585231423622028380768E-1,  4.483285822873995129957E-1,
    7.687932158124475434091E-2,  2.449868241021887685904E-2,
    8.832165941927796567926E-4,  4.590952299511353531215E-4,
    -4.729848351866523044863E-6, 2.665195537390710170105E-6,
};
#endif
#if DEC
static short A5[32] = {
    0140257, 0142605, 0076335, 0113632, 0140065, 0056535, 0161231, 0074311,
    0040312, 0053741, 0004357, 0076405, 0137605, 0031142, 0165503, 0136705,
    0036765, 0051341, 0053573, 0007602, 0135610, 0010143, 0027643, 0110522,
    0034467, 0052762, 0062024, 0120161, 0033417, 0135620, 0036500, 0062647,
};
static short B[32] = {
    /* 0040200,0000000,0000000,0000000, */
    0040133, 0144054, 0031516, 0004100, 0037745, 0105522, 0166622, 0123146,
    0037235, 0071347, 0157560, 0157464, 0036710, 0130565, 0173747, 0041670,
    0035547, 0103651, 0106243, 0101240, 0035360, 0131267, 0176263, 0140257,
    0133636, 0132426, 0102537, 0102531, 0033462, 0155665, 0167503, 0176350,
};
#endif
#if IBMPC
static short A5[32] = {
    0xb2f3, 0xaf9b, 0xf8b0, 0xbff5, 0x2f19, 0xbc53, 0xabab, 0xbfe6,
    0xefa1, 0x211d, 0x4afc, 0x3ff9, 0x77b9, 0x5d68, 0xa64c, 0xbfd0,
    0x61f0, 0x2aef, 0xaa5c, 0x3f9e, 0x722a, 0x65f4, 0x020c, 0xbf51,
    0x940e, 0x4c82, 0xeabe, 0x3f06, 0x0cb5, 0x07a8, 0xf772, 0x3ec1,
};
static short B5[32] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0xc108, 0x8669, 0x7905, 0x3feb, 0x54cd, 0x5db2, 0xb16a, 0x3fdc,
    0x1be7, 0xfbee, 0xae5c, 0x3fb3, 0xe877, 0xbefc, 0x162e, 0x3f99,
    0x7054, 0x3194, 0xf0f5, 0x3f4c, 0x7816, 0xff96, 0x1656, 0x3f3e,
    0xf0ab, 0xd0ab, 0xd6a2, 0xbed3, 0x7f9d, 0xbde8, 0x5b76, 0x3ec6,
};
#endif
#if MIEEE
static short A5[32] = {
    0xbff5, 0xf8b0, 0xaf9b, 0xb2f3, 0xbfe6, 0xabab, 0xbc53, 0x2f19,
    0x3ff9, 0x4afc, 0x211d, 0xefa1, 0xbfd0, 0xa64c, 0x5d68, 0x77b9,
    0x3f9e, 0xaa5c, 0x2aef, 0x61f0, 0xbf51, 0x020c, 0x65f4, 0x722a,
    0x3f06, 0xeabe, 0x4c82, 0x940e, 0x3ec1, 0xf772, 0x07a8, 0x0cb5,
};
static short B5[32] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0x3feb, 0x7905, 0x8669, 0xc108, 0x3fdc, 0xb16a, 0x5db2, 0x54cd,
    0x3fb3, 0xae5c, 0xfbee, 0x1be7, 0x3f99, 0x162e, 0xbefc, 0xe877,
    0x3f4c, 0xf0f5, 0x3194, 0x7054, 0x3f3e, 0x1656, 0xff96, 0x7816,
    0xbed3, 0xd6a2, 0xd0ab, 0xf0ab, 0x3ec6, 0x5b76, 0xbde8, 0x7f9d,
};
#endif
/* 2 <= x <= 4
   x exp(-x) Ei(x) - 1  =  1/x A6(1/x) / B6(1/x)
   Theoretical absolute error = 4.89e-17  */
#if UNK
static double A6[8] = {
    1.981808503259689673238E-2,  -1.271645625984917501326E0,
    -2.088160335681228318920E0,  2.755544509187936721172E0,
    -4.409507048701600257171E-1, 4.665623805935891391017E-2,
    -1.545042679673485262580E-3, 7.059980605299617478514E-5,
};
static double B6[7] = {
    /* 1.000000000000000000000E0, */
    1.476498670914921440652E0,  5.629177174822436244827E-1,
    1.699017897879307263248E-1, 2.291647179034212017463E-2,
    4.450150439728752875043E-3, 1.727439612206521482874E-4,
    3.953167195549672482304E-5,
};
#endif
#if DEC
static short A6[32] = {
    0036642, 0054611, 0061263, 0000140, 0140242, 0142510, 0125732, 0072035,
    0140405, 0122153, 0037643, 0104527, 0040460, 0055327, 0055550, 0116240,
    0137741, 0142112, 0070441, 0103510, 0037077, 0015234, 0104750, 0146765,
    0135712, 0101407, 0107554, 0020253, 0034624, 0007373, 0072621, 0063735,
};
static short B6[28] = {
    /* 0040200,0000000,0000000,0000000, */
    0040274, 0176750, 0110025, 0061006, 0040020, 0015540, 0021354,
    0155050, 0037455, 0175274, 0015257, 0021112, 0036673, 0135523,
    0016042, 0117203, 0036221, 0151221, 0046352, 0144174, 0035065,
    0021232, 0117727, 0152432, 0034445, 0147317, 0037300, 0067123,
};
#endif
#if IBMPC
static short A6[32] = {
    0x600c, 0x2c56, 0x4b31, 0x3f94, 0x4e84, 0x157b, 0x58a9, 0xbff4,
    0x712b, 0x67f4, 0xb48d, 0xc000, 0x1394, 0xeb6d, 0x0b5a, 0x4006,
    0x30e9, 0x4e24, 0x3889, 0xbfdc, 0x19bf, 0x913d, 0xe353, 0x3fa7,
    0x8415, 0xf1ed, 0x5060, 0xbf59, 0x2cfc, 0x6eb2, 0x81df, 0x3f12,
};
static short B6[28] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0xac41, 0x1202, 0x9fbd, 0x3ff7, 0x9b45, 0x045d, 0x036c,
    0x3fe2, 0xe449, 0x8355, 0xbf57, 0x3fc5, 0x53d0, 0x6384,
    0x776a, 0x3f97, 0x590f, 0x299d, 0x3a52, 0x3f72, 0xfaa3,
    0x53fa, 0xa453, 0x3f26, 0x0dca, 0xe7d8, 0xb9d9, 0x3f04,
};
#endif
#if MIEEE
static short A6[32] = {
    0x3f94, 0x4b31, 0x2c56, 0x600c, 0xbff4, 0x58a9, 0x157b, 0x4e84,
    0xc000, 0xb48d, 0x67f4, 0x712b, 0x4006, 0x0b5a, 0xeb6d, 0x1394,
    0xbfdc, 0x3889, 0x4e24, 0x30e9, 0x3fa7, 0xe353, 0x913d, 0x19bf,
    0xbf59, 0x5060, 0xf1ed, 0x8415, 0x3f12, 0x81df, 0x6eb2, 0x2cfc,
};
static short B6[28] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0x3ff7, 0x9fbd, 0x1202, 0xac41, 0x3fe2, 0x036c, 0x045d,
    0x9b45, 0x3fc5, 0xbf57, 0x8355, 0xe449, 0x3f97, 0x776a,
    0x6384, 0x53d0, 0x3f72, 0x3a52, 0x299d, 0x590f, 0x3f26,
    0xa453, 0x53fa, 0xfaa3, 0x3f04, 0xb9d9, 0xe7d8, 0x0dca,
};
#endif
/* 32 <= x <= 64
   x exp(-x) Ei(x) - 1  =  1/x A7(1/x) / B7(1/x)
   Theoretical absolute error = 7.71e-18  */
#if UNK
static double A7[6] = {
    1.212561118105456670844E-1, -5.823133179043894485122E-1,
    2.348887314557016779211E-1, -3.040034318113248237280E-2,
    1.510082146865190661777E-3, -2.523137095499571377122E-5,
};
static double B7[5] = {
    /* 1.000000000000000000000E0, */
    -1.002252150365854016662E0,  2.928709694872224144953E-1,
    -3.337004338674007801307E-2, 1.560544881127388842819E-3,
    -2.523137093603234562648E-5,
};
#endif
#if DEC
static short A7[24] = {
    0037370, 0052437, 0152524, 0150125, 0140025, 0011174, 0050154, 0131330,
    0037560, 0103253, 0167464, 0062245, 0136771, 0005043, 0174001, 0023345,
    0035705, 0166762, 0157300, 0016451, 0134323, 0123764, 0157767, 0134477,
};
static short B7[20] = {
    /* 0040200,0000000,0000000,0000000, */
    0140200, 0044714, 0064025, 0060324, 0037625, 0171457, 0003712,
    0073131, 0137010, 0127406, 0150061, 0141746, 0035714, 0105462,
    0072356, 0103712, 0134323, 0123764, 0156514, 0077414,
};
#endif
#if IBMPC
static short A7[24] = {
    0x9a0b, 0xfaaa, 0x0aa3, 0x3fbf, 0x965b, 0x8a0d, 0xa24f, 0xbfe2,
    0x8c95, 0x7de6, 0x10d5, 0x3fce, 0x24dd, 0x7f00, 0x2144, 0xbf9f,
    0x03a5, 0x5bd8, 0xbdbe, 0x3f58, 0xf728, 0x9bfe, 0x74fe, 0xbefa,
};
static short B7[20] = {
    /* 0x0000,0x0000,0x0000,0x3ff0, */
    0xac1a, 0x8d02, 0x0939, 0xbff0, 0x4ecb, 0xe0f9, 0xbe65,
    0x3fd2, 0x387d, 0xda06, 0x15e0, 0xbfa1, 0xd0f9, 0x4e9d,
    0x9166, 0x3f59, 0x8fe2, 0x9ba9, 0x74fe, 0xbefa,
};
#endif
#if MIEEE
static short A7[24] = {
    0x3fbf, 0x0aa3, 0xfaaa, 0x9a0b, 0xbfe2, 0xa24f, 0x8a0d, 0x965b,
    0x3fce, 0x10d5, 0x7de6, 0x8c95, 0xbf9f, 0x2144, 0x7f00, 0x24dd,
    0x3f58, 0xbdbe, 0x5bd8, 0x03a5, 0xbefa, 0x74fe, 0x9bfe, 0xf728,
};
static short B7[20] = {
    /* 0x3ff0,0x0000,0x0000,0x0000, */
    0xbff0, 0x0939, 0x8d02, 0xac1a, 0x3fd2, 0xbe65, 0xe0f9,
    0x4ecb, 0xbfa1, 0x15e0, 0xda06, 0x387d, 0x3f59, 0x9166,
    0x4e9d, 0xd0f9, 0xbefa, 0x74fe, 0x9ba9, 0x8fe2,
};
#endif

double ei(x) double x;
{
  double f, w;

  if (x <= 0.0) {
    mtherr("ei", DOMAIN);
    return 0.0;
  } else if (x < 2.0) {
    /* Power series.
                              inf    n
                               -    x
       Ei(x) = EUL + ln x  +   >   ----
                               -   n n!
                              n=1
    */
    f = polevl(x, A, 5) / p1evl(x, B, 6);
    /*      f = polevl(x,A,6) / p1evl(x,B,7); */
    /*      f = polevl(x,A,8) / p1evl(x,B,9); */
    return (EUL + log(x) + x * f);
  } else if (x < 4.0) {
    /* Asymptotic expansion.
                              1       2       6
      x exp(-x) Ei(x) =  1 + ---  +  ---  +  ---- + ...
                              x        2       3
                                      x       x
    */
    w = 1.0 / x;
    f = polevl(w, A6, 7) / p1evl(w, B6, 7);
    return (exp(x) * w * (1.0 + w * f));
  } else if (x < 8.0) {
    w = 1.0 / x;
    f = polevl(w, A5, 7) / p1evl(w, B5, 8);
    return (exp(x) * w * (1.0 + w * f));
  } else if (x < 16.0) {
    w = 1.0 / x;
    f = polevl(w, A2, 9) / p1evl(w, B2, 9);
    return (exp(x) * w * (1.0 + w * f));
  } else if (x < 32.0) {
    w = 1.0 / x;
    f = polevl(w, A4, 7) / p1evl(w, B4, 8);
    return (exp(x) * w * (1.0 + w * f));
  } else if (x < 64.0) {
    w = 1.0 / x;
    f = polevl(w, A7, 5) / p1evl(w, B7, 5);
    return (exp(x) * w * (1.0 + w * f));
  } else {
    w = 1.0 / x;
    f = polevl(w, A3, 8) / p1evl(w, B3, 9);
    return (exp(x) * w * (1.0 + w * f));
  }
}