123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375 |
- /* shichi.c
- *
- * Hyperbolic sine and cosine integrals
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, Chi, Shi, shichi();
- *
- * shichi( x, &Chi, &Shi );
- *
- *
- * DESCRIPTION:
- *
- * Approximates the integrals
- *
- * x
- * -
- * | | cosh t - 1
- * Chi(x) = eul + ln x + | ----------- dt,
- * | | t
- * -
- * 0
- *
- * x
- * -
- * | | sinh t
- * Shi(x) = | ------ dt
- * | | t
- * -
- * 0
- *
- * where eul = 0.57721566490153286061 is Euler's constant.
- * The integrals are evaluated by power series for x < 8
- * and by Chebyshev expansions for x between 8 and 88.
- * For large x, both functions approach exp(x)/2x.
- * Arguments greater than 88 in magnitude return MAXNUM.
- *
- *
- * ACCURACY:
- *
- * Test interval 0 to 88.
- * Relative error:
- * arithmetic function # trials peak rms
- * DEC Shi 3000 9.1e-17
- * IEEE Shi 30000 6.9e-16 1.6e-16
- * Absolute error, except relative when |Chi| > 1:
- * DEC Chi 2500 9.3e-17
- * IEEE Chi 30000 8.4e-16 1.4e-16
- */
- /*
- Cephes Math Library Release 2.8: June, 2000
- Copyright 1984, 1987, 2000 by Stephen L. Moshier
- */
- #include "mconf.h"
- #ifdef UNK
- /* x exp(-x) shi(x), inverted interval 8 to 18 */
- static double S1[] = {1.83889230173399459482E-17, -9.55485532279655569575E-17,
- 2.04326105980879882648E-16, 1.09896949074905343022E-15,
- -1.31313534344092599234E-14, 5.93976226264314278932E-14,
- -3.47197010497749154755E-14, -1.40059764613117131000E-12,
- 9.49044626224223543299E-12, -1.61596181145435454033E-11,
- -1.77899784436430310321E-10, 1.35455469767246947469E-9,
- -1.03257121792819495123E-9, -3.56699611114982536845E-8,
- 1.44818877384267342057E-7, 7.82018215184051295296E-7,
- -5.39919118403805073710E-6, -3.12458202168959833422E-5,
- 8.90136741950727517826E-5, 2.02558474743846862168E-3,
- 2.96064440855633256972E-2, 1.11847751047257036625E0};
- /* x exp(-x) shi(x), inverted interval 18 to 88 */
- static double S2[] = {-1.05311574154850938805E-17, 2.62446095596355225821E-17,
- 8.82090135625368160657E-17, -3.38459811878103047136E-16,
- -8.30608026366935789136E-16, 3.93397875437050071776E-15,
- 1.01765565969729044505E-14, -4.21128170307640802703E-14,
- -1.60818204519802480035E-13, 3.34714954175994481761E-13,
- 2.72600352129153073807E-12, 1.66894954752839083608E-12,
- -3.49278141024730899554E-11, -1.58580661666482709598E-10,
- -1.79289437183355633342E-10, 1.76281629144264523277E-9,
- 1.69050228879421288846E-8, 1.25391771228487041649E-7,
- 1.16229947068677338732E-6, 1.61038260117376323993E-5,
- 3.49810375601053973070E-4, 1.28478065259647610779E-2,
- 1.03665722588798326712E0};
- #endif
- #ifdef DEC
- static unsigned short S1[] = {
- 0022251, 0115635, 0165120, 0006574, 0122734, 0050751, 0020305, 0101356,
- 0023153, 0111154, 0011103, 0177462, 0023636, 0060321, 0060253, 0124246,
- 0124554, 0106655, 0152525, 0166400, 0025205, 0140145, 0171006, 0106556,
- 0125034, 0056427, 0004205, 0176022, 0126305, 0016731, 0025011, 0134453,
- 0027046, 0172453, 0112604, 0116235, 0127216, 0022071, 0116600, 0137667,
- 0130103, 0115126, 0071104, 0052535, 0030672, 0025450, 0010071, 0141414,
- 0130615, 0165136, 0132137, 0177737, 0132031, 0031611, 0074436, 0175407,
- 0032433, 0077602, 0104345, 0060076, 0033121, 0165741, 0167177, 0172433,
- 0133665, 0025262, 0174621, 0022612, 0134403, 0006761, 0124566, 0145405,
- 0034672, 0126332, 0034737, 0116744, 0036004, 0137654, 0037332, 0131766,
- 0036762, 0104466, 0121445, 0124326, 0040217, 0025105, 0062145, 0042640};
- static unsigned short S2[] = {
- 0122102, 0041774, 0016051, 0055137, 0022362, 0010125, 0007651, 0015773,
- 0022713, 0062551, 0040227, 0071645, 0123303, 0015732, 0025731, 0146570,
- 0123557, 0064016, 0002067, 0067711, 0024215, 0136214, 0132374, 0124234,
- 0024467, 0051425, 0071066, 0064210, 0125075, 0124305, 0135123, 0024170,
- 0125465, 0010261, 0005560, 0034232, 0025674, 0066602, 0030724, 0174557,
- 0026477, 0151520, 0051510, 0067250, 0026352, 0161076, 0113154, 0116271,
- 0127431, 0116470, 0177465, 0127274, 0130056, 0056174, 0170315, 0013321,
- 0130105, 0020575, 0075327, 0036710, 0030762, 0043625, 0113046, 0125035,
- 0031621, 0033211, 0154354, 0022077, 0032406, 0121555, 0074270, 0041141,
- 0033234, 0000116, 0041611, 0173743, 0034207, 0013263, 0174715, 0115563,
- 0035267, 0063300, 0175753, 0117266, 0036522, 0077633, 0033255, 0136200,
- 0040204, 0130457, 0014454, 0166254};
- #endif
- #ifdef IBMPC
- static unsigned short S1[] = {
- 0x01b0, 0xbd4a, 0x3373, 0x3c75, 0xb05e, 0x2418, 0x8a3d, 0xbc9b, 0x7fe6,
- 0x8248, 0x724d, 0x3cad, 0x7515, 0x2c15, 0xcc1a, 0x3cd3, 0xbda0, 0xbaaa,
- 0x91b5, 0xbd0d, 0xd1ae, 0xbe40, 0xb80c, 0x3d30, 0xbf82, 0xe110, 0x8ba2,
- 0xbd23, 0x3725, 0x2541, 0xa3bb, 0xbd78, 0x9394, 0x72b0, 0xdea5, 0x3da4,
- 0x17f7, 0x33b0, 0xc487, 0xbdb1, 0x8aac, 0xce48, 0x734a, 0xbde8, 0x3862,
- 0x0207, 0x4565, 0x3e17, 0xfffc, 0xd68b, 0xbd4b, 0xbe11, 0xdf61, 0x2f23,
- 0x2671, 0xbe63, 0xac08, 0x511c, 0x6ff0, 0x3e83, 0xfea3, 0x3dcf, 0x3d7c,
- 0x3eaa, 0x24b1, 0x5f32, 0xa556, 0xbed6, 0xd961, 0x352e, 0x61be, 0xbf00,
- 0xf3bd, 0x473b, 0x559b, 0x3f17, 0x567f, 0x87db, 0x97f5, 0x3f60, 0xb51b,
- 0xd464, 0x5126, 0x3f9e, 0xa8b4, 0xac8c, 0xe548, 0x3ff1};
- static unsigned short S2[] = {
- 0x2b4c, 0x8385, 0x487f, 0xbc68, 0x237f, 0xa1f5, 0x420a, 0x3c7e, 0xee75,
- 0x2812, 0x6cad, 0x3c99, 0x39af, 0x457b, 0x637b, 0xbcb8, 0xedf9, 0xc086,
- 0xed01, 0xbccd, 0x9513, 0x969f, 0xb791, 0x3cf1, 0xcd11, 0xae46, 0xea62,
- 0x3d06, 0x650f, 0xb74a, 0xb518, 0xbd27, 0x0713, 0x216e, 0xa216, 0xbd46,
- 0x9f2e, 0x463a, 0x8db0, 0x3d57, 0x0dd5, 0x0a69, 0xfa6a, 0x3d87, 0x9397,
- 0xd2cd, 0x5c47, 0x3d7d, 0xb5d8, 0x1fe6, 0x33a7, 0xbdc3, 0xa2da, 0x9e19,
- 0xcb8f, 0xbde5, 0xe7b9, 0xaf5a, 0xa42f, 0xbde8, 0xd544, 0xb2c4, 0x48f2,
- 0x3e1e, 0x8488, 0x3b1d, 0x26d1, 0x3e52, 0x084c, 0xaf17, 0xd46d, 0x3e80,
- 0x3efc, 0xc871, 0x8009, 0x3eb3, 0xb36e, 0x7f39, 0xe2d6, 0x3ef0, 0x73d7,
- 0x1f7d, 0xecd8, 0x3f36, 0xb790, 0x66d5, 0x4ff3, 0x3f8a, 0x9d96, 0xe325,
- 0x9625, 0x3ff0};
- #endif
- #ifdef MIEEE
- static unsigned short S1[] = {
- 0x3c75, 0x3373, 0xbd4a, 0x01b0, 0xbc9b, 0x8a3d, 0x2418, 0xb05e, 0x3cad,
- 0x724d, 0x8248, 0x7fe6, 0x3cd3, 0xcc1a, 0x2c15, 0x7515, 0xbd0d, 0x91b5,
- 0xbaaa, 0xbda0, 0x3d30, 0xb80c, 0xbe40, 0xd1ae, 0xbd23, 0x8ba2, 0xe110,
- 0xbf82, 0xbd78, 0xa3bb, 0x2541, 0x3725, 0x3da4, 0xdea5, 0x72b0, 0x9394,
- 0xbdb1, 0xc487, 0x33b0, 0x17f7, 0xbde8, 0x734a, 0xce48, 0x8aac, 0x3e17,
- 0x4565, 0x0207, 0x3862, 0xbe11, 0xbd4b, 0xd68b, 0xfffc, 0xbe63, 0x2671,
- 0x2f23, 0xdf61, 0x3e83, 0x6ff0, 0x511c, 0xac08, 0x3eaa, 0x3d7c, 0x3dcf,
- 0xfea3, 0xbed6, 0xa556, 0x5f32, 0x24b1, 0xbf00, 0x61be, 0x352e, 0xd961,
- 0x3f17, 0x559b, 0x473b, 0xf3bd, 0x3f60, 0x97f5, 0x87db, 0x567f, 0x3f9e,
- 0x5126, 0xd464, 0xb51b, 0x3ff1, 0xe548, 0xac8c, 0xa8b4};
- static unsigned short S2[] = {
- 0xbc68, 0x487f, 0x8385, 0x2b4c, 0x3c7e, 0x420a, 0xa1f5, 0x237f, 0x3c99,
- 0x6cad, 0x2812, 0xee75, 0xbcb8, 0x637b, 0x457b, 0x39af, 0xbccd, 0xed01,
- 0xc086, 0xedf9, 0x3cf1, 0xb791, 0x969f, 0x9513, 0x3d06, 0xea62, 0xae46,
- 0xcd11, 0xbd27, 0xb518, 0xb74a, 0x650f, 0xbd46, 0xa216, 0x216e, 0x0713,
- 0x3d57, 0x8db0, 0x463a, 0x9f2e, 0x3d87, 0xfa6a, 0x0a69, 0x0dd5, 0x3d7d,
- 0x5c47, 0xd2cd, 0x9397, 0xbdc3, 0x33a7, 0x1fe6, 0xb5d8, 0xbde5, 0xcb8f,
- 0x9e19, 0xa2da, 0xbde8, 0xa42f, 0xaf5a, 0xe7b9, 0x3e1e, 0x48f2, 0xb2c4,
- 0xd544, 0x3e52, 0x26d1, 0x3b1d, 0x8488, 0x3e80, 0xd46d, 0xaf17, 0x084c,
- 0x3eb3, 0x8009, 0xc871, 0x3efc, 0x3ef0, 0xe2d6, 0x7f39, 0xb36e, 0x3f36,
- 0xecd8, 0x1f7d, 0x73d7, 0x3f8a, 0x4ff3, 0x66d5, 0xb790, 0x3ff0, 0x9625,
- 0xe325, 0x9d96};
- #endif
- #ifdef UNK
- /* x exp(-x) chin(x), inverted interval 8 to 18 */
- static double C1[] = {-8.12435385225864036372E-18, 2.17586413290339214377E-17,
- 5.22624394924072204667E-17, -9.48812110591690559363E-16,
- 5.35546311647465209166E-15, -1.21009970113732918701E-14,
- -6.00865178553447437951E-14, 7.16339649156028587775E-13,
- -2.93496072607599856104E-12, -1.40359438136491256904E-12,
- 8.76302288609054966081E-11, -4.40092476213282340617E-10,
- -1.87992075640569295479E-10, 1.31458150989474594064E-8,
- -4.75513930924765465590E-8, -2.21775018801848880741E-7,
- 1.94635531373272490962E-6, 4.33505889257316408893E-6,
- -6.13387001076494349496E-5, -3.13085477492997465138E-4,
- 4.97164789823116062801E-4, 2.64347496031374526641E-2,
- 1.11446150876699213025E0};
- /* x exp(-x) chin(x), inverted interval 18 to 88 */
- static double C2[] = {8.06913408255155572081E-18, -2.08074168180148170312E-17,
- -5.98111329658272336816E-17, 2.68533951085945765591E-16,
- 4.52313941698904694774E-16, -3.10734917335299464535E-15,
- -4.42823207332531972288E-15, 3.49639695410806959872E-14,
- 6.63406731718911586609E-14, -3.71902448093119218395E-13,
- -1.27135418132338309016E-12, 2.74851141935315395333E-12,
- 2.33781843985453438400E-11, 2.71436006377612442764E-11,
- -2.56600180000355990529E-10, -1.61021375163803438552E-9,
- -4.72543064876271773512E-9, -3.00095178028681682282E-9,
- 7.79387474390914922337E-8, 1.06942765566401507066E-6,
- 1.59503164802313196374E-5, 3.49592575153777996871E-4,
- 1.28475387530065247392E-2, 1.03665693917934275131E0};
- #endif
- #ifdef DEC
- static unsigned short C1[] = {
- 0122025, 0157055, 0021702, 0021427, 0022310, 0130043, 0123265, 0022340,
- 0022561, 0002231, 0017746, 0013043, 0123610, 0136375, 0002352, 0024467,
- 0024300, 0171555, 0141300, 0000446, 0124531, 0176777, 0126210, 0035616,
- 0125207, 0046604, 0167760, 0077132, 0026111, 0120666, 0026606, 0064143,
- 0126516, 0103615, 0054127, 0005436, 0126305, 0104721, 0025415, 0004134,
- 0027700, 0131556, 0164725, 0157553, 0130361, 0170602, 0077274, 0055406,
- 0130116, 0131420, 0125472, 0017231, 0031541, 0153747, 0177312, 0056304,
- 0132114, 0035517, 0041545, 0043151, 0132556, 0020415, 0110044, 0172442,
- 0033402, 0117041, 0031152, 0010364, 0033621, 0072737, 0050647, 0013720,
- 0134600, 0121366, 0140010, 0063265, 0135244, 0022637, 0013756, 0044742,
- 0035402, 0052052, 0006523, 0043564, 0036730, 0106660, 0020277, 0162146,
- 0040216, 0123254, 0135147, 0005724};
- static unsigned short C2[] = {
- 0022024, 0154550, 0104311, 0144257, 0122277, 0165037, 0133443, 0155601,
- 0122611, 0165102, 0157053, 0055252, 0023232, 0146235, 0153511, 0113222,
- 0023402, 0057340, 0145304, 0010471, 0124137, 0164171, 0113071, 0100002,
- 0124237, 0105473, 0056130, 0022022, 0025035, 0073266, 0056746, 0164433,
- 0025225, 0061313, 0055600, 0165407, 0125721, 0056312, 0107613, 0051215,
- 0126262, 0166534, 0115336, 0066653, 0026501, 0064307, 0127442, 0065573,
- 0027315, 0121375, 0142020, 0045356, 0027356, 0140764, 0070641, 0046570,
- 0130215, 0010503, 0146335, 0177737, 0130735, 0047134, 0015215, 0163665,
- 0131242, 0056523, 0155276, 0050053, 0131116, 0034515, 0050707, 0163512,
- 0032247, 0057507, 0107545, 0032007, 0033217, 0104501, 0021706, 0025047,
- 0034205, 0146413, 0033746, 0076562, 0035267, 0044605, 0065355, 0002772,
- 0036522, 0077173, 0130716, 0170304, 0040204, 0130454, 0130571, 0027270};
- #endif
- #ifdef IBMPC
- static unsigned short C1[] = {
- 0x4463, 0xa478, 0xbbc5, 0xbc62, 0xa49c, 0x74d6, 0x1604, 0x3c79, 0xc2c4,
- 0x23fc, 0x2093, 0x3c8e, 0x4527, 0xa09d, 0x179f, 0xbcd1, 0x0025, 0xb858,
- 0x1e6d, 0x3cf8, 0x0772, 0xf591, 0x3fbf, 0xbd0b, 0x0fcb, 0x9dfe, 0xe9b0,
- 0xbd30, 0xcd0c, 0xc5b0, 0x3436, 0x3d69, 0xe164, 0xab0a, 0xd0f1, 0xbd89,
- 0xa10c, 0x2561, 0xb13a, 0xbd78, 0xbbed, 0xdd3a, 0x166d, 0x3dd8, 0x8b61,
- 0x4fd7, 0x3e30, 0xbdfe, 0x43d3, 0x1567, 0xd662, 0xbde9, 0x4b98, 0xffd9,
- 0x3afc, 0x3e4c, 0xa8cd, 0xe86c, 0x8769, 0xbe69, 0x9ea4, 0xb204, 0xc421,
- 0xbe8d, 0x421f, 0x264d, 0x53c4, 0x3ec0, 0xe2fa, 0xea34, 0x2ebb, 0x3ed2,
- 0x0cd7, 0xd801, 0x145e, 0xbf10, 0xc93c, 0xe2fd, 0x84b3, 0xbf34, 0x68ef,
- 0x41aa, 0x4a85, 0x3f40, 0xfc8d, 0x0417, 0x11b6, 0x3f9b, 0xe17b, 0x974c,
- 0xd4d5, 0x3ff1};
- static unsigned short C2[] = {
- 0x3916, 0x1119, 0x9b2d, 0x3c62, 0x7b70, 0xf6e4, 0xfd43, 0xbc77, 0x6b55,
- 0x5bc5, 0x3d48, 0xbc91, 0x32d2, 0xbae9, 0x5993, 0x3cb3, 0x8227, 0x1958,
- 0x4bdc, 0x3cc0, 0x3000, 0x32c7, 0xfd0f, 0xbceb, 0x0482, 0x6b8b, 0xf167,
- 0xbcf3, 0xdd23, 0xcbbc, 0xaed6, 0x3d23, 0x1d61, 0x6b70, 0xac59, 0x3d32,
- 0x6a52, 0x51f1, 0x2b99, 0xbd5a, 0xcdb5, 0x935b, 0x5dab, 0xbd76, 0x4d6f,
- 0xf5e4, 0x2d18, 0x3d88, 0x095e, 0xb882, 0xb45f, 0x3db9, 0x29af, 0x8e34,
- 0xd83e, 0x3dbd, 0xbffc, 0x799b, 0xa228, 0xbdf1, 0xbcf7, 0x8351, 0xa9cb,
- 0xbe1b, 0xca05, 0x7b57, 0x4baa, 0xbe34, 0xfce9, 0xaa38, 0xc729, 0xbe29,
- 0xa681, 0xf1ec, 0xebe8, 0x3e74, 0xc545, 0x2478, 0xf128, 0x3eb1, 0xcfae,
- 0x66fc, 0xb9a1, 0x3ef0, 0xa0bf, 0xad5d, 0xe930, 0x3f36, 0xde19, 0x7639,
- 0x4fcf, 0x3f8a, 0x25d7, 0x962f, 0x9625, 0x3ff0};
- #endif
- #ifdef MIEEE
- static unsigned short C1[] = {
- 0xbc62, 0xbbc5, 0xa478, 0x4463, 0x3c79, 0x1604, 0x74d6, 0xa49c, 0x3c8e,
- 0x2093, 0x23fc, 0xc2c4, 0xbcd1, 0x179f, 0xa09d, 0x4527, 0x3cf8, 0x1e6d,
- 0xb858, 0x0025, 0xbd0b, 0x3fbf, 0xf591, 0x0772, 0xbd30, 0xe9b0, 0x9dfe,
- 0x0fcb, 0x3d69, 0x3436, 0xc5b0, 0xcd0c, 0xbd89, 0xd0f1, 0xab0a, 0xe164,
- 0xbd78, 0xb13a, 0x2561, 0xa10c, 0x3dd8, 0x166d, 0xdd3a, 0xbbed, 0xbdfe,
- 0x3e30, 0x4fd7, 0x8b61, 0xbde9, 0xd662, 0x1567, 0x43d3, 0x3e4c, 0x3afc,
- 0xffd9, 0x4b98, 0xbe69, 0x8769, 0xe86c, 0xa8cd, 0xbe8d, 0xc421, 0xb204,
- 0x9ea4, 0x3ec0, 0x53c4, 0x264d, 0x421f, 0x3ed2, 0x2ebb, 0xea34, 0xe2fa,
- 0xbf10, 0x145e, 0xd801, 0x0cd7, 0xbf34, 0x84b3, 0xe2fd, 0xc93c, 0x3f40,
- 0x4a85, 0x41aa, 0x68ef, 0x3f9b, 0x11b6, 0x0417, 0xfc8d, 0x3ff1, 0xd4d5,
- 0x974c, 0xe17b};
- static unsigned short C2[] = {
- 0x3c62, 0x9b2d, 0x1119, 0x3916, 0xbc77, 0xfd43, 0xf6e4, 0x7b70, 0xbc91,
- 0x3d48, 0x5bc5, 0x6b55, 0x3cb3, 0x5993, 0xbae9, 0x32d2, 0x3cc0, 0x4bdc,
- 0x1958, 0x8227, 0xbceb, 0xfd0f, 0x32c7, 0x3000, 0xbcf3, 0xf167, 0x6b8b,
- 0x0482, 0x3d23, 0xaed6, 0xcbbc, 0xdd23, 0x3d32, 0xac59, 0x6b70, 0x1d61,
- 0xbd5a, 0x2b99, 0x51f1, 0x6a52, 0xbd76, 0x5dab, 0x935b, 0xcdb5, 0x3d88,
- 0x2d18, 0xf5e4, 0x4d6f, 0x3db9, 0xb45f, 0xb882, 0x095e, 0x3dbd, 0xd83e,
- 0x8e34, 0x29af, 0xbdf1, 0xa228, 0x799b, 0xbffc, 0xbe1b, 0xa9cb, 0x8351,
- 0xbcf7, 0xbe34, 0x4baa, 0x7b57, 0xca05, 0xbe29, 0xc729, 0xaa38, 0xfce9,
- 0x3e74, 0xebe8, 0xf1ec, 0xa681, 0x3eb1, 0xf128, 0x2478, 0xc545, 0x3ef0,
- 0xb9a1, 0x66fc, 0xcfae, 0x3f36, 0xe930, 0xad5d, 0xa0bf, 0x3f8a, 0x4fcf,
- 0x7639, 0xde19, 0x3ff0, 0x9625, 0x962f, 0x25d7};
- #endif
- /* Sine and cosine integrals */
- #ifdef ANSIPROT
- extern double log(double);
- extern double exp(double);
- extern double fabs(double);
- extern double chbevl(double, void *, int);
- #else
- double log(), exp(), fabs(), chbevl();
- #endif
- #define EUL 0.57721566490153286061
- extern double MACHEP, MAXNUM, PIO2;
- int shichi(x, si, ci) double x;
- double *si, *ci;
- {
- double k, z, c, s, a;
- short sign;
- if (x < 0.0) {
- sign = -1;
- x = -x;
- } else
- sign = 0;
- if (x == 0.0) {
- *si = 0.0;
- *ci = -MAXNUM;
- return (0);
- }
- if (x >= 8.0)
- goto chb;
- z = x * x;
- /* Direct power series expansion */
- a = 1.0;
- s = 1.0;
- c = 0.0;
- k = 2.0;
- do {
- a *= z / k;
- c += a / k;
- k += 1.0;
- a /= k;
- s += a / k;
- k += 1.0;
- } while (fabs(a / s) > MACHEP);
- s *= x;
- goto done;
- chb:
- if (x < 18.0) {
- a = (576.0 / x - 52.0) / 10.0;
- k = exp(x) / x;
- s = k * chbevl(a, S1, 22);
- c = k * chbevl(a, C1, 23);
- goto done;
- }
- if (x <= 88.0) {
- a = (6336.0 / x - 212.0) / 70.0;
- k = exp(x) / x;
- s = k * chbevl(a, S2, 23);
- c = k * chbevl(a, C2, 24);
- goto done;
- } else {
- if (sign)
- *si = -MAXNUM;
- else
- *si = MAXNUM;
- *ci = MAXNUM;
- return (0);
- }
- done:
- if (sign)
- s = -s;
- *si = s;
- *ci = EUL + log(x) + c;
- return (0);
- }
|