htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/gamma_inc.c

   1 /* specfunc/gamma_inc.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /* Author:  G. Jungman */
  22
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_erf.h>
  27 #include <gsl/gsl_sf_exp.h>
  28 #include <gsl/gsl_sf_log.h>
  29 #include <gsl/gsl_sf_gamma.h>
  30 #include <gsl/gsl_sf_expint.h>
  31
  32 #include "error.h"
  33
  34 /* The dominant part,
  35  * D(a,x) := x^a e^(-x) / Gamma(a+1)
  36  */
  37 static
  38 int
  39 gamma_inc_D(const double a, const double x, gsl_sf_result * result)
  40 {
  41   if(a < 10.0) {
  42     double lnr;
  43     gsl_sf_result lg;
  44     gsl_sf_lngamma_e(a+1.0, &lg);
  45     lnr = a * log(x) - x - lg.val;
  46     result->val = exp(lnr);
  47     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
  48     return GSL_SUCCESS;
  49   }
  50   else {
  51     gsl_sf_result gstar;
  52     gsl_sf_result ln_term;
  53     double term1;
  54     if (x < 0.5*a) {
  55       double u = x/a;
  56       double ln_u = log(u);
  57       ln_term.val = ln_u - u + 1.0;
  58       ln_term.err = (fabs(ln_u) + fabs(u) + 1.0) * GSL_DBL_EPSILON;
  59     } else {
  60       double mu = (x-a)/a;
  61       gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */
  62     };
  63     gsl_sf_gammastar_e(a, &gstar);
  64     term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
  65     result->val  = term1/gstar.val;
  66     result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
  67     result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
  68     return GSL_SUCCESS;
  69   }
  70
  71 }
  72
  73
  74 /* P series representation.
  75  */
  76 static
  77 int
  78 gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
  79 {
  80   const int nmax = 5000;
  81
  82   gsl_sf_result D;
  83   int stat_D = gamma_inc_D(a, x, &D);
  84
  85   double sum  = 1.0;
  86   double term = 1.0;
  87   int n;
  88   for(n=1; n<nmax; n++) {
  89     term *= x/(a+n);
  90     sum  += term;
  91     if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  92   }
  93
  94   result->val  = D.val * sum;
  95   result->err  = D.err * fabs(sum);
  96   result->err += (1.0 + n) * GSL_DBL_EPSILON * fabs(result->val);
  97
  98   if(n == nmax)
  99     GSL_ERROR ("error", GSL_EMAXITER);
 100   else
 101     return stat_D;
 102 }
 103
 104
 105 /* Q large x asymptotic
 106  */
 107 static
 108 int
 109 gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
 110 {
 111   const int nmax = 5000;
 112
 113   gsl_sf_result D;
 114   const int stat_D = gamma_inc_D(a, x, &D);
 115
 116   double sum  = 1.0;
 117   double term = 1.0;
 118   double last = 1.0;
 119   int n;
 120   for(n=1; n<nmax; n++) {
 121     term *= (a-n)/x;
 122     if(fabs(term/last) > 1.0) break;
 123     if(fabs(term/sum)  < GSL_DBL_EPSILON) break;
 124     sum  += term;
 125     last  = term;
 126   }
 127
 128   result->val  = D.val * (a/x) * sum;
 129   result->err  = D.err * fabs((a/x) * sum);
 130   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 131
 132   if(n == nmax)
 133     GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER);
 134   else
 135     return stat_D;
 136 }
 137
 138
 139 /* Uniform asymptotic for x near a, a and x large.
 140  * See [Temme, p. 285]
 141  */
 142 static
 143 int
 144 gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
 145 {
 146   const double rta = sqrt(a);
 147   const double eps = (x-a)/a;
 148
 149   gsl_sf_result ln_term;
 150   const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term);  /* log(1+eps) - eps */
 151   const double eta  = GSL_SIGN(eps) * sqrt(-2.0*ln_term.val);
 152
 153   gsl_sf_result erfc;
 154
 155   double R;
 156   double c0, c1;
 157
 158   /* This used to say erfc(eta*M_SQRT2*rta), which is wrong.
 159    * The sqrt(2) is in the denominator. Oops.
 160    * Fixed: [GJ] Mon Nov 15 13:25:32 MST 2004
 161    */
 162   gsl_sf_erfc_e(eta*rta/M_SQRT2, &erfc);
 163
 164   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
 165     c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0)));
 166     c1 = -1.0/540.0 - eps/288.0;
 167   }
 168   else {
 169     const double rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
 170     const double lam = x/a;
 171     c0 = (1.0 - 1.0/rt_term)/eps;
 172     c1 = -(eta*eta*eta * (lam*lam + 10.0*lam + 1.0) - 12.0 * eps*eps*eps) / (12.0 * eta*eta*eta*eps*eps*eps);
 173   }
 174
 175   R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);
 176
 177   result->val  = 0.5 * erfc.val + R;
 178   result->err  = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
 179   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 180
 181   return stat_ln;
 182 }
 183
 184
 185 /* Continued fraction which occurs in evaluation
 186  * of Q(a,x) or Gamma(a,x).
 187  *
 188  *              1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x
 189  *   F(a,x) =  ---- ------- ----- -------- ----- -------- ...
 190  *             1 +   1 +     1 +   1 +      1 +   1 +
 191  *
 192  * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no).
 193  *
 194  * Split out from gamma_inc_Q_CF() by GJ [Tue Apr  1 13:16:41 MST 2003].
 195  * See gamma_inc_Q_CF() below.
 196  *
 197  */
 198 static int
 199 gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result)
 200 {
 201   const int    nmax  =  5000;
 202   const double small =  gsl_pow_3 (GSL_DBL_EPSILON);
 203
 204   double hn = 1.0;           /* convergent */
 205   double Cn = 1.0 / small;
 206   double Dn = 1.0;
 207   int n;
 208
 209   /* n == 1 has a_1, b_1, b_0 independent of a,x,
 210      so that has been done by hand                */
 211   for ( n = 2 ; n < nmax ; n++ )
 212   {
 213     double an;
 214     double delta;
 215
 216     if(GSL_IS_ODD(n))
 217       an = 0.5*(n-1)/x;
 218     else
 219       an = (0.5*n-a)/x;
 220
 221     Dn = 1.0 + an * Dn;
 222     if ( fabs(Dn) < small )
 223       Dn = small;
 224     Cn = 1.0 + an/Cn;
 225     if ( fabs(Cn) < small )
 226       Cn = small;
 227     Dn = 1.0 / Dn;
 228     delta = Cn * Dn;
 229     hn *= delta;
 230     if(fabs(delta-1.0) < GSL_DBL_EPSILON) break;
 231   }
 232
 233   result->val = hn;
 234   result->err = 2.0*GSL_DBL_EPSILON * fabs(hn);
 235   result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);
 236
 237   if(n == nmax)
 238     GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
 239   else
 240     return GSL_SUCCESS;
 241 }
 242
 243
 244 /* Continued fraction for Q.
 245  *
 246  * Q(a,x) = D(a,x) a/x F(a,x)
 247  *
 248  * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):
 249  *
 250  * Since the Gautschi equivalent series method for CF evaluation may lead
 251  * to singularities, I have replaced it with the modified Lentz algorithm
 252  * given in
 253  *
 254  * I J Thompson and A R Barnett
 255  * Coulomb and Bessel Functions of Complex Arguments and Order
 256  * J Computational Physics 64:490-509 (1986)
 257  *
 258  * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been
 259  * removed.
 260  *
 261  * Identification of terms between the above equation for F(a, x) and
 262  * the first equation in the appendix of Thompson&Barnett is as follows:
 263  *
 264  *    b_0 = 0, b_n = 1 for all n > 0
 265  *
 266  *    a_1 = 1
 267  *    a_n = (n/2-a)/x    for n even
 268  *    a_n = (n-1)/(2x)   for n odd
 269  *
 270  */
 271 static
 272 int
 273 gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
 274 {
 275   gsl_sf_result D;
 276   gsl_sf_result F;
 277   const int stat_D = gamma_inc_D(a, x, &D);
 278   const int stat_F = gamma_inc_F_CF(a, x, &F);
 279
 280   result->val  = D.val * (a/x) * F.val;
 281   result->err  = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err);
 282
 283   return GSL_ERROR_SELECT_2(stat_F, stat_D);
 284 }
 285
 286
 287 /* Useful for small a and x. Handles the subtraction analytically.
 288  */
 289 static
 290 int
 291 gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
 292 {
 293   double term1;  /* 1 - x^a/Gamma(a+1) */
 294   double sum;    /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
 295   int stat_sum;
 296   double term2;  /* a temporary variable used at the end */
 297
 298   {
 299     /* Evaluate series for 1 - x^a/Gamma(a+1), small a
 300      */
 301     const double pg21 = -2.404113806319188570799476;  /* PolyGamma[2,1] */
 302     const double lnx  = log(x);
 303     const double el   = M_EULER+lnx;
 304     const double c1 = -el;
 305     const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
 306     const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
 307     const double c4 = -0.04166666666666666667
 308                        * (-1.758243446661483480 + lnx)
 309                        * (-0.764428657272716373 + lnx)
 310                        * ( 0.723980571623507657 + lnx)
 311                        * ( 4.107554191916823640 + lnx);
 312     const double c5 = -0.0083333333333333333
 313                        * (-2.06563396085715900 + lnx)
 314                        * (-1.28459889470864700 + lnx)
 315                        * (-0.27583535756454143 + lnx)
 316                        * ( 1.33677371336239618 + lnx)
 317                        * ( 5.17537282427561550 + lnx);
 318     const double c6 = -0.0013888888888888889
 319                        * (-2.30814336454783200 + lnx)
 320                        * (-1.65846557706987300 + lnx)
 321                        * (-0.88768082560020400 + lnx)
 322                        * ( 0.17043847751371778 + lnx)
 323                        * ( 1.92135970115863890 + lnx)
 324                        * ( 6.22578557795474900 + lnx);
 325     const double c7 = -0.00019841269841269841
 326                        * (-2.5078657901291800 + lnx)
 327                        * (-1.9478900888958200 + lnx)
 328                        * (-1.3194837322612730 + lnx)
 329                        * (-0.5281322700249279 + lnx)
 330                        * ( 0.5913834939078759 + lnx)
 331                        * ( 2.4876819633378140 + lnx)
 332                        * ( 7.2648160783762400 + lnx);
 333     const double c8 = -0.00002480158730158730
 334                        * (-2.677341544966400 + lnx)
 335                        * (-2.182810448271700 + lnx)
 336                        * (-1.649350342277400 + lnx)
 337                        * (-1.014099048290790 + lnx)
 338                        * (-0.191366955370652 + lnx)
 339                        * ( 0.995403817918724 + lnx)
 340                        * ( 3.041323283529310 + lnx)
 341                        * ( 8.295966556941250 + lnx);
 342     const double c9 = -2.75573192239859e-6
 343                        * (-2.8243487670469080 + lnx)
 344                        * (-2.3798494322701120 + lnx)
 345                        * (-1.9143674728689960 + lnx)
 346                        * (-1.3814529102920370 + lnx)
 347                        * (-0.7294312810261694 + lnx)
 348                        * ( 0.1299079285269565 + lnx)
 349                        * ( 1.3873333251885240 + lnx)
 350                        * ( 3.5857258865210760 + lnx)
 351                        * ( 9.3214237073814600 + lnx);
 352     const double c10 = -2.75573192239859e-7
 353                        * (-2.9540329644556910 + lnx)
 354                        * (-2.5491366926991850 + lnx)
 355                        * (-2.1348279229279880 + lnx)
 356                        * (-1.6741881076349450 + lnx)
 357                        * (-1.1325949616098420 + lnx)
 358                        * (-0.4590034650618494 + lnx)
 359                        * ( 0.4399352987435699 + lnx)
 360                        * ( 1.7702236517651670 + lnx)
 361                        * ( 4.1231539047474080 + lnx)
 362                        * ( 10.342627908148680 + lnx);
 363
 364     term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
 365   }
 366
 367   {
 368     /* Evaluate the sum.
 369      */
 370     const int nmax = 5000;
 371     double t = 1.0;
 372     int n;
 373     sum = 1.0;
 374
 375     for(n=1; n<nmax; n++) {
 376       t *= -x/(n+1.0);
 377       sum += (a+1.0)/(a+n+1.0)*t;
 378       if(fabs(t/sum) < GSL_DBL_EPSILON) break;
 379     }
 380
 381     if(n == nmax)
 382       stat_sum = GSL_EMAXITER;
 383     else
 384       stat_sum = GSL_SUCCESS;
 385   }
 386
 387   term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
 388   result->val  = term1 + term2;
 389   result->err  = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
 390   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 391   return stat_sum;
 392 }
 393
 394
 395 /* series for small a and x, but not defined for a == 0 */
 396 static int
 397 gamma_inc_series(double a, double x, gsl_sf_result * result)
 398 {
 399   gsl_sf_result Q;
 400   gsl_sf_result G;
 401   const int stat_Q = gamma_inc_Q_series(a, x, &Q);
 402   const int stat_G = gsl_sf_gamma_e(a, &G);
 403   result->val = Q.val * G.val;
 404   result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val);
 405   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 406
 407   return GSL_ERROR_SELECT_2(stat_Q, stat_G);
 408 }
 409
 410
 411 static int
 412 gamma_inc_a_gt_0(double a, double x, gsl_sf_result * result)
 413 {
 414   /* x > 0 and a > 0; use result for Q */
 415   gsl_sf_result Q;
 416   gsl_sf_result G;
 417   const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q);
 418   const int stat_G = gsl_sf_gamma_e(a, &G);
 419
 420   result->val = G.val * Q.val;
 421   result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val);
 422   result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val);
 423
 424   return GSL_ERROR_SELECT_2(stat_G, stat_Q);
 425 }
 426
 427
 428 static int
 429 gamma_inc_CF(double a, double x, gsl_sf_result * result)
 430 {
 431   gsl_sf_result F;
 432   gsl_sf_result pre;
 433   const double am1lgx = (a-1.0)*log(x);
 434   const int stat_F = gamma_inc_F_CF(a, x, &F);
 435   const int stat_E = gsl_sf_exp_err_e(am1lgx - x, GSL_DBL_EPSILON*fabs(am1lgx), &pre);
 436
 437   result->val = F.val * pre.val;
 438   result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
 439   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 440
 441   return GSL_ERROR_SELECT_2(stat_F, stat_E);
 442 }
 443
 444
 445 /* evaluate Gamma(0,x), x > 0 */
 446 #define GAMMA_INC_A_0(x, result) gsl_sf_expint_E1_e(x, result)
 447
 448
 449 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 450
 451 int
 452 gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
 453 {
 454   if(a < 0.0 || x < 0.0) {
 455     DOMAIN_ERROR(result);
 456   }
 457   else if(x == 0.0) {
 458     result->val = 1.0;
 459     result->err = 0.0;
 460     return GSL_SUCCESS;
 461   }
 462   else if(a == 0.0)
 463   {
 464     result->val = 0.0;
 465     result->err = 0.0;
 466     return GSL_SUCCESS;
 467   }
 468   else if(x <= 0.5*a) {
 469     /* If the series is quick, do that. It is
 470      * robust and simple.
 471      */
 472     gsl_sf_result P;
 473     int stat_P = gamma_inc_P_series(a, x, &P);
 474     result->val  = 1.0 - P.val;
 475     result->err  = P.err;
 476     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 477     return stat_P;
 478   }
 479   else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
 480     /* Then try the difficult asymptotic regime.
 481      * This is the only way to do this region.
 482      */
 483     return gamma_inc_Q_asymp_unif(a, x, result);
 484   }
 485   else if(a < 0.2 && x < 5.0) {
 486     /* Cancellations at small a must be handled
 487      * analytically; x should not be too big
 488      * either since the series terms grow
 489      * with x and log(x).
 490      */
 491     return gamma_inc_Q_series(a, x, result);
 492   }
 493   else if(a <= x) {
 494     if(x <= 1.0e+06) {
 495       /* Continued fraction is excellent for x >~ a.
 496        * We do not let x be too large when x > a since
 497        * it is somewhat pointless to try this there;
 498        * the function is rapidly decreasing for
 499        * x large and x > a, and it will just
 500        * underflow in that region anyway. We
 501        * catch that case in the standard
 502        * large-x method.
 503        */
 504       return gamma_inc_Q_CF(a, x, result);
 505     }
 506     else {
 507       return gamma_inc_Q_large_x(a, x, result);
 508     }
 509   }
 510   else {
 511     if(x > a - sqrt(a)) {
 512       /* Continued fraction again. The convergence
 513        * is a little slower here, but that is fine.
 514        * We have to trade that off against the slow
 515        * convergence of the series, which is the
 516        * only other option.
 517        */
 518       return gamma_inc_Q_CF(a, x, result);
 519     }
 520     else {
 521       gsl_sf_result P;
 522       int stat_P = gamma_inc_P_series(a, x, &P);
 523       result->val  = 1.0 - P.val;
 524       result->err  = P.err;
 525       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 526       return stat_P;
 527     }
 528   }
 529 }
 530
 531
 532 int
 533 gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
 534 {
 535   if(a <= 0.0 || x < 0.0) {
 536     DOMAIN_ERROR(result);
 537   }
 538   else if(x == 0.0) {
 539     result->val = 0.0;
 540     result->err = 0.0;
 541     return GSL_SUCCESS;
 542   }
 543   else if(x < 20.0 || x < 0.5*a) {
 544     /* Do the easy series cases. Robust and quick.
 545      */
 546     return gamma_inc_P_series(a, x, result);
 547   }
 548   else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
 549     /* Crossover region. Note that Q and P are
 550      * roughly the same order of magnitude here,
 551      * so the subtraction is stable.
 552      */
 553     gsl_sf_result Q;
 554     int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
 555     result->val  = 1.0 - Q.val;
 556     result->err  = Q.err;
 557     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 558     return stat_Q;
 559   }
 560   else if(a <= x) {
 561     /* Q <~ P in this area, so the
 562      * subtractions are stable.
 563      */
 564     gsl_sf_result Q;
 565     int stat_Q;
 566     if(a > 0.2*x) {
 567       stat_Q = gamma_inc_Q_CF(a, x, &Q);
 568     }
 569     else {
 570       stat_Q = gamma_inc_Q_large_x(a, x, &Q);
 571     }
 572     result->val  = 1.0 - Q.val;
 573     result->err  = Q.err;
 574     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 575     return stat_Q;
 576   }
 577   else {
 578     if((x-a)*(x-a) < a) {
 579       /* This condition is meant to insure
 580        * that Q is not very close to 1,
 581        * so the subtraction is stable.
 582        */
 583       gsl_sf_result Q;
 584       int stat_Q = gamma_inc_Q_CF(a, x, &Q);
 585       result->val  = 1.0 - Q.val;
 586       result->err  = Q.err;
 587       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 588       return stat_Q;
 589     }
 590     else {
 591       return gamma_inc_P_series(a, x, result);
 592     }
 593   }
 594 }
 595
 596
 597 int
 598 gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result)
 599 {
 600   if(x < 0.0) {
 601     DOMAIN_ERROR(result);
 602   }
 603   else if(x == 0.0) {
 604     return gsl_sf_gamma_e(a, result);
 605   }
 606   else if(a == 0.0)
 607   {
 608     return GAMMA_INC_A_0(x, result);
 609   }
 610   else if(a > 0.0)
 611   {
 612     return gamma_inc_a_gt_0(a, x, result);
 613   }
 614   else if(x > 0.25)
 615   {
 616     /* continued fraction seems to fail for x too small; otherwise
 617        it is ok, independent of the value of |x/a|, because of the
 618        non-oscillation in the expansion, i.e. the CF is
 619        un-conditionally convergent for a < 0 and x > 0
 620      */
 621     return gamma_inc_CF(a, x, result);
 622   }
 623   else if(fabs(a) < 0.5)
 624   {
 625     return gamma_inc_series(a, x, result);
 626   }
 627   else
 628   {
 629     /* a = fa + da; da >= 0 */
 630     const double fa = floor(a);
 631     const double da = a - fa;
 632
 633     gsl_sf_result g_da;
 634     const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da)
 635                                      : GAMMA_INC_A_0(x, &g_da));
 636
 637     double alpha = da;
 638     double gax = g_da.val;
 639
 640     /* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */
 641     do
 642     {
 643       const double shift = exp(-x + (alpha-1.0)*log(x));
 644       gax = (gax - shift) / (alpha - 1.0);
 645       alpha -= 1.0;
 646     } while(alpha > a);
 647
 648     result->val = gax;
 649     result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax);
 650     return stat_g_da;
 651   }
 652
 653 }
 654
 655
 656 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 657
 658 #include "eval.h"
 659
 660 double gsl_sf_gamma_inc_P(const double a, const double x)
 661 {
 662   EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
 663 }
 664
 665 double gsl_sf_gamma_inc_Q(const double a, const double x)
 666 {
 667   EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
 668 }
 669
 670 double gsl_sf_gamma_inc(const double a, const double x)
 671 {
 672    EVAL_RESULT(gsl_sf_gamma_inc_e(a, x, &result));
 673 }