lib/libm/hyp2f1.c

   1 /*                                                      hyp2f1.c
   2  *
   3  *      Gauss hypergeometric function   F
   4  *                                     2 1
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double a, b, c, x, y, hyp2f1();
  10  *
  11  * y = hyp2f1( a, b, c, x );
  12  *
  13  *
  14  * DESCRIPTION:
  15  *
  16  *
  17  *  hyp2f1( a, b, c, x )  =   F ( a, b; c; x )
  18  *                           2 1
  19  *
  20  *           inf.
  21  *            -   a(a+1)...(a+k) b(b+1)...(b+k)   k+1
  22  *   =  1 +   >   -----------------------------  x   .
  23  *            -         c(c+1)...(c+k) (k+1)!
  24  *          k = 0
  25  *
  26  *  Cases addressed are
  27  *      Tests and escapes for negative integer a, b, or c
  28  *      Linear transformation if c - a or c - b negative integer
  29  *      Special case c = a or c = b
  30  *      Linear transformation for  x near +1
  31  *      Transformation for x < -0.5
  32  *      Psi function expansion if x > 0.5 and c - a - b integer
  33  *      Conditionally, a recurrence on c to make c-a-b > 0
  34  *
  35  * |x| > 1 is rejected.
  36  *
  37  * The parameters a, b, c are considered to be integer
  38  * valued if they are within 1.0e-14 of the nearest integer
  39  * (1.0e-13 for IEEE arithmetic).
  40  *
  41  * ACCURACY:
  42  *
  43  *
  44  *               Relative error (-1 < x < 1):
  45  * arithmetic   domain     # trials      peak         rms
  46  *    IEEE      -1,7        230000      1.2e-11     5.2e-14
  47  *
  48  * Several special cases also tested with a, b, c in
  49  * the range -7 to 7.
  50  *
  51  * ERROR MESSAGES:
  52  *
  53  * A "partial loss of precision" message is printed if
  54  * the internally estimated relative error exceeds 1^-12.
  55  * A "singularity" message is printed on overflow or
  56  * in cases not addressed (such as x < -1).
  57  */
  58 \f
  59 /*                                                      hyp2f1  */
  60
  61
  62 /*
  63 Cephes Math Library Release 2.8:  June, 2000
  64 Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
  65 */
  66
  67
  68 #include "mconf.h"
  69
  70 #ifdef DEC
  71 #define EPS 1.0e-14
  72 #define EPS2 1.0e-11
  73 #endif
  74
  75 #ifdef IBMPC
  76 #define EPS 1.0e-13
  77 #define EPS2 1.0e-10
  78 #endif
  79
  80 #ifdef MIEEE
  81 #define EPS 1.0e-13
  82 #define EPS2 1.0e-10
  83 #endif
  84
  85 #ifdef UNK
  86 #define EPS 1.0e-13
  87 #define EPS2 1.0e-10
  88 #endif
  89
  90 #define ETHRESH 1.0e-12
  91
  92 #ifdef ANSIPROT
  93 extern double fabs ( double );
  94 extern double pow ( double, double );
  95 extern double round ( double );
  96 extern double gamma ( double );
  97 extern double log ( double );
  98 extern double exp ( double );
  99 extern double psi ( double );
 100 static double hyt2f1(double, double, double, double, double *);
 101 static double hys2f1(double, double, double, double, double *);
 102 double hyp2f1(double, double, double, double);
 103 #else
 104 double fabs(), pow(), round(), gamma(), log(), exp(), psi();
 105 static double hyt2f1();
 106 static double hys2f1();
 107 double hyp2f1();
 108 #endif
 109 extern double MAXNUM, MACHEP;
 110
 111 double hyp2f1( a, b, c, x )
 112 double a, b, c, x;
 113 {
 114 double d, d1, d2, e;
 115 double p, q, r, s, y, ax;
 116 double ia, ib, ic, id, err;
 117 int flag, i, aid;
 118
 119 err = 0.0;
 120 ax = fabs(x);
 121 s = 1.0 - x;
 122 flag = 0;
 123 ia = round(a); /* nearest integer to a */
 124 ib = round(b);
 125
 126 if( a <= 0 )
 127         {
 128         if( fabs(a-ia) < EPS )          /* a is a negative integer */
 129                 flag |= 1;
 130         }
 131
 132 if( b <= 0 )
 133         {
 134         if( fabs(b-ib) < EPS )          /* b is a negative integer */
 135                 flag |= 2;
 136         }
 137
 138 if( ax < 1.0 )
 139         {
 140         if( fabs(b-c) < EPS )           /* b = c */
 141                 {
 142                 y = pow( s, -a );       /* s to the -a power */
 143                 goto hypdon;
 144                 }
 145         if( fabs(a-c) < EPS )           /* a = c */
 146                 {
 147                 y = pow( s, -b );       /* s to the -b power */
 148                 goto hypdon;
 149                 }
 150         }
 151
 152
 153
 154 if( c <= 0.0 )
 155         {
 156         ic = round(c);  /* nearest integer to c */
 157         if( fabs(c-ic) < EPS )          /* c is a negative integer */
 158                 {
 159                 /* check if termination before explosion */
 160                 if( (flag & 1) && (ia > ic) )
 161                         goto hypok;
 162                 if( (flag & 2) && (ib > ic) )
 163                         goto hypok;
 164                 goto hypdiv;
 165                 }
 166         }
 167
 168 if( flag )                      /* function is a polynomial */
 169         goto hypok;
 170
 171 if( ax > 1.0 )                  /* series diverges      */
 172         goto hypdiv;
 173
 174 p = c - a;
 175 ia = round(p); /* nearest integer to c-a */
 176 if( (ia <= 0.0) && (fabs(p-ia) < EPS) ) /* negative int c - a */
 177         flag |= 4;
 178
 179 r = c - b;
 180 ib = round(r); /* nearest integer to c-b */
 181 if( (ib <= 0.0) && (fabs(r-ib) < EPS) ) /* negative int c - b */
 182         flag |= 8;
 183
 184 d = c - a - b;
 185 id = round(d); /* nearest integer to d */
 186 q = fabs(d-id);
 187
 188 /* Thanks to Christian Burger <BURGER@DMRHRZ11.HRZ.Uni-Marburg.DE>
 189  * for reporting a bug here.  */
 190 if( fabs(ax-1.0) < EPS )                        /* |x| == 1.0   */
 191         {
 192         if( x > 0.0 )
 193                 {
 194                 if( flag & 12 ) /* negative int c-a or c-b */
 195                         {
 196                         if( d >= 0.0 )
 197                                 goto hypf;
 198                         else
 199                                 goto hypdiv;
 200                         }
 201                 if( d <= 0.0 )
 202                         goto hypdiv;
 203                 y = gamma(c)*gamma(d)/(gamma(p)*gamma(r));
 204                 goto hypdon;
 205                 }
 206
 207         if( d <= -1.0 )
 208                 goto hypdiv;
 209
 210         }
 211
 212 /* Conditionally make d > 0 by recurrence on c
 213  * AMS55 #15.2.27
 214  */
 215 if( d < 0.0 )
 216         {
 217 /* Try the power series first */
 218         y = hyt2f1( a, b, c, x, &err );
 219         if( err < ETHRESH )
 220                 goto hypdon;
 221 /* Apply the recurrence if power series fails */
 222         err = 0.0;
 223         aid = 2 - id;
 224         e = c + aid;
 225         d2 = hyp2f1(a,b,e,x);
 226         d1 = hyp2f1(a,b,e+1.0,x);
 227         q = a + b + 1.0;
 228         for( i=0; i<aid; i++ )
 229                 {
 230                 r = e - 1.0;
 231                 y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s);
 232                 e = r;
 233                 d1 = d2;
 234                 d2 = y;
 235                 }
 236         goto hypdon;
 237         }
 238
 239
 240 if( flag & 12 )
 241         goto hypf; /* negative integer c-a or c-b */
 242
 243 hypok:
 244 y = hyt2f1( a, b, c, x, &err );
 245
 246
 247 hypdon:
 248 if( err > ETHRESH )
 249         {
 250         mtherr( "hyp2f1", PLOSS );
 251 /*      printf( "Estimated err = %.2e\n", err ); */
 252         }
 253 return(y);
 254
 255 /* The transformation for c-a or c-b negative integer
 256  * AMS55 #15.3.3
 257  */
 258 hypf:
 259 y = pow( s, d ) * hys2f1( c-a, c-b, c, x, &err );
 260 goto hypdon;
 261
 262 /* The alarm exit */
 263 hypdiv:
 264 mtherr( "hyp2f1", OVERFLOW );
 265 return( MAXNUM );
 266 }
 267
 268
 269
 270
 271
 272
 273 /* Apply transformations for |x| near 1
 274  * then call the power series
 275  */
 276 static double hyt2f1( a, b, c, x, loss )
 277 double a, b, c, x;
 278 double *loss;
 279 {
 280 double p, q, r, s, t, y, d, err, err1;
 281 double ax, id, d1, d2, e, y1;
 282 int i, aid;
 283
 284 err = 0.0;
 285 s = 1.0 - x;
 286 if( x < -0.5 )
 287         {
 288         if( b > a )
 289                 y = pow( s, -a ) * hys2f1( a, c-b, c, -x/s, &err );
 290
 291         else
 292                 y = pow( s, -b ) * hys2f1( c-a, b, c, -x/s, &err );
 293
 294         goto done;
 295         }
 296
 297 d = c - a - b;
 298 id = round(d);  /* nearest integer to d */
 299
 300 if( x > 0.9 )
 301 {
 302 if( fabs(d-id) > EPS ) /* test for integer c-a-b */
 303         {
 304 /* Try the power series first */
 305         y = hys2f1( a, b, c, x, &err );
 306         if( err < ETHRESH )
 307                 goto done;
 308 /* If power series fails, then apply AMS55 #15.3.6 */
 309         q = hys2f1( a, b, 1.0-d, s, &err );
 310         q *= gamma(d) /(gamma(c-a) * gamma(c-b));
 311         r = pow(s,d) * hys2f1( c-a, c-b, d+1.0, s, &err1 );
 312         r *= gamma(-d)/(gamma(a) * gamma(b));
 313         y = q + r;
 314
 315         q = fabs(q); /* estimate cancellation error */
 316         r = fabs(r);
 317         if( q > r )
 318                 r = q;
 319         err += err1 + (MACHEP*r)/y;
 320
 321         y *= gamma(c);
 322         goto done;
 323         }
 324 else
 325         {
 326 /* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */
 327         if( id >= 0.0 )
 328                 {
 329                 e = d;
 330                 d1 = d;
 331                 d2 = 0.0;
 332                 aid = id;
 333                 }
 334         else
 335                 {
 336                 e = -d;
 337                 d1 = 0.0;
 338                 d2 = d;
 339                 aid = -id;
 340                 }
 341
 342         ax = log(s);
 343
 344         /* sum for t = 0 */
 345         y = psi(1.0) + psi(1.0+e) - psi(a+d1) - psi(b+d1) - ax;
 346         y /= gamma(e+1.0);
 347
 348         p = (a+d1) * (b+d1) * s / gamma(e+2.0); /* Poch for t=1 */
 349         t = 1.0;
 350         do
 351                 {
 352                 r = psi(1.0+t) + psi(1.0+t+e) - psi(a+t+d1)
 353                         - psi(b+t+d1) - ax;
 354                 q = p * r;
 355                 y += q;
 356                 p *= s * (a+t+d1) / (t+1.0);
 357                 p *= (b+t+d1) / (t+1.0+e);
 358                 t += 1.0;
 359                 }
 360         while( fabs(q/y) > EPS );
 361
 362
 363         if( id == 0.0 )
 364                 {
 365                 y *= gamma(c)/(gamma(a)*gamma(b));
 366                 goto psidon;
 367                 }
 368
 369         y1 = 1.0;
 370
 371         if( aid == 1 )
 372                 goto nosum;
 373
 374         t = 0.0;
 375         p = 1.0;
 376         for( i=1; i<aid; i++ )
 377                 {
 378                 r = 1.0-e+t;
 379                 p *= s * (a+t+d2) * (b+t+d2) / r;
 380                 t += 1.0;
 381                 p /= t;
 382                 y1 += p;
 383                 }
 384 nosum:
 385         p = gamma(c);
 386         y1 *= gamma(e) * p / (gamma(a+d1) * gamma(b+d1));
 387
 388         y *= p / (gamma(a+d2) * gamma(b+d2));
 389         if( (aid & 1) != 0 )
 390                 y = -y;
 391
 392         q = pow( s, id );       /* s to the id power */
 393         if( id > 0.0 )
 394                 y *= q;
 395         else
 396                 y1 *= q;
 397
 398         y += y1;
 399 psidon:
 400         goto done;
 401         }
 402
 403 }
 404
 405 /* Use defining power series if no special cases */
 406 y = hys2f1( a, b, c, x, &err );
 407
 408 done:
 409 *loss = err;
 410 return(y);
 411 }
 412
 413
 414
 415
 416
 417 /* Defining power series expansion of Gauss hypergeometric function */
 418
 419 static double hys2f1( a, b, c, x, loss )
 420 double a, b, c, x;
 421 double *loss; /* estimates loss of significance */
 422 {
 423 double f, g, h, k, m, s, u, umax;
 424 int i;
 425
 426 i = 0;
 427 umax = 0.0;
 428 f = a;
 429 g = b;
 430 h = c;
 431 s = 1.0;
 432 u = 1.0;
 433 k = 0.0;
 434 do
 435         {
 436         if( fabs(h) < EPS )
 437                 {
 438                 *loss = 1.0;
 439                 return( MAXNUM );
 440                 }
 441         m = k + 1.0;
 442         u = u * ((f+k) * (g+k) * x / ((h+k) * m));
 443         s += u;
 444         k = fabs(u);  /* remember largest term summed */
 445         if( k > umax )
 446                 umax = k;
 447         k = m;
 448         if( ++i > 10000 ) /* should never happen */
 449                 {
 450                 *loss = 1.0;
 451                 return(s);
 452                 }
 453         }
 454 while( fabs(u/s) > MACHEP );
 455
 456 /* return estimated relative error */
 457 *loss = (MACHEP*umax)/fabs(s) + (MACHEP*i);
 458
 459 return(s);
 460 }