lib/libm/hypergf.c

   1 /*                                                      hypergf.c
   2  *
   3  *      Confluent hypergeometric function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float a, b, x, y, hypergf();
  10  *
  11  * y = hypergf( a, b, x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Computes the confluent hypergeometric function
  18  *
  19  *                          1           2
  20  *                       a x    a(a+1) x
  21  *   F ( a,b;x )  =  1 + ---- + --------- + ...
  22  *  1 1                  b 1!   b(b+1) 2!
  23  *
  24  * Many higher transcendental functions are special cases of
  25  * this power series.
  26  *
  27  * As is evident from the formula, b must not be a negative
  28  * integer or zero unless a is an integer with 0 >= a > b.
  29  *
  30  * The routine attempts both a direct summation of the series
  31  * and an asymptotic expansion.  In each case error due to
  32  * roundoff, cancellation, and nonconvergence is estimated.
  33  * The result with smaller estimated error is returned.
  34  *
  35  *
  36  *
  37  * ACCURACY:
  38  *
  39  * Tested at random points (a, b, x), all three variables
  40  * ranging from 0 to 30.
  41  *                      Relative error:
  42  * arithmetic   domain     # trials      peak         rms
  43  *    IEEE      0,5         10000       6.6e-7      1.3e-7
  44  *    IEEE      0,30        30000       1.1e-5      6.5e-7
  45  *
  46  * Larger errors can be observed when b is near a negative
  47  * integer or zero.  Certain combinations of arguments yield
  48  * serious cancellation error in the power series summation
  49  * and also are not in the region of near convergence of the
  50  * asymptotic series.  An error message is printed if the
  51  * self-estimated relative error is greater than 1.0e-3.
  52  *
  53  */
  54 \f
  55 /*                                                      hyperg.c */
  56
  57
  58 /*
  59 Cephes Math Library Release 2.1:  November, 1988
  60 Copyright 1984, 1987, 1988 by Stephen L. Moshier
  61 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  62 */
  63
  64 #include "mconf.h"
  65
  66 extern float MAXNUMF, MACHEPF;
  67
  68 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  69
  70 #ifdef ANSIC
  71 float expf(float);
  72 float hyp2f0f(float, float, float, int, float *);
  73 static float hy1f1af(float, float, float, float *);
  74 static float hy1f1pf(float, float, float, float *);
  75 float logf(float), gammaf(float), lgamf(float);
  76 #else
  77 float expf(), hyp2f0f();
  78 float logf(), gammaf(), lgamf();
  79 static float hy1f1pf(), hy1f1af();
  80 #endif
  81
  82 #ifdef ANSIC
  83 float hypergf( float aa, float bb, float xx )
  84 #else
  85 float hypergf( aa, bb, xx)
  86 double aa, bb, xx;
  87 #endif
  88 {
  89 float a, b, x, asum, psum, acanc, pcanc, temp;
  90
  91
  92 a = aa;
  93 b = bb;
  94 x = xx;
  95 /* See if a Kummer transformation will help */
  96 temp = b - a;
  97 if( fabsf(temp) < 0.001 * fabsf(a) )
  98         return( expf(x) * hypergf( temp, b, -x )  );
  99
 100 psum = hy1f1pf( a, b, x, &pcanc );
 101 if( pcanc < 1.0e-6 )
 102         goto done;
 103
 104
 105 /* try asymptotic series */
 106
 107 asum = hy1f1af( a, b, x, &acanc );
 108
 109
 110 /* Pick the result with less estimated error */
 111
 112 if( acanc < pcanc )
 113         {
 114         pcanc = acanc;
 115         psum = asum;
 116         }
 117
 118 done:
 119 if( pcanc > 1.0e-3 )
 120         mtherr( "hyperg", PLOSS );
 121
 122 return( psum );
 123 }
 124
 125
 126
 127
 128 /* Power series summation for confluent hypergeometric function         */
 129
 130
 131 #ifdef ANSIC
 132 static float hy1f1pf( float aa, float bb, float xx, float *err )
 133 #else
 134 static float hy1f1pf( aa, bb, xx, err )
 135 double aa, bb, xx;
 136 float *err;
 137 #endif
 138 {
 139 float a, b, x, n, a0, sum, t, u, temp;
 140 float an, bn, maxt, pcanc;
 141
 142 a = aa;
 143 b = bb;
 144 x = xx;
 145 /* set up for power series summation */
 146 an = a;
 147 bn = b;
 148 a0 = 1.0;
 149 sum = 1.0;
 150 n = 1.0;
 151 t = 1.0;
 152 maxt = 0.0;
 153
 154
 155 while( t > MACHEPF )
 156         {
 157         if( bn == 0 )                   /* check bn first since if both */
 158                 {
 159                 mtherr( "hypergf", SING );
 160                 return( MAXNUMF );      /* an and bn are zero it is     */
 161                 }
 162         if( an == 0 )                   /* a singularity                */
 163                 return( sum );
 164         if( n > 200 )
 165                 goto pdone;
 166         u = x * ( an / (bn * n) );
 167
 168         /* check for blowup */
 169         temp = fabsf(u);
 170         if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) )
 171                 {
 172                 pcanc = 1.0;    /* estimate 100% error */
 173                 goto blowup;
 174                 }
 175
 176         a0 *= u;
 177         sum += a0;
 178         t = fabsf(a0);
 179         if( t > maxt )
 180                 maxt = t;
 181 /*
 182         if( (maxt/fabsf(sum)) > 1.0e17 )
 183                 {
 184                 pcanc = 1.0;
 185                 goto blowup;
 186                 }
 187 */
 188         an += 1.0;
 189         bn += 1.0;
 190         n += 1.0;
 191         }
 192
 193 pdone:
 194
 195 /* estimate error due to roundoff and cancellation */
 196 if( sum != 0.0 )
 197         maxt /= fabsf(sum);
 198 maxt *= MACHEPF;        /* this way avoids multiply overflow */
 199 pcanc = fabsf( MACHEPF * n  +  maxt );
 200
 201 blowup:
 202
 203 *err = pcanc;
 204
 205 return( sum );
 206 }
 207
 208
 209 /*                                                      hy1f1a()        */
 210 /* asymptotic formula for hypergeometric function:
 211  *
 212  *        (    -a
 213  *  --    ( |z|
 214  * |  (b) ( -------- 2f0( a, 1+a-b, -1/x )
 215  *        (  --
 216  *        ( |  (b-a)
 217  *
 218  *
 219  *                                x    a-b                     )
 220  *                               e  |x|                        )
 221  *                             + -------- 2f0( b-a, 1-a, 1/x ) )
 222  *                                --                           )
 223  *                               |  (a)                        )
 224  */
 225
 226 #ifdef ANSIC
 227 static float hy1f1af( float aa, float bb, float xx, float *err )
 228 #else
 229 static float hy1f1af( aa, bb, xx, err )
 230 double aa, bb, xx;
 231 float *err;
 232 #endif
 233 {
 234 float a, b, x, h1, h2, t, u, temp, acanc, asum, err1, err2;
 235
 236 a = aa;
 237 b = bb;
 238 x = xx;
 239 if( x == 0 )
 240         {
 241         acanc = 1.0;
 242         asum = MAXNUMF;
 243         goto adone;
 244         }
 245 temp = logf( fabsf(x) );
 246 t = x + temp * (a-b);
 247 u = -temp * a;
 248
 249 if( b > 0 )
 250         {
 251         temp = lgamf(b);
 252         t += temp;
 253         u += temp;
 254         }
 255
 256 h1 = hyp2f0f( a, a-b+1, -1.0/x, 1, &err1 );
 257
 258 temp = expf(u) / gammaf(b-a);
 259 h1 *= temp;
 260 err1 *= temp;
 261
 262 h2 = hyp2f0f( b-a, 1.0-a, 1.0/x, 2, &err2 );
 263
 264 if( a < 0 )
 265         temp = expf(t) / gammaf(a);
 266 else
 267         temp = expf( t - lgamf(a) );
 268
 269 h2 *= temp;
 270 err2 *= temp;
 271
 272 if( x < 0.0 )
 273         asum = h1;
 274 else
 275         asum = h2;
 276
 277 acanc = fabsf(err1) + fabsf(err2);
 278
 279
 280 if( b < 0 )
 281         {
 282         temp = gammaf(b);
 283         asum *= temp;
 284         acanc *= fabsf(temp);
 285         }
 286
 287
 288 if( asum != 0.0 )
 289         acanc /= fabsf(asum);
 290
 291 acanc *= 30.0;  /* fudge factor, since error of asymptotic formula
 292                  * often seems this much larger than advertised */
 293
 294 adone:
 295
 296
 297 *err = acanc;
 298 return( asum );
 299 }
 300 \f
 301 /*                                                      hyp2f0()        */
 302
 303 #ifdef ANSIC
 304 float hyp2f0f(float aa, float bb, float xx, int type, float *err)
 305 #else
 306 float hyp2f0f( aa, bb, xx, type, err )
 307 double aa, bb, xx;
 308 int type;       /* determines what converging factor to use */
 309 float *err;
 310 #endif
 311 {
 312 float a, b, x, a0, alast, t, tlast, maxt;
 313 float n, an, bn, u, sum, temp;
 314
 315 a = aa;
 316 b = bb;
 317 x = xx;
 318 an = a;
 319 bn = b;
 320 a0 = 1.0;
 321 alast = 1.0;
 322 sum = 0.0;
 323 n = 1.0;
 324 t = 1.0;
 325 tlast = 1.0e9;
 326 maxt = 0.0;
 327
 328 do
 329         {
 330         if( an == 0 )
 331                 goto pdone;
 332         if( bn == 0 )
 333                 goto pdone;
 334
 335         u = an * (bn * x / n);
 336
 337         /* check for blowup */
 338         temp = fabsf(u);
 339         if( (temp > 1.0 ) && (maxt > (MAXNUMF/temp)) )
 340                 goto error;
 341
 342         a0 *= u;
 343         t = fabsf(a0);
 344
 345         /* terminating condition for asymptotic series */
 346         if( t > tlast )
 347                 goto ndone;
 348
 349         tlast = t;
 350         sum += alast;   /* the sum is one term behind */
 351         alast = a0;
 352
 353         if( n > 200 )
 354                 goto ndone;
 355
 356         an += 1.0;
 357         bn += 1.0;
 358         n += 1.0;
 359         if( t > maxt )
 360                 maxt = t;
 361         }
 362 while( t > MACHEPF );
 363
 364
 365 pdone:  /* series converged! */
 366
 367 /* estimate error due to roundoff and cancellation */
 368 *err = fabsf(  MACHEPF * (n + maxt)  );
 369
 370 alast = a0;
 371 goto done;
 372
 373 ndone:  /* series did not converge */
 374
 375 /* The following "Converging factors" are supposed to improve accuracy,
 376  * but do not actually seem to accomplish very much. */
 377
 378 n -= 1.0;
 379 x = 1.0/x;
 380
 381 switch( type )  /* "type" given as subroutine argument */
 382 {
 383 case 1:
 384         alast *= ( 0.5 + (0.125 + 0.25*b - 0.5*a + 0.25*x - 0.25*n)/x );
 385         break;
 386
 387 case 2:
 388         alast *= 2.0/3.0 - b + 2.0*a + x - n;
 389         break;
 390
 391 default:
 392         ;
 393 }
 394
 395 /* estimate error due to roundoff, cancellation, and nonconvergence */
 396 *err = MACHEPF * (n + maxt)  +  fabsf( a0 );
 397
 398
 399 done:
 400 sum += alast;
 401 return( sum );
 402
 403 /* series blew up: */
 404 error:
 405 *err = MAXNUMF;
 406 mtherr( "hypergf", TLOSS );
 407 return( sum );
 408 }