123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835 |
- /* 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));
- }
- }
|