Quelle  k0.c

/* 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 "mconf.h"

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