From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/k0.c | 333 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 333 insertions(+) create mode 100644 libm/double/k0.c (limited to 'libm/double/k0.c') diff --git a/libm/double/k0.c b/libm/double/k0.c new file mode 100644 index 000000000..7d09cb4a1 --- /dev/null +++ b/libm/double/k0.c @@ -0,0 +1,333 @@ +/* k0.c + * + * Modified Bessel function, third kind, order zero + * + * + * + * SYNOPSIS: + * + * double x, y, k0(); + * + * y = k0( x ); + * + * + * + * DESCRIPTION: + * + * Returns modified Bessel function of the third kind + * of order zero of the argument. + * + * The range is partitioned into the two intervals [0,8] and + * (8, infinity). Chebyshev polynomial expansions are employed + * in each interval. + * + * + * + * ACCURACY: + * + * Tested at 2000 random points between 0 and 8. Peak absolute + * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. + * Relative error: + * arithmetic domain # trials peak rms + * DEC 0, 30 3100 1.3e-16 2.1e-17 + * IEEE 0, 30 30000 1.2e-15 1.6e-16 + * + * ERROR MESSAGES: + * + * message condition value returned + * K0 domain x <= 0 MAXNUM + * + */ + /* k0e() + * + * Modified Bessel function, third kind, order zero, + * exponentially scaled + * + * + * + * SYNOPSIS: + * + * double x, y, k0e(); + * + * y = k0e( x ); + * + * + * + * DESCRIPTION: + * + * Returns exponentially scaled modified Bessel function + * of the third kind of order zero of the argument. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0, 30 30000 1.4e-15 1.4e-16 + * See k0(). + * + */ + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 2000 by Stephen L. Moshier +*/ + +#include + +/* Chebyshev coefficients for K0(x) + log(x/2) I0(x) + * in the interval [0,2]. The odd order coefficients are all + * zero; only the even order coefficients are listed. + * + * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL. + */ + +#ifdef UNK +static double A[] = +{ + 1.37446543561352307156E-16, + 4.25981614279661018399E-14, + 1.03496952576338420167E-11, + 1.90451637722020886025E-9, + 2.53479107902614945675E-7, + 2.28621210311945178607E-5, + 1.26461541144692592338E-3, + 3.59799365153615016266E-2, + 3.44289899924628486886E-1, +-5.35327393233902768720E-1 +}; +#endif + +#ifdef DEC +static unsigned short A[] = { +0023036,0073417,0032477,0165673, +0025077,0154126,0016046,0012517, +0027066,0011342,0035211,0005041, +0031002,0160233,0037454,0050224, +0032610,0012747,0037712,0173741, +0034277,0144007,0172147,0162375, +0035645,0140563,0125431,0165626, +0037023,0057662,0125124,0102051, +0037660,0043304,0004411,0166707, +0140011,0005467,0047227,0130370 +}; +#endif + +#ifdef IBMPC +static unsigned short A[] = { +0xfd77,0xe6a7,0xcee1,0x3ca3, +0xc2aa,0xc384,0xfb0a,0x3d27, +0x2144,0x4751,0xc25c,0x3da6, +0x8a13,0x67e5,0x5c13,0x3e20, +0x5efc,0xe7f9,0x02bc,0x3e91, +0xfca0,0xfe8c,0xf900,0x3ef7, +0x3d73,0x7563,0xb82e,0x3f54, +0x9085,0x554a,0x6bf6,0x3fa2, +0x3db9,0x8121,0x08d8,0x3fd6, +0xf61f,0xe9d2,0x2166,0xbfe1 +}; +#endif + +#ifdef MIEEE +static unsigned short A[] = { +0x3ca3,0xcee1,0xe6a7,0xfd77, +0x3d27,0xfb0a,0xc384,0xc2aa, +0x3da6,0xc25c,0x4751,0x2144, +0x3e20,0x5c13,0x67e5,0x8a13, +0x3e91,0x02bc,0xe7f9,0x5efc, +0x3ef7,0xf900,0xfe8c,0xfca0, +0x3f54,0xb82e,0x7563,0x3d73, +0x3fa2,0x6bf6,0x554a,0x9085, +0x3fd6,0x08d8,0x8121,0x3db9, +0xbfe1,0x2166,0xe9d2,0xf61f +}; +#endif + + + +/* Chebyshev coefficients for exp(x) sqrt(x) K0(x) + * in the inverted interval [2,infinity]. + * + * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2). + */ + +#ifdef UNK +static double B[] = { + 5.30043377268626276149E-18, +-1.64758043015242134646E-17, + 5.21039150503902756861E-17, +-1.67823109680541210385E-16, + 5.51205597852431940784E-16, +-1.84859337734377901440E-15, + 6.34007647740507060557E-15, +-2.22751332699166985548E-14, + 8.03289077536357521100E-14, +-2.98009692317273043925E-13, + 1.14034058820847496303E-12, +-4.51459788337394416547E-12, + 1.85594911495471785253E-11, +-7.95748924447710747776E-11, + 3.57739728140030116597E-10, +-1.69753450938905987466E-9, + 8.57403401741422608519E-9, +-4.66048989768794782956E-8, + 2.76681363944501510342E-7, +-1.83175552271911948767E-6, + 1.39498137188764993662E-5, +-1.28495495816278026384E-4, + 1.56988388573005337491E-3, +-3.14481013119645005427E-2, + 2.44030308206595545468E0 +}; +#endif + +#ifdef DEC +static unsigned short B[] = { +0021703,0106456,0076144,0173406, +0122227,0173144,0116011,0030033, +0022560,0044562,0006506,0067642, +0123101,0076243,0123273,0131013, +0023436,0157713,0056243,0141331, +0124005,0032207,0063726,0164664, +0024344,0066342,0051756,0162300, +0124710,0121365,0154053,0077022, +0025264,0161166,0066246,0077420, +0125647,0141671,0006443,0103212, +0026240,0076431,0077147,0160445, +0126636,0153741,0174002,0105031, +0027243,0040102,0035375,0163073, +0127656,0176256,0113476,0044653, +0030304,0125544,0006377,0130104, +0130751,0047257,0110537,0127324, +0031423,0046400,0014772,0012164, +0132110,0025240,0155247,0112570, +0032624,0105314,0007437,0021574, +0133365,0155243,0174306,0116506, +0034152,0004776,0061643,0102504, +0135006,0136277,0036104,0175023, +0035715,0142217,0162474,0115022, +0137000,0147671,0065177,0134356, +0040434,0026754,0175163,0044070 +}; +#endif + +#ifdef IBMPC +static unsigned short B[] = { +0x9ee1,0xcf8c,0x71a5,0x3c58, +0x2603,0x9381,0xfecc,0xbc72, +0xcdf4,0x41a8,0x092e,0x3c8e, +0x7641,0x74d7,0x2f94,0xbca8, +0x785b,0x6b94,0xdbf9,0x3cc3, +0xdd36,0xecfa,0xa690,0xbce0, +0xdc98,0x4a7d,0x8d9c,0x3cfc, +0x6fc2,0xbb05,0x145e,0xbd19, +0xcfe2,0xcd94,0x9c4e,0x3d36, +0x70d1,0x21a4,0xf877,0xbd54, +0xfc25,0x2fcc,0x0fa3,0x3d74, +0x5143,0x3f00,0xdafc,0xbd93, +0xbcc7,0x475f,0x6808,0x3db4, +0xc935,0xd2e7,0xdf95,0xbdd5, +0xf608,0x819f,0x956c,0x3df8, +0xf5db,0xf22b,0x29d5,0xbe1d, +0x428e,0x033f,0x69a0,0x3e42, +0xf2af,0x1b54,0x0554,0xbe69, +0xe46f,0x81e3,0x9159,0x3e92, +0xd3a9,0x7f18,0xbb54,0xbebe, +0x70a9,0xcc74,0x413f,0x3eed, +0x9f42,0xe788,0xd797,0xbf20, +0x9342,0xfca7,0xb891,0x3f59, +0xf71e,0x2d4f,0x19f7,0xbfa0, +0x6907,0x9f4e,0x85bd,0x4003 +}; +#endif + +#ifdef MIEEE +static unsigned short B[] = { +0x3c58,0x71a5,0xcf8c,0x9ee1, +0xbc72,0xfecc,0x9381,0x2603, +0x3c8e,0x092e,0x41a8,0xcdf4, +0xbca8,0x2f94,0x74d7,0x7641, +0x3cc3,0xdbf9,0x6b94,0x785b, +0xbce0,0xa690,0xecfa,0xdd36, +0x3cfc,0x8d9c,0x4a7d,0xdc98, +0xbd19,0x145e,0xbb05,0x6fc2, +0x3d36,0x9c4e,0xcd94,0xcfe2, +0xbd54,0xf877,0x21a4,0x70d1, +0x3d74,0x0fa3,0x2fcc,0xfc25, +0xbd93,0xdafc,0x3f00,0x5143, +0x3db4,0x6808,0x475f,0xbcc7, +0xbdd5,0xdf95,0xd2e7,0xc935, +0x3df8,0x956c,0x819f,0xf608, +0xbe1d,0x29d5,0xf22b,0xf5db, +0x3e42,0x69a0,0x033f,0x428e, +0xbe69,0x0554,0x1b54,0xf2af, +0x3e92,0x9159,0x81e3,0xe46f, +0xbebe,0xbb54,0x7f18,0xd3a9, +0x3eed,0x413f,0xcc74,0x70a9, +0xbf20,0xd797,0xe788,0x9f42, +0x3f59,0xb891,0xfca7,0x9342, +0xbfa0,0x19f7,0x2d4f,0xf71e, +0x4003,0x85bd,0x9f4e,0x6907 +}; +#endif + +/* k0.c */ +#ifdef ANSIPROT +extern double chbevl ( double, void *, int ); +extern double exp ( double ); +extern double i0 ( double ); +extern double log ( double ); +extern double sqrt ( double ); +#else +double chbevl(), exp(), i0(), log(), sqrt(); +#endif +extern double PI; +extern double MAXNUM; + +double k0(x) +double x; +{ +double y, z; + +if( x <= 0.0 ) + { + mtherr( "k0", DOMAIN ); + return( MAXNUM ); + } + +if( x <= 2.0 ) + { + y = x * x - 2.0; + y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x); + return( y ); + } +z = 8.0/x - 2.0; +y = exp(-x) * chbevl( z, B, 25 ) / sqrt(x); +return(y); +} + + + + +double k0e( x ) +double x; +{ +double y; + +if( x <= 0.0 ) + { + mtherr( "k0e", DOMAIN ); + return( MAXNUM ); + } + +if( x <= 2.0 ) + { + y = x * x - 2.0; + y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x); + return( y * exp(x) ); + } + +y = chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x); +return(y); +} -- cgit v1.2.3