lib/libm/k0f.c

   1 /*                                                      k0f.c
   2  *
   3  *      Modified Bessel function, third kind, order zero
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, k0f();
  10  *
  11  * y = k0f( 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  *    IEEE      0, 30       30000       7.8e-7      8.5e-8
  33  *
  34  * ERROR MESSAGES:
  35  *
  36  *   message         condition      value returned
  37  *  K0 domain          x <= 0          MAXNUM
  38  *
  39  */
  40 \f/*                                                     k0ef()
  41  *
  42  *      Modified Bessel function, third kind, order zero,
  43  *      exponentially scaled
  44  *
  45  *
  46  *
  47  * SYNOPSIS:
  48  *
  49  * float x, y, k0ef();
  50  *
  51  * y = k0ef( x );
  52  *
  53  *
  54  *
  55  * DESCRIPTION:
  56  *
  57  * Returns exponentially scaled modified Bessel function
  58  * of the third kind of order zero of the argument.
  59  *
  60  *
  61  *
  62  * ACCURACY:
  63  *
  64  *                      Relative error:
  65  * arithmetic   domain     # trials      peak         rms
  66  *    IEEE      0, 30       30000       8.1e-7      7.8e-8
  67  * See k0().
  68  *
  69  */
  70 \f
  71 /*
  72 Cephes Math Library Release 2.0:  April, 1987
  73 Copyright 1984, 1987 by Stephen L. Moshier
  74 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  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 static float A[] =
  87 {
  88  1.90451637722020886025E-9f,
  89  2.53479107902614945675E-7f,
  90  2.28621210311945178607E-5f,
  91  1.26461541144692592338E-3f,
  92  3.59799365153615016266E-2f,
  93  3.44289899924628486886E-1f,
  94 -5.35327393233902768720E-1f
  95 };
  96
  97
  98
  99 /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
 100  * in the inverted interval [2,infinity].
 101  *
 102  * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
 103  */
 104
 105 static float B[] = {
 106 -1.69753450938905987466E-9f,
 107  8.57403401741422608519E-9f,
 108 -4.66048989768794782956E-8f,
 109  2.76681363944501510342E-7f,
 110 -1.83175552271911948767E-6f,
 111  1.39498137188764993662E-5f,
 112 -1.28495495816278026384E-4f,
 113  1.56988388573005337491E-3f,
 114 -3.14481013119645005427E-2f,
 115  2.44030308206595545468E0f
 116 };
 117
 118 /*                                                      k0.c    */
 119
 120 extern float MAXNUMF;
 121
 122 #ifdef ANSIC
 123 float chbevlf(float, float *, int);
 124 float expf(float), i0f(float), logf(float), sqrtf(float);
 125 #else
 126 float chbevlf(), expf(), i0f(), logf(), sqrtf();
 127 #endif
 128
 129
 130 #ifdef ANSIC
 131 float k0f( float xx )
 132 #else
 133 float k0f(xx)
 134 double xx;
 135 #endif
 136 {
 137 float x, y, z;
 138
 139 x = xx;
 140 if( x <= 0.0f )
 141         {
 142         mtherr( "k0f", DOMAIN );
 143         return( MAXNUMF );
 144         }
 145
 146 if( x <= 2.0f )
 147         {
 148         y = x * x - 2.0f;
 149         y = chbevlf( y, A, 7 ) - logf( 0.5f * x ) * i0f(x);
 150         return( y );
 151         }
 152 z = 8.0f/x - 2.0f;
 153 y = expf(-x) * chbevlf( z, B, 10 ) / sqrtf(x);
 154 return(y);
 155 }
 156
 157
 158
 159 #ifdef ANSIC
 160 float k0ef( float xx )
 161 #else
 162 float k0ef( xx )
 163 double xx;
 164 #endif
 165 {
 166 float x, y;
 167
 168
 169 x = xx;
 170 if( x <= 0.0f )
 171         {
 172         mtherr( "k0ef", DOMAIN );
 173         return( MAXNUMF );
 174         }
 175
 176 if( x <= 2.0f )
 177         {
 178         y = x * x - 2.0f;
 179         y = chbevlf( y, A, 7 ) - logf( 0.5f * x ) * i0f(x);
 180         return( y * expf(x) );
 181         }
 182
 183 y = chbevlf( 8.0f/x - 2.0f, B, 10 ) / sqrtf(x);
 184 return(y);
 185 }