diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-05-10 00:40:28 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-05-10 00:40:28 +0000 |
commit | 1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch) | |
tree | 579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/ei.c | |
parent | 22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff) |
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/ei.c')
-rw-r--r-- | libm/double/ei.c | 1062 |
1 files changed, 1062 insertions, 0 deletions
diff --git a/libm/double/ei.c b/libm/double/ei.c new file mode 100644 index 000000000..4994fa99c --- /dev/null +++ b/libm/double/ei.c @@ -0,0 +1,1062 @@ +/* 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 <math.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)); + } +} |