lib/libm/hyp2f1f.c

   1 /*                                                      hyp2f1f.c
   2  *
   3  *      Gauss hypergeometric function   F
   4  *                                     2 1
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float a, b, c, x, y, hyp2f1f();
  10  *
  11  * y = hyp2f1f( 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-6 of the nearest integer.
  39  *
  40  * ACCURACY:
  41  *
  42  *                      Relative error (-1 < x < 1):
  43  * arithmetic   domain     # trials      peak         rms
  44  *    IEEE      0,3         30000       5.8e-4      4.3e-6
  45  */
  46 \f
  47 /*                                                      hyp2f1  */
  48
  49
  50 /*
  51 Cephes Math Library Release 2.2:  November, 1992
  52 Copyright 1984, 1987, 1992 by Stephen L. Moshier
  53 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54 */
  55
  56
  57 #include "mconf.h"
  58
  59 #define EPS 1.0e-5
  60 #define EPS2 1.0e-5
  61 #define ETHRESH 1.0e-5
  62
  63 extern float MAXNUMF, MACHEPF;
  64
  65 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  66
  67 #ifdef ANSIC
  68 float powf(float, float);
  69 static float hys2f1f(float, float, float, float, float *);
  70 static float hyt2f1f(float, float, float, float, float *);
  71 float gammaf(float), logf(float), expf(float), psif(float);
  72 float floorf(float);
  73 #else
  74 float powf(), gammaf(), logf(), expf(), psif();
  75 float floorf();
  76 static float hyt2f1f(), hys2f1f();
  77 #endif
  78
  79 #define roundf(x) (floorf((x)+(float )0.5))
  80
  81
  82
  83
  84 #ifdef ANSIC
  85 float hyp2f1f( float aa, float bb, float cc, float xx )
  86 #else
  87 float hyp2f1f( aa, bb, cc, xx )
  88 double aa, bb, cc, xx;
  89 #endif
  90 {
  91 float a, b, c, x;
  92 float d, d1, d2, e;
  93 float p, q, r, s, y, ax;
  94 float ia, ib, ic, id, err;
  95 int flag, i, aid;
  96
  97 a = aa;
  98 b = bb;
  99 c = cc;
 100 x = xx;
 101 err = 0.0;
 102 ax = fabsf(x);
 103 s = 1.0 - x;
 104 flag = 0;
 105 ia = roundf(a); /* nearest integer to a */
 106 ib = roundf(b);
 107
 108 if( a <= 0 )
 109         {
 110         if( fabsf(a-ia) < EPS )         /* a is a negative integer */
 111                 flag |= 1;
 112         }
 113
 114 if( b <= 0 )
 115         {
 116         if( fabsf(b-ib) < EPS )         /* b is a negative integer */
 117                 flag |= 2;
 118         }
 119
 120 if( ax < 1.0 )
 121         {
 122         if( fabsf(b-c) < EPS )          /* b = c */
 123                 {
 124                 y = powf( s, -a );      /* s to the -a power */
 125                 goto hypdon;
 126                 }
 127         if( fabsf(a-c) < EPS )          /* a = c */
 128                 {
 129                 y = powf( s, -b );      /* s to the -b power */
 130                 goto hypdon;
 131                 }
 132         }
 133
 134
 135
 136 if( c <= 0.0 )
 137         {
 138         ic = roundf(c);         /* nearest integer to c */
 139         if( fabsf(c-ic) < EPS )         /* c is a negative integer */
 140                 {
 141                 /* check if termination before explosion */
 142                 if( (flag & 1) && (ia > ic) )
 143                         goto hypok;
 144                 if( (flag & 2) && (ib > ic) )
 145                         goto hypok;
 146                 goto hypdiv;
 147                 }
 148         }
 149
 150 if( flag )                      /* function is a polynomial */
 151         goto hypok;
 152
 153 if( ax > 1.0 )                  /* series diverges      */
 154         goto hypdiv;
 155
 156 p = c - a;
 157 ia = roundf(p);
 158 if( (ia <= 0.0) && (fabsf(p-ia) < EPS) )        /* negative int c - a */
 159         flag |= 4;
 160
 161 r = c - b;
 162 ib = roundf(r); /* nearest integer to r */
 163 if( (ib <= 0.0) && (fabsf(r-ib) < EPS) )        /* negative int c - b */
 164         flag |= 8;
 165
 166 d = c - a - b;
 167 id = roundf(d); /* nearest integer to d */
 168 q = fabsf(d-id);
 169
 170 if( fabsf(ax-1.0) < EPS )                       /* |x| == 1.0   */
 171         {
 172         if( x > 0.0 )
 173                 {
 174                 if( flag & 12 ) /* negative int c-a or c-b */
 175                         {
 176                         if( d >= 0.0 )
 177                                 goto hypf;
 178                         else
 179                                 goto hypdiv;
 180                         }
 181                 if( d <= 0.0 )
 182                         goto hypdiv;
 183                 y = gammaf(c)*gammaf(d)/(gammaf(p)*gammaf(r));
 184                 goto hypdon;
 185                 }
 186
 187         if( d <= -1.0 )
 188                 goto hypdiv;
 189         }
 190
 191 /* Conditionally make d > 0 by recurrence on c
 192  * AMS55 #15.2.27
 193  */
 194 if( d < 0.0 )
 195         {
 196 /* Try the power series first */
 197         y = hyt2f1f( a, b, c, x, &err );
 198         if( err < ETHRESH )
 199                 goto hypdon;
 200 /* Apply the recurrence if power series fails */
 201         err = 0.0;
 202         aid = 2 - id;
 203         e = c + aid;
 204         d2 = hyp2f1f(a,b,e,x);
 205         d1 = hyp2f1f(a,b,e+1.0,x);
 206         q = a + b + 1.0;
 207         for( i=0; i<aid; i++ )
 208                 {
 209                 r = e - 1.0;
 210                 y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s);
 211                 e = r;
 212                 d1 = d2;
 213                 d2 = y;
 214                 }
 215         goto hypdon;
 216         }
 217
 218
 219 if( flag & 12 )
 220         goto hypf; /* negative integer c-a or c-b */
 221
 222 hypok:
 223 y = hyt2f1f( a, b, c, x, &err );
 224
 225 hypdon:
 226 if( err > ETHRESH )
 227         {
 228         mtherr( "hyp2f1", PLOSS );
 229 /*      printf( "Estimated err = %.2e\n", err );*/
 230         }
 231 return(y);
 232
 233 /* The transformation for c-a or c-b negative integer
 234  * AMS55 #15.3.3
 235  */
 236 hypf:
 237 y = powf( s, d ) * hys2f1f( c-a, c-b, c, x, &err );
 238 goto hypdon;
 239
 240 /* The alarm exit */
 241 hypdiv:
 242 mtherr( "hyp2f1f", OVERFLOW );
 243 return( MAXNUMF );
 244 }
 245
 246
 247
 248
 249 /* Apply transformations for |x| near 1
 250  * then call the power series
 251  */
 252 #ifdef ANSIC
 253 static float hyt2f1f( float aa, float bb, float cc, float xx, float *loss )
 254 #else
 255 static float hyt2f1f( aa, bb, cc, xx, loss )
 256 double aa, bb, cc, xx;
 257 float *loss;
 258 #endif
 259 {
 260 float a, b, c, x;
 261 float p, q, r, s, t, y, d, err, err1;
 262 float ax, id, d1, d2, e, y1;
 263 int i, aid;
 264
 265 a = aa;
 266 b = bb;
 267 c = cc;
 268 x = xx;
 269 err = 0.0;
 270 s = 1.0 - x;
 271 if( x < -0.5 )
 272         {
 273         if( b > a )
 274                 y = powf( s, -a ) * hys2f1f( a, c-b, c, -x/s, &err );
 275
 276         else
 277                 y = powf( s, -b ) * hys2f1f( c-a, b, c, -x/s, &err );
 278
 279         goto done;
 280         }
 281
 282
 283
 284 d = c - a - b;
 285 id = roundf(d); /* nearest integer to d */
 286
 287 if( x > 0.8 )
 288 {
 289
 290 if( fabsf(d-id) > EPS2 ) /* test for integer c-a-b */
 291         {
 292 /* Try the power series first */
 293         y = hys2f1f( a, b, c, x, &err );
 294         if( err < ETHRESH )
 295                 goto done;
 296 /* If power series fails, then apply AMS55 #15.3.6 */
 297         q = hys2f1f( a, b, 1.0-d, s, &err );
 298         q *= gammaf(d) /(gammaf(c-a) * gammaf(c-b));
 299         r = powf(s,d) * hys2f1f( c-a, c-b, d+1.0, s, &err1 );
 300         r *= gammaf(-d)/(gammaf(a) * gammaf(b));
 301         y = q + r;
 302
 303         q = fabsf(q); /* estimate cancellation error */
 304         r = fabsf(r);
 305         if( q > r )
 306                 r = q;
 307         err += err1 + (MACHEPF*r)/y;
 308
 309         y *= gammaf(c);
 310         goto done;
 311         }
 312 else
 313         {
 314 /* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */
 315         if( id >= 0.0 )
 316                 {
 317                 e = d;
 318                 d1 = d;
 319                 d2 = 0.0;
 320                 aid = id;
 321                 }
 322         else
 323                 {
 324                 e = -d;
 325                 d1 = 0.0;
 326                 d2 = d;
 327                 aid = -id;
 328                 }
 329
 330         ax = logf(s);
 331
 332         /* sum for t = 0 */
 333         y = psif(1.0) + psif(1.0+e) - psif(a+d1) - psif(b+d1) - ax;
 334         y /= gammaf(e+1.0);
 335
 336         p = (a+d1) * (b+d1) * s / gammaf(e+2.0);        /* Poch for t=1 */
 337         t = 1.0;
 338         do
 339                 {
 340                 r = psif(1.0+t) + psif(1.0+t+e) - psif(a+t+d1)
 341                         - psif(b+t+d1) - ax;
 342                 q = p * r;
 343                 y += q;
 344                 p *= s * (a+t+d1) / (t+1.0);
 345                 p *= (b+t+d1) / (t+1.0+e);
 346                 t += 1.0;
 347                 }
 348         while( fabsf(q/y) > EPS );
 349
 350
 351         if( id == 0.0 )
 352                 {
 353                 y *= gammaf(c)/(gammaf(a)*gammaf(b));
 354                 goto psidon;
 355                 }
 356
 357         y1 = 1.0;
 358
 359         if( aid == 1 )
 360                 goto nosum;
 361
 362         t = 0.0;
 363         p = 1.0;
 364         for( i=1; i<aid; i++ )
 365                 {
 366                 r = 1.0-e+t;
 367                 p *= s * (a+t+d2) * (b+t+d2) / r;
 368                 t += 1.0;
 369                 p /= t;
 370                 y1 += p;
 371                 }
 372
 373
 374 nosum:
 375         p = gammaf(c);
 376         y1 *= gammaf(e) * p / (gammaf(a+d1) * gammaf(b+d1));
 377         y *= p / (gammaf(a+d2) * gammaf(b+d2));
 378         if( (aid & 1) != 0 )
 379                 y = -y;
 380
 381         q = powf( s, id );      /* s to the id power */
 382         if( id > 0.0 )
 383                 y *= q;
 384         else
 385                 y1 *= q;
 386
 387         y += y1;
 388 psidon:
 389         goto done;
 390         }
 391 }
 392
 393
 394 /* Use defining power series if no special cases */
 395 y = hys2f1f( a, b, c, x, &err );
 396
 397 done:
 398 *loss = err;
 399 return(y);
 400 }
 401
 402
 403
 404
 405
 406 /* Defining power series expansion of Gauss hypergeometric function */
 407
 408 #if 1
 409 #ifdef ANSIC
 410 static float hys2f1f( float aa, float bb, float cc, float xx, float *loss )
 411 #else
 412 static float hys2f1f( aa, bb, cc, xx, loss )
 413 double aa, bb, cc, xx;
 414 float *loss;
 415 #endif
 416 {
 417 int i;
 418 float a, b, c, x;
 419 float f, g, h, k, m, s, u, umax;
 420
 421
 422 a = aa;
 423 b = bb;
 424 c = cc;
 425 x = xx;
 426 i = 0;
 427 umax = 0.0;
 428 f = a;
 429 g = b;
 430 h = c;
 431 k = 0.0;
 432 s = 1.0;
 433 u = 1.0;
 434
 435 do
 436         {
 437         if( fabsf(h) < EPS )
 438                 return( MAXNUMF );
 439         m = k + 1.0;
 440         u = u * ((f+k) * (g+k) * x / ((h+k) * m));
 441         s += u;
 442         k = fabsf(u);  /* remember largest term summed */
 443         if( k > umax )
 444                 umax = k;
 445         k = m;
 446         if( ++i > 10000 ) /* should never happen */
 447                 {
 448                 *loss = 1.0;
 449                 return(s);
 450                 }
 451         }
 452 while( fabsf(u/s) > MACHEPF );
 453
 454 /* return estimated relative error */
 455 *loss = (MACHEPF*umax)/fabsf(s) + (MACHEPF*i);
 456
 457 return(s);
 458 }
 459
 460
 461 #else /* 0 */
 462
 463 extern double MACHEP;
 464
 465 #ifdef ANSIC
 466 static float hys2f1f( float aa, float bb, float cc, float xx, float *loss )
 467 #else
 468 static float hys2f1f( aa, bb, cc, xx, loss )
 469 double aa, bb, cc, xx;
 470 float *loss;
 471 #endif
 472 {
 473 int i;
 474 double a, b, c, x;
 475 double f, g, h, k, m, s, u, umax;
 476
 477
 478 a = aa;
 479 b = bb;
 480 c = cc;
 481 x = xx;
 482 i = 0;
 483 umax = 0.0;
 484 f = a;
 485 g = b;
 486 h = c;
 487 k = 0.0;
 488 s = 1.0;
 489 u = 1.0;
 490
 491 do
 492         {
 493         if( fabsf(h) < EPS )
 494                 {
 495                 *loss = 1.0;
 496                 return( MAXNUMF );
 497                 }
 498         m = k + 1.0;
 499         u = u * ((f+k) * (g+k) * x / ((h+k) * m));
 500         s += u;
 501         k = fabsf(u);  /* remember largest term summed */
 502         if( k > umax )
 503                 umax = k;
 504         k = m;
 505         if( ++i > 10000 ) /* should never happen */
 506                 {
 507                 *loss = 1.0;
 508                 return(s);
 509                 }
 510         }
 511 while( fabsf(u/s) > MACHEP );
 512
 513 /* return estimated relative error */
 514 *loss = (MACHEPF*umax)/fabsf(s) + (MACHEPF*i);
 515
 516 return(s);
 517 }
 518 #endif