lib/libm/i1f.c

   1 /*                                                      i1f.c
   2  *
   3  *      Modified Bessel function of order one
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, i1f();
  10  *
  11  * y = i1f( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns modified Bessel function of order one of the
  18  * argument.
  19  *
  20  * The function is defined as i1(x) = -i j1( ix ).
  21  *
  22  * The range is partitioned into the two intervals [0,8] and
  23  * (8, infinity).  Chebyshev polynomial expansions are employed
  24  * in each interval.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    IEEE      0, 30       100000      1.5e-6      1.6e-7
  33  *
  34  *
  35  */
  36 \f/*                                                     i1ef.c
  37  *
  38  *      Modified Bessel function of order one,
  39  *      exponentially scaled
  40  *
  41  *
  42  *
  43  * SYNOPSIS:
  44  *
  45  * float x, y, i1ef();
  46  *
  47  * y = i1ef( x );
  48  *
  49  *
  50  *
  51  * DESCRIPTION:
  52  *
  53  * Returns exponentially scaled modified Bessel function
  54  * of order one of the argument.
  55  *
  56  * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
  57  *
  58  *
  59  *
  60  * ACCURACY:
  61  *
  62  *                      Relative error:
  63  * arithmetic   domain     # trials      peak         rms
  64  *    IEEE      0, 30       30000       1.5e-6      1.5e-7
  65  * See i1().
  66  *
  67  */
  68 \f
  69 /*                                                      i1.c 2          */
  70
  71
  72 /*
  73 Cephes Math Library Release 2.0:  March, 1987
  74 Copyright 1985, 1987 by Stephen L. Moshier
  75 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  76 */
  77
  78 #include "mconf.h"
  79
  80 /* Chebyshev coefficients for exp(-x) I1(x) / x
  81  * in the interval [0,8].
  82  *
  83  * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
  84  */
  85
  86 static float A[] =
  87 {
  88  9.38153738649577178388E-9f,
  89 -4.44505912879632808065E-8f,
  90  2.00329475355213526229E-7f,
  91 -8.56872026469545474066E-7f,
  92  3.47025130813767847674E-6f,
  93 -1.32731636560394358279E-5f,
  94  4.78156510755005422638E-5f,
  95 -1.61760815825896745588E-4f,
  96  5.12285956168575772895E-4f,
  97 -1.51357245063125314899E-3f,
  98  4.15642294431288815669E-3f,
  99 -1.05640848946261981558E-2f,
 100  2.47264490306265168283E-2f,
 101 -5.29459812080949914269E-2f,
 102  1.02643658689847095384E-1f,
 103 -1.76416518357834055153E-1f,
 104  2.52587186443633654823E-1f
 105 };
 106
 107
 108 /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
 109  * in the inverted interval [8,infinity].
 110  *
 111  * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
 112  */
 113
 114 static float B[] =
 115 {
 116 -3.83538038596423702205E-9f,
 117 -2.63146884688951950684E-8f,
 118 -2.51223623787020892529E-7f,
 119 -3.88256480887769039346E-6f,
 120 -1.10588938762623716291E-4f,
 121 -9.76109749136146840777E-3f,
 122  7.78576235018280120474E-1f
 123 };
 124 \f
 125 /*                                                      i1.c    */
 126
 127 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
 128
 129 #ifdef ANSIC
 130 float chbevlf(float, float *, int);
 131 float expf(float), sqrtf(float);
 132 #else
 133 float chbevlf(), expf(), sqrtf();
 134 #endif
 135
 136
 137 #ifdef ANSIC
 138 float i1f(float xx)
 139 #else
 140 float i1f(xx)
 141 double xx;
 142 #endif
 143 {
 144 float x, y, z;
 145
 146 x = xx;
 147 z = fabsf(x);
 148 if( z <= 8.0f )
 149         {
 150         y = 0.5f*z - 2.0f;
 151         z = chbevlf( y, A, 17 ) * z * expf(z);
 152         }
 153 else
 154         {
 155         z = expf(z) * chbevlf( 32.0f/z - 2.0f, B, 7 ) / sqrtf(z);
 156         }
 157 if( x < 0.0f )
 158         z = -z;
 159 return( z );
 160 }
 161 \f
 162 /*                                                      i1e()   */
 163
 164 #ifdef ANSIC
 165 float i1ef( float xx )
 166 #else
 167 float i1ef( xx )
 168 double xx;
 169 #endif
 170 {
 171 float x, y, z;
 172
 173 x = xx;
 174 z = fabsf(x);
 175 if( z <= 8.0f )
 176         {
 177         y = 0.5f*z - 2.0f;
 178         z = chbevlf( y, A, 17 ) * z;
 179         }
 180 else
 181         {
 182         z = chbevlf( 32.0f/z - 2.0f, B, 7 ) / sqrtf(z);
 183         }
 184 if( x < 0.0f )
 185         z = -z;
 186 return( z );
 187 }