lib/libm/k1f.c

   1 /*                                                      k1f.c
   2  *
   3  *      Modified Bessel function, third kind, order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, k1f();
  10  *
  11  * y = k1f( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Computes the modified Bessel function of the third kind
  18  * of order one of the argument.
  19  *
  20  * The range is partitioned into the two intervals [0,2] and
  21  * (2, infinity).  Chebyshev polynomial expansions are employed
  22  * in each interval.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *                      Relative error:
  29  * arithmetic   domain     # trials      peak         rms
  30  *    IEEE      0, 30       30000       4.6e-7      7.6e-8
  31  *
  32  * ERROR MESSAGES:
  33  *
  34  *   message         condition      value returned
  35  * k1 domain          x <= 0          MAXNUM
  36  *
  37  */
  38 \f/*                                                     k1ef.c
  39  *
  40  *      Modified Bessel function, third kind, order one,
  41  *      exponentially scaled
  42  *
  43  *
  44  *
  45  * SYNOPSIS:
  46  *
  47  * float x, y, k1ef();
  48  *
  49  * y = k1ef( x );
  50  *
  51  *
  52  *
  53  * DESCRIPTION:
  54  *
  55  * Returns exponentially scaled modified Bessel function
  56  * of the third kind of order one of the argument:
  57  *
  58  *      k1e(x) = exp(x) * k1(x).
  59  *
  60  *
  61  *
  62  * ACCURACY:
  63  *
  64  *                      Relative error:
  65  * arithmetic   domain     # trials      peak         rms
  66  *    IEEE      0, 30       30000       4.9e-7      6.7e-8
  67  * See k1().
  68  *
  69  */
  70 \f
  71 /*
  72 Cephes Math Library Release 2.2: June, 1992
  73 Copyright 1984, 1987, 1992 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 x(K1(x) - log(x/2) I1(x))
  80  * in the interval [0,2].
  81  *
  82  * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
  83  */
  84
  85 #define MINNUMF 6.0e-39
  86 static float A[] =
  87 {
  88 -2.21338763073472585583E-8f,
  89 -2.43340614156596823496E-6f,
  90 -1.73028895751305206302E-4f,
  91 -6.97572385963986435018E-3f,
  92 -1.22611180822657148235E-1f,
  93 -3.53155960776544875667E-1f,
  94  1.52530022733894777053E0f
  95 };
  96
  97
  98
  99
 100 /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
 101  * in the interval [2,infinity].
 102  *
 103  * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
 104  */
 105
 106 static float B[] =
 107 {
 108  2.01504975519703286596E-9f,
 109 -1.03457624656780970260E-8f,
 110  5.74108412545004946722E-8f,
 111 -3.50196060308781257119E-7f,
 112  2.40648494783721712015E-6f,
 113 -1.93619797416608296024E-5f,
 114  1.95215518471351631108E-4f,
 115 -2.85781685962277938680E-3f,
 116  1.03923736576817238437E-1f,
 117  2.72062619048444266945E0f
 118 };
 119
 120
 121
 122 extern float MAXNUMF;
 123 #ifdef ANSIC
 124 float chbevlf(float, float *, int);
 125 float expf(float), i1f(float), logf(float), sqrtf(float);
 126 #else
 127 float chbevlf(), expf(), i1f(), logf(), sqrtf();
 128 #endif
 129
 130 #ifdef ANSIC
 131 float k1f(float xx)
 132 #else
 133 float k1f(xx)
 134 double xx;
 135 #endif
 136 {
 137 float x, y;
 138
 139 x = xx;
 140 if( x <= MINNUMF )
 141         {
 142         mtherr( "k1f", DOMAIN );
 143         return( MAXNUMF );
 144         }
 145
 146 if( x <= 2.0f )
 147         {
 148         y = x * x - 2.0f;
 149         y =  logf( 0.5f * x ) * i1f(x)  +  chbevlf( y, A, 7 ) / x;
 150         return( y );
 151         }
 152
 153 return(  expf(-x) * chbevlf( 8.0f/x - 2.0f, B, 10 ) / sqrtf(x) );
 154
 155 }
 156
 157
 158
 159 #ifdef ANSIC
 160 float k1ef( float xx )
 161 #else
 162 float k1ef( xx )
 163 double xx;
 164 #endif
 165 {
 166 float x, y;
 167
 168 x = xx;
 169 if( x <= 0.0f )
 170         {
 171         mtherr( "k1ef", DOMAIN );
 172         return( MAXNUMF );
 173         }
 174
 175 if( x <= 2.0f )
 176         {
 177         y = x * x - 2.0f;
 178         y =  logf( 0.5f * x ) * i1f(x)  +  chbevlf( y, A, 7 ) / x;
 179         return( y * expf(x) );
 180         }
 181
 182 return(  chbevlf( 8.0f/x - 2.0f, B, 10 ) / sqrtf(x) );
 183
 184 }