lib/libm/struvef.c

   1 /*                                                      struvef.c
   2  *
   3  *      Struve function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float v, x, y, struvef();
  10  *
  11  * y = struvef( v, x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Computes the Struve function Hv(x) of order v, argument x.
  18  * Negative x is rejected unless v is an integer.
  19  *
  20  * This module also contains the hypergeometric functions 1F2
  21  * and 3F0 and a routine for the Bessel function Yv(x) with
  22  * noninteger v.
  23  *
  24  *
  25  *
  26  * ACCURACY:
  27  *
  28  *  v varies from 0 to 10.
  29  *    Absolute error (relative error when |Hv(x)| > 1):
  30  * arithmetic   domain     # trials      peak         rms
  31  *    IEEE      -10,10      100000      9.0e-5      4.0e-6
  32  *
  33  */
  34 \f
  35
  36 /*
  37 Cephes Math Library Release 2.2:  July, 1992
  38 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  39 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  40 */
  41
  42 #include "mconf.h"
  43 #define DEBUG 0
  44
  45 extern float MACHEPF, MAXNUMF, PIF;
  46
  47 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  48
  49 #ifdef ANSIC
  50 float gammaf(float), powf(float, float), sqrtf(float);
  51 float yvf(float, float);
  52 float floorf(float), ynf(int, float);
  53 float jvf(float, float);
  54 float sinf(float), cosf(float);
  55 #else
  56 float gammaf(), powf(), sqrtf(), yvf();
  57 float floorf(), ynf(), jvf(), sinf(), cosf();
  58 #endif
  59
  60 #ifdef ANSIC
  61 float onef2f( float aa, float bb, float cc, float xx, float *err )
  62 #else
  63 float onef2f( aa, bb, cc, xx, err )
  64 double aa, bb, cc, xx;
  65 float *err;
  66 #endif
  67 {
  68 float a, b, c, x, n, a0, sum, t;
  69 float an, bn, cn, max, z;
  70
  71 a = aa;
  72 b = bb;
  73 c = cc;
  74 x = xx;
  75 an = a;
  76 bn = b;
  77 cn = c;
  78 a0 = 1.0;
  79 sum = 1.0;
  80 n = 1.0;
  81 t = 1.0;
  82 max = 0.0;
  83
  84 do
  85         {
  86         if( an == 0 )
  87                 goto done;
  88         if( bn == 0 )
  89                 goto error;
  90         if( cn == 0 )
  91                 goto error;
  92         if( (a0 > 1.0e34) || (n > 200) )
  93                 goto error;
  94         a0 *= (an * x) / (bn * cn * n);
  95         sum += a0;
  96         an += 1.0;
  97         bn += 1.0;
  98         cn += 1.0;
  99         n += 1.0;
 100         z = fabsf( a0 );
 101         if( z > max )
 102                 max = z;
 103         if( sum != 0 )
 104                 t = fabsf( a0 / sum );
 105         else
 106                 t = z;
 107         }
 108 while( t > MACHEPF );
 109
 110 done:
 111
 112 *err = fabsf( MACHEPF*max /sum );
 113
 114 #if DEBUG
 115         printf(" onef2f cancellation error %.5E\n", *err );
 116 #endif
 117
 118 goto xit;
 119
 120 error:
 121 #if DEBUG
 122 printf("onef2f does not converge\n");
 123 #endif
 124 *err = MAXNUMF;
 125
 126 xit:
 127
 128 #if DEBUG
 129 printf("onef2( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
 130 #endif
 131 return(sum);
 132 }
 133
 134
 135
 136 #ifdef ANSIC
 137 float threef0f( float aa, float bb, float cc, float xx, float *err )
 138 #else
 139 float threef0f( aa, bb, cc, xx, err )
 140 double aa, bb, cc, xx;
 141 float *err;
 142 #endif
 143 {
 144 float a, b, c, x, n, a0, sum, t, conv, conv1;
 145 float an, bn, cn, max, z;
 146
 147 a = aa;
 148 b = bb;
 149 c = cc;
 150 x = xx;
 151 an = a;
 152 bn = b;
 153 cn = c;
 154 a0 = 1.0;
 155 sum = 1.0;
 156 n = 1.0;
 157 t = 1.0;
 158 max = 0.0;
 159 conv = 1.0e38;
 160 conv1 = conv;
 161
 162 do
 163         {
 164         if( an == 0.0 )
 165                 goto done;
 166         if( bn == 0.0 )
 167                 goto done;
 168         if( cn == 0.0 )
 169                 goto done;
 170         if( (a0 > 1.0e34) || (n > 200) )
 171                 goto error;
 172         a0 *= (an * bn * cn * x) / n;
 173         an += 1.0;
 174         bn += 1.0;
 175         cn += 1.0;
 176         n += 1.0;
 177         z = fabsf( a0 );
 178         if( z > max )
 179                 max = z;
 180         if( z >= conv )
 181                 {
 182                 if( (z < max) && (z > conv1) )
 183                         goto done;
 184                 }
 185         conv1 = conv;
 186         conv = z;
 187         sum += a0;
 188         if( sum != 0 )
 189                 t = fabsf( a0 / sum );
 190         else
 191                 t = z;
 192         }
 193 while( t > MACHEPF );
 194
 195 done:
 196
 197 t = fabsf( MACHEPF*max/sum );
 198 #if DEBUG
 199         printf(" threef0f cancellation error %.5E\n", t );
 200 #endif
 201
 202 max = fabsf( conv/sum );
 203 if( max > t )
 204         t = max;
 205 #if DEBUG
 206         printf(" threef0f convergence %.5E\n", max );
 207 #endif
 208
 209 goto xit;
 210
 211 error:
 212 #if DEBUG
 213 printf("threef0f does not converge\n");
 214 #endif
 215 t = MAXNUMF;
 216
 217 xit:
 218
 219 #if DEBUG
 220 printf("threef0f( %.2E %.2E %.2E %.5E ) =  %.3E  %.6E\n", a, b, c, x, n, sum);
 221 #endif
 222
 223 *err = t;
 224 return(sum);
 225 }
 226
 227
 228
 229
 230 #ifdef ANSIC
 231 float struvef( float vv, float xx )
 232 #else
 233 float struvef( vv, xx )
 234 double vv, xx;
 235 #endif
 236 {
 237 float v, x, y, ya, f, g, h, t;
 238 float onef2err, threef0err;
 239
 240 v = vv;
 241 x = xx;
 242 f = floorf(v);
 243 if( (v < 0) && ( v-f == 0.5 ) )
 244         {
 245         y = jvf( -v, x );
 246         f = 1.0 - f;
 247         g =  2.0 * floorf(0.5*f);
 248         if( g != f )
 249                 y = -y;
 250         return(y);
 251         }
 252 t = 0.25*x*x;
 253 f = fabsf(x);
 254 g = 1.5 * fabsf(v);
 255 if( (f > 30.0) && (f > g) )
 256         {
 257         onef2err = MAXNUMF;
 258         y = 0.0;
 259         }
 260 else
 261         {
 262         y = onef2f( 1.0, 1.5, 1.5+v, -t, &onef2err );
 263         }
 264
 265 if( (f < 18.0) || (x < 0.0) )
 266         {
 267         threef0err = MAXNUMF;
 268         ya = 0.0;
 269         }
 270 else
 271         {
 272         ya = threef0f( 1.0, 0.5, 0.5-v, -1.0/t, &threef0err );
 273         }
 274
 275 f = sqrtf( PIF );
 276 h = powf( 0.5*x, v-1.0 );
 277
 278 if( onef2err <= threef0err )
 279         {
 280         g = gammaf( v + 1.5 );
 281         y = y * h * t / ( 0.5 * f * g );
 282         return(y);
 283         }
 284 else
 285         {
 286         g = gammaf( v + 0.5 );
 287         ya = ya * h / ( f * g );
 288         ya = ya + yvf( v, x );
 289         return(ya);
 290         }
 291 }
 292
 293
 294
 295
 296 /* Bessel function of noninteger order
 297  */
 298
 299 #ifdef ANSIC
 300 float yvf( float vv, float xx )
 301 #else
 302 float yvf( vv, xx )
 303 double vv, xx;
 304 #endif
 305 {
 306 float v, x,  y, t;
 307 int n;
 308
 309 v = vv;
 310 x = xx;
 311 y = floorf( v );
 312 if( y == v )
 313         {
 314         n = v;
 315         y = ynf( n, x );
 316         return( y );
 317         }
 318 t = PIF * v;
 319 y = (cosf(t) * jvf( v, x ) - jvf( -v, x ))/sinf(t);
 320 return( y );
 321 }
 322
 323 /* Crossover points between ascending series and asymptotic series
 324  * for Struve function
 325  *
 326  *       v       x
 327  *
 328  *       0      19.2
 329  *       1      18.95
 330  *       2      19.15
 331  *       3      19.3
 332  *       5      19.7
 333  *      10      21.35
 334  *      20      26.35
 335  *      30      32.31
 336  *      40      40.0
 337  */