lib/libm/hyperg.c

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