lib/libm/k0.c

   1 /*                                                      k0.c
   2  *
   3  *      Modified Bessel function, third kind, order zero
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, k0();
  10  *
  11  * y = k0( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns modified Bessel function of the third kind
  18  * of order zero of the argument.
  19  *
  20  * The range is partitioned into the two intervals [0,8] and
  21  * (8, infinity).  Chebyshev polynomial expansions are employed
  22  * in each interval.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  * Tested at 2000 random points between 0 and 8.  Peak absolute
  29  * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    DEC       0, 30        3100       1.3e-16     2.1e-17
  33  *    IEEE      0, 30       30000       1.2e-15     1.6e-16
  34  *
  35  * ERROR MESSAGES:
  36  *
  37  *   message         condition      value returned
  38  *  K0 domain          x <= 0          MAXNUM
  39  *
  40  */
  41 \f/*                                                     k0e()
  42  *
  43  *      Modified Bessel function, third kind, order zero,
  44  *      exponentially scaled
  45  *
  46  *
  47  *
  48  * SYNOPSIS:
  49  *
  50  * double x, y, k0e();
  51  *
  52  * y = k0e( x );
  53  *
  54  *
  55  *
  56  * DESCRIPTION:
  57  *
  58  * Returns exponentially scaled modified Bessel function
  59  * of the third kind of order zero of the argument.
  60  *
  61  *
  62  *
  63  * ACCURACY:
  64  *
  65  *                      Relative error:
  66  * arithmetic   domain     # trials      peak         rms
  67  *    IEEE      0, 30       30000       1.4e-15     1.4e-16
  68  * See k0().
  69  *
  70  */
  71 \f
  72 /*
  73 Cephes Math Library Release 2.8:  June, 2000
  74 Copyright 1984, 1987, 2000 by Stephen L. Moshier
  75 */
  76
  77 #include "mconf.h"
  78
  79 /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
  80  * in the interval [0,2].  The odd order coefficients are all
  81  * zero; only the even order coefficients are listed.
  82  *
  83  * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
  84  */
  85
  86 #ifdef UNK
  87 static double A[] =
  88 {
  89  1.37446543561352307156E-16,
  90  4.25981614279661018399E-14,
  91  1.03496952576338420167E-11,
  92  1.90451637722020886025E-9,
  93  2.53479107902614945675E-7,
  94  2.28621210311945178607E-5,
  95  1.26461541144692592338E-3,
  96  3.59799365153615016266E-2,
  97  3.44289899924628486886E-1,
  98 -5.35327393233902768720E-1
  99 };
 100 #endif
 101
 102 #ifdef DEC
 103 static unsigned short A[] = {
 104 0023036,0073417,0032477,0165673,
 105 0025077,0154126,0016046,0012517,
 106 0027066,0011342,0035211,0005041,
 107 0031002,0160233,0037454,0050224,
 108 0032610,0012747,0037712,0173741,
 109 0034277,0144007,0172147,0162375,
 110 0035645,0140563,0125431,0165626,
 111 0037023,0057662,0125124,0102051,
 112 0037660,0043304,0004411,0166707,
 113 0140011,0005467,0047227,0130370
 114 };
 115 #endif
 116
 117 #ifdef IBMPC
 118 static unsigned short A[] = {
 119 0xfd77,0xe6a7,0xcee1,0x3ca3,
 120 0xc2aa,0xc384,0xfb0a,0x3d27,
 121 0x2144,0x4751,0xc25c,0x3da6,
 122 0x8a13,0x67e5,0x5c13,0x3e20,
 123 0x5efc,0xe7f9,0x02bc,0x3e91,
 124 0xfca0,0xfe8c,0xf900,0x3ef7,
 125 0x3d73,0x7563,0xb82e,0x3f54,
 126 0x9085,0x554a,0x6bf6,0x3fa2,
 127 0x3db9,0x8121,0x08d8,0x3fd6,
 128 0xf61f,0xe9d2,0x2166,0xbfe1
 129 };
 130 #endif
 131
 132 #ifdef MIEEE
 133 static unsigned short A[] = {
 134 0x3ca3,0xcee1,0xe6a7,0xfd77,
 135 0x3d27,0xfb0a,0xc384,0xc2aa,
 136 0x3da6,0xc25c,0x4751,0x2144,
 137 0x3e20,0x5c13,0x67e5,0x8a13,
 138 0x3e91,0x02bc,0xe7f9,0x5efc,
 139 0x3ef7,0xf900,0xfe8c,0xfca0,
 140 0x3f54,0xb82e,0x7563,0x3d73,
 141 0x3fa2,0x6bf6,0x554a,0x9085,
 142 0x3fd6,0x08d8,0x8121,0x3db9,
 143 0xbfe1,0x2166,0xe9d2,0xf61f
 144 };
 145 #endif
 146
 147
 148
 149 /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
 150  * in the inverted interval [2,infinity].
 151  *
 152  * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
 153  */
 154
 155 #ifdef UNK
 156 static double B[] = {
 157  5.30043377268626276149E-18,
 158 -1.64758043015242134646E-17,
 159  5.21039150503902756861E-17,
 160 -1.67823109680541210385E-16,
 161  5.51205597852431940784E-16,
 162 -1.84859337734377901440E-15,
 163  6.34007647740507060557E-15,
 164 -2.22751332699166985548E-14,
 165  8.03289077536357521100E-14,
 166 -2.98009692317273043925E-13,
 167  1.14034058820847496303E-12,
 168 -4.51459788337394416547E-12,
 169  1.85594911495471785253E-11,
 170 -7.95748924447710747776E-11,
 171  3.57739728140030116597E-10,
 172 -1.69753450938905987466E-9,
 173  8.57403401741422608519E-9,
 174 -4.66048989768794782956E-8,
 175  2.76681363944501510342E-7,
 176 -1.83175552271911948767E-6,
 177  1.39498137188764993662E-5,
 178 -1.28495495816278026384E-4,
 179  1.56988388573005337491E-3,
 180 -3.14481013119645005427E-2,
 181  2.44030308206595545468E0
 182 };
 183 #endif
 184
 185 #ifdef DEC
 186 static unsigned short B[] = {
 187 0021703,0106456,0076144,0173406,
 188 0122227,0173144,0116011,0030033,
 189 0022560,0044562,0006506,0067642,
 190 0123101,0076243,0123273,0131013,
 191 0023436,0157713,0056243,0141331,
 192 0124005,0032207,0063726,0164664,
 193 0024344,0066342,0051756,0162300,
 194 0124710,0121365,0154053,0077022,
 195 0025264,0161166,0066246,0077420,
 196 0125647,0141671,0006443,0103212,
 197 0026240,0076431,0077147,0160445,
 198 0126636,0153741,0174002,0105031,
 199 0027243,0040102,0035375,0163073,
 200 0127656,0176256,0113476,0044653,
 201 0030304,0125544,0006377,0130104,
 202 0130751,0047257,0110537,0127324,
 203 0031423,0046400,0014772,0012164,
 204 0132110,0025240,0155247,0112570,
 205 0032624,0105314,0007437,0021574,
 206 0133365,0155243,0174306,0116506,
 207 0034152,0004776,0061643,0102504,
 208 0135006,0136277,0036104,0175023,
 209 0035715,0142217,0162474,0115022,
 210 0137000,0147671,0065177,0134356,
 211 0040434,0026754,0175163,0044070
 212 };
 213 #endif
 214
 215 #ifdef IBMPC
 216 static unsigned short B[] = {
 217 0x9ee1,0xcf8c,0x71a5,0x3c58,
 218 0x2603,0x9381,0xfecc,0xbc72,
 219 0xcdf4,0x41a8,0x092e,0x3c8e,
 220 0x7641,0x74d7,0x2f94,0xbca8,
 221 0x785b,0x6b94,0xdbf9,0x3cc3,
 222 0xdd36,0xecfa,0xa690,0xbce0,
 223 0xdc98,0x4a7d,0x8d9c,0x3cfc,
 224 0x6fc2,0xbb05,0x145e,0xbd19,
 225 0xcfe2,0xcd94,0x9c4e,0x3d36,
 226 0x70d1,0x21a4,0xf877,0xbd54,
 227 0xfc25,0x2fcc,0x0fa3,0x3d74,
 228 0x5143,0x3f00,0xdafc,0xbd93,
 229 0xbcc7,0x475f,0x6808,0x3db4,
 230 0xc935,0xd2e7,0xdf95,0xbdd5,
 231 0xf608,0x819f,0x956c,0x3df8,
 232 0xf5db,0xf22b,0x29d5,0xbe1d,
 233 0x428e,0x033f,0x69a0,0x3e42,
 234 0xf2af,0x1b54,0x0554,0xbe69,
 235 0xe46f,0x81e3,0x9159,0x3e92,
 236 0xd3a9,0x7f18,0xbb54,0xbebe,
 237 0x70a9,0xcc74,0x413f,0x3eed,
 238 0x9f42,0xe788,0xd797,0xbf20,
 239 0x9342,0xfca7,0xb891,0x3f59,
 240 0xf71e,0x2d4f,0x19f7,0xbfa0,
 241 0x6907,0x9f4e,0x85bd,0x4003
 242 };
 243 #endif
 244
 245 #ifdef MIEEE
 246 static unsigned short B[] = {
 247 0x3c58,0x71a5,0xcf8c,0x9ee1,
 248 0xbc72,0xfecc,0x9381,0x2603,
 249 0x3c8e,0x092e,0x41a8,0xcdf4,
 250 0xbca8,0x2f94,0x74d7,0x7641,
 251 0x3cc3,0xdbf9,0x6b94,0x785b,
 252 0xbce0,0xa690,0xecfa,0xdd36,
 253 0x3cfc,0x8d9c,0x4a7d,0xdc98,
 254 0xbd19,0x145e,0xbb05,0x6fc2,
 255 0x3d36,0x9c4e,0xcd94,0xcfe2,
 256 0xbd54,0xf877,0x21a4,0x70d1,
 257 0x3d74,0x0fa3,0x2fcc,0xfc25,
 258 0xbd93,0xdafc,0x3f00,0x5143,
 259 0x3db4,0x6808,0x475f,0xbcc7,
 260 0xbdd5,0xdf95,0xd2e7,0xc935,
 261 0x3df8,0x956c,0x819f,0xf608,
 262 0xbe1d,0x29d5,0xf22b,0xf5db,
 263 0x3e42,0x69a0,0x033f,0x428e,
 264 0xbe69,0x0554,0x1b54,0xf2af,
 265 0x3e92,0x9159,0x81e3,0xe46f,
 266 0xbebe,0xbb54,0x7f18,0xd3a9,
 267 0x3eed,0x413f,0xcc74,0x70a9,
 268 0xbf20,0xd797,0xe788,0x9f42,
 269 0x3f59,0xb891,0xfca7,0x9342,
 270 0xbfa0,0x19f7,0x2d4f,0xf71e,
 271 0x4003,0x85bd,0x9f4e,0x6907
 272 };
 273 #endif
 274
 275 /*                                                      k0.c    */
 276 #ifdef ANSIPROT
 277 extern double chbevl ( double, void *, int );
 278 extern double exp ( double );
 279 extern double i0 ( double );
 280 extern double log ( double );
 281 extern double sqrt ( double );
 282 #else
 283 double chbevl(), exp(), i0(), log(), sqrt();
 284 #endif
 285 extern double PI;
 286 extern double MAXNUM;
 287
 288 double k0(x)
 289 double x;
 290 {
 291 double y, z;
 292
 293 if( x <= 0.0 )
 294         {
 295         mtherr( "k0", DOMAIN );
 296         return( MAXNUM );
 297         }
 298
 299 if( x <= 2.0 )
 300         {
 301         y = x * x - 2.0;
 302         y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x);
 303         return( y );
 304         }
 305 z = 8.0/x - 2.0;
 306 y = exp(-x) * chbevl( z, B, 25 ) / sqrt(x);
 307 return(y);
 308 }
 309
 310
 311
 312
 313 double k0e( x )
 314 double x;
 315 {
 316 double y;
 317
 318 if( x <= 0.0 )
 319         {
 320         mtherr( "k0e", DOMAIN );
 321         return( MAXNUM );
 322         }
 323
 324 if( x <= 2.0 )
 325         {
 326         y = x * x - 2.0;
 327         y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x);
 328         return( y * exp(x) );
 329         }
 330
 331 y = chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x);
 332 return(y);
 333 }