lib/libm/j0f.c

   1 /*                                                      j0f.c
   2  *
   3  *      Bessel function of order zero
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, y, j0f();
  10  *
  11  * y = j0f( x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns Bessel function of order zero of the argument.
  18  *
  19  * The domain is divided into the intervals [0, 2] and
  20  * (2, infinity). In the first interval the following polynomial
  21  * approximation is used:
  22  *
  23  *
  24  *        2         2         2
  25  * (w - r  ) (w - r  ) (w - r  ) P(w)
  26  *       1         2         3
  27  *
  28  *            2
  29  * where w = x  and the three r's are zeros of the function.
  30  *
  31  * In the second interval, the modulus and phase are approximated
  32  * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
  33  * and Phase(x) = x + 1/x R(1/x^2) - pi/4.  The function is
  34  *
  35  *   j0(x) = Modulus(x) cos( Phase(x) ).
  36  *
  37  *
  38  *
  39  * ACCURACY:
  40  *
  41  *                      Absolute error:
  42  * arithmetic   domain     # trials      peak         rms
  43  *    IEEE      0, 2        100000      1.3e-7      3.6e-8
  44  *    IEEE      2, 32       100000      1.9e-7      5.4e-8
  45  *
  46  */
  47 \f/*                                                     y0f.c
  48  *
  49  *      Bessel function of the second kind, order zero
  50  *
  51  *
  52  *
  53  * SYNOPSIS:
  54  *
  55  * float x, y, y0f();
  56  *
  57  * y = y0f( x );
  58  *
  59  *
  60  *
  61  * DESCRIPTION:
  62  *
  63  * Returns Bessel function of the second kind, of order
  64  * zero, of the argument.
  65  *
  66  * The domain is divided into the intervals [0, 2] and
  67  * (2, infinity). In the first interval a rational approximation
  68  * R(x) is employed to compute
  69  *
  70  *                  2         2         2
  71  * y0(x)  =  (w - r  ) (w - r  ) (w - r  ) R(x)  +  2/pi ln(x) j0(x).
  72  *                 1         2         3
  73  *
  74  * Thus a call to j0() is required.  The three zeros are removed
  75  * from R(x) to improve its numerical stability.
  76  *
  77  * In the second interval, the modulus and phase are approximated
  78  * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
  79  * and Phase(x) = x + 1/x S(1/x^2) - pi/4.  Then the function is
  80  *
  81  *   y0(x) = Modulus(x) sin( Phase(x) ).
  82  *
  83  *
  84  *
  85  *
  86  * ACCURACY:
  87  *
  88  *  Absolute error, when y0(x) < 1; else relative error:
  89  *
  90  * arithmetic   domain     # trials      peak         rms
  91  *    IEEE      0,  2       100000      2.4e-7      3.4e-8
  92  *    IEEE      2, 32       100000      1.8e-7      5.3e-8
  93  *
  94  */
  95 \f
  96 /*
  97 Cephes Math Library Release 2.2:  June, 1992
  98 Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
  99 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 100 */
 101
 102
 103 #include "mconf.h"
 104
 105 static float MO[8] = {
 106 -6.838999669318810E-002f,
 107  1.864949361379502E-001f,
 108 -2.145007480346739E-001f,
 109  1.197549369473540E-001f,
 110 -3.560281861530129E-003f,
 111 -4.969382655296620E-002f,
 112 -3.355424622293709E-006f,
 113  7.978845717621440E-001f
 114 };
 115
 116 static float PH[8] = {
 117  3.242077816988247E+001f,
 118 -3.630592630518434E+001f,
 119  1.756221482109099E+001f,
 120 -4.974978466280903E+000f,
 121  1.001973420681837E+000f,
 122 -1.939906941791308E-001f,
 123  6.490598792654666E-002f,
 124 -1.249992184872738E-001f
 125 };
 126
 127 static float YP[5] = {
 128  9.454583683980369E-008f,
 129 -9.413212653797057E-006f,
 130  5.344486707214273E-004f,
 131 -1.584289289821316E-002f,
 132  1.707584643733568E-001f
 133 };
 134
 135 float YZ1 =  0.43221455686510834878f;
 136 float YZ2 = 22.401876406482861405f;
 137 float YZ3 = 64.130620282338755553f;
 138
 139 static float DR1 =  5.78318596294678452118f;
 140 /*
 141 static float DR2 = 30.4712623436620863991;
 142 static float DR3 = 74.887006790695183444889;
 143 */
 144
 145 static float JP[5] = {
 146 -6.068350350393235E-008f,
 147  6.388945720783375E-006f,
 148 -3.969646342510940E-004f,
 149  1.332913422519003E-002f,
 150 -1.729150680240724E-001f
 151 };
 152 extern float PIO4F;
 153
 154
 155 #ifdef ANSIC
 156 float polevlf(float, float *, int);
 157 float logf(float), sinf(float), cosf(float), sqrtf(float);
 158
 159 float j0f( float xx )
 160 #else
 161 float polevlf(), logf(), sinf(), cosf(), sqrtf();
 162
 163 float j0f(xx)
 164 double xx;
 165 #endif
 166 {
 167 float x, w, z, p, q, xn;
 168
 169
 170 if( xx < 0 )
 171         x = -xx;
 172 else
 173         x = xx;
 174
 175 if( x <= 2.0f )
 176         {
 177         z = x * x;
 178         if( x < 1.0e-3f )
 179                 return( 1.0f - 0.25f*z );
 180
 181         p = (z-DR1) * polevlf( z, JP, 4);
 182         return( p );
 183         }
 184
 185 q = 1.0f/x;
 186 w = sqrtf(q);
 187
 188 p = w * polevlf( q, MO, 7);
 189 w = q*q;
 190 xn = q * polevlf( w, PH, 7) - PIO4F;
 191 p = p * cosf(xn + x);
 192 return(p);
 193 }
 194 \f
 195 /*                                                      y0() 2  */
 196 /* Bessel function of second kind, order zero   */
 197
 198 /* Rational approximation coefficients YP[] are used for x < 6.5.
 199  * The function computed is  y0(x)  -  2 ln(x) j0(x) / pi,
 200  * whose value at x = 0 is  2 * ( log(0.5) + EUL ) / pi
 201  * = 0.073804295108687225 , EUL is Euler's constant.
 202  */
 203
 204 static float TWOOPI =  0.636619772367581343075535f; /* 2/pi */
 205 extern float MAXNUMF;
 206
 207 #ifdef ANSIC
 208 float y0f( float xx )
 209 #else
 210 float y0f(xx)
 211 double xx;
 212 #endif
 213 {
 214 float x, w, z, p, q, xn;
 215
 216
 217 x = xx;
 218 if( x <= 2.0f )
 219         {
 220         if( x <= 0.0f )
 221                 {
 222                 mtherr( "y0f", DOMAIN );
 223                 return( -MAXNUMF );
 224                 }
 225         z = x * x;
 226 /*      w = (z-YZ1)*(z-YZ2)*(z-YZ3) * polevlf( z, YP, 4);*/
 227         w = (z-YZ1) * polevlf( z, YP, 4);
 228         w += TWOOPI * logf(x) * j0f(x);
 229         return( w );
 230         }
 231
 232 q = 1.0f/x;
 233 w = sqrtf(q);
 234
 235 p = w * polevlf( q, MO, 7);
 236 w = q*q;
 237 xn = q * polevlf( w, PH, 7) - PIO4F;
 238 p = p * sinf(xn + x);
 239 return( p );
 240 }