/[gsl]/gsl/specfunc/gamma_inc.c
ViewVC logotype

Diff of /gsl/specfunc/gamma_inc.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.1.2 by bjg, Sun Sep 8 09:41:16 2002 UTC revision 1.1.1.2.8.1 by bjg, Sat Jun 28 18:08:40 2003 UTC
# Line 26  Line 26 
26  #include <gsl/gsl_sf_exp.h>  #include <gsl/gsl_sf_exp.h>
27  #include <gsl/gsl_sf_log.h>  #include <gsl/gsl_sf_log.h>
28  #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
29    #include <gsl/gsl_sf_expint.h>
30    
31  #include "error.h"  #include "error.h"
32    
# Line 46  gamma_inc_D(const double a, const double Line 47  gamma_inc_D(const double a, const double
47      return GSL_SUCCESS;      return GSL_SUCCESS;
48    }    }
49    else {    else {
     double mu = (x-a)/a;  
     double term1;  
50      gsl_sf_result gstar;      gsl_sf_result gstar;
51      gsl_sf_result ln_term;      gsl_sf_result ln_term;
52      gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */      double term1;
53        if (x < a) {
54          double u = x/a;
55          ln_term.val = log(u) - u + 1.0;
56          ln_term.err = ln_term.val * GSL_DBL_EPSILON;
57        } else {
58          double mu = (x-a)/a;
59          gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */
60        };
61      gsl_sf_gammastar_e(a, &gstar);      gsl_sf_gammastar_e(a, &gstar);
62      term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);      term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
63      result->val  = term1/gstar.val;      result->val  = term1/gstar.val;
# Line 58  gamma_inc_D(const double a, const double Line 65  gamma_inc_D(const double a, const double
65      result->err += gstar.err/fabs(gstar.val) * fabs(result->val);      result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
66      return GSL_SUCCESS;      return GSL_SUCCESS;
67    }    }
68        
69  }  }
70    
71    
# Line 120  gamma_inc_Q_large_x(const double a, cons Line 128  gamma_inc_Q_large_x(const double a, cons
128    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
129    
130    if(n == nmax)    if(n == nmax)
131      GSL_ERROR ("error", GSL_EMAXITER);      GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER);
132    else    else
133      return stat_D;      return stat_D;
134  }  }
# Line 169  gamma_inc_Q_asymp_unif(const double a, c Line 177  gamma_inc_Q_asymp_unif(const double a, c
177  }  }
178    
179    
180  /* Continued fraction for Q.  /* Continued fraction which occurs in evaluation
181   *   * of Q(a,x) or Gamma(a,x).
  * Q(a,x) = D(a,x) a/x F(a,x)  
  *             1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x  
  *  F(a,x) =  ---- ------- ----- -------- ----- -------- ...  
  *            1 +   1 +     1 +   1 +      1 +   1 +  
182   *   *
183   * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):   *              1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x
184     *   F(a,x) =  ---- ------- ----- -------- ----- -------- ...
185     *             1 +   1 +     1 +   1 +      1 +   1 +
186   *   *
187   * Since the Gautschi equivalent series method for CF evaluation may lead   * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no).
  * to singularities, I have replaced it with the modified Lentz algorithm  
  * given in  
188   *   *
189   * I J Thompson and A R Barnett   * Split out from gamma_inc_Q_CF() by GJ [Tue Apr  1 13:16:41 MST 2003].
190   * Coulomb and Bessel Functions of Complex Arguments and Order   * See gamma_inc_Q_CF() below.
  * J Computational Physics 64:490-509 (1986)  
  *  
  * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been  
  * removed.  
  *  
  * Identification of terms between the above equation for F(a, x) and  
  * the first equation in the appendix of Thompson&Barnett is as follows:  
  *  
  *    b_0 = 0, b_n = 1 for all n > 0  
  *  
  *    a_1 = 1  
  *    a_n = (n/2-a)/x    for n even  
  *    a_n = (n-1)/(2x)   for n odd  
191   *   *
192   */   */
193    static int
194  static  gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result)
 int  
 gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)  
195  {  {
196    const int    nmax  =  5000;    const int    nmax  =  5000;
197    const double small =  gsl_pow_3 (GSL_DBL_EPSILON);    const double small =  gsl_pow_3 (GSL_DBL_EPSILON);
198    
   gsl_sf_result D;  
   const int stat_D = gamma_inc_D(a, x, &D);  
   
199    double hn = 1.0;           /* convergent */    double hn = 1.0;           /* convergent */
200    double Cn = 1.0 / small;    double Cn = 1.0 / small;
201    double Dn = 1.0;    double Dn = 1.0;
# Line 217  gamma_inc_Q_CF(const double a, const dou Line 203  gamma_inc_Q_CF(const double a, const dou
203    
204    /* n == 1 has a_1, b_1, b_0 independent of a,x,    /* n == 1 has a_1, b_1, b_0 independent of a,x,
205       so that has been done by hand                */       so that has been done by hand                */
206    for ( n = 2 ; n < nmax ; n++ ) {    for ( n = 2 ; n < nmax ; n++ )
207      {
208      double an;      double an;
209      double delta;      double delta;
210    
# Line 235  gamma_inc_Q_CF(const double a, const dou Line 222  gamma_inc_Q_CF(const double a, const dou
222      Dn = 1.0 / Dn;      Dn = 1.0 / Dn;
223      delta = Cn * Dn;      delta = Cn * Dn;
224      hn *= delta;      hn *= delta;
225      if(fabs(delta-1) < GSL_DBL_EPSILON) break;      if(fabs(delta-1.0) < GSL_DBL_EPSILON) break;
226    }    }
227    
228    result->val  = D.val * (a/x) * hn;    result->val = hn;
229    result->err  = D.err * fabs((a/x) * hn);    result->err = 2.0*GSL_DBL_EPSILON * fabs(hn);
230    result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);    result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);
231    
232    if(n == nmax)    if(n == nmax)
233      GSL_ERROR ("error", GSL_EMAXITER);      GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
234    else    else
235      return stat_D;      return GSL_SUCCESS;
236    }
237    
238    
239    /* Continued fraction for Q.
240     *
241     * Q(a,x) = D(a,x) a/x F(a,x)
242     *
243     * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):
244     *
245     * Since the Gautschi equivalent series method for CF evaluation may lead
246     * to singularities, I have replaced it with the modified Lentz algorithm
247     * given in
248     *
249     * I J Thompson and A R Barnett
250     * Coulomb and Bessel Functions of Complex Arguments and Order
251     * J Computational Physics 64:490-509 (1986)
252     *
253     * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been
254     * removed.
255     *
256     * Identification of terms between the above equation for F(a, x) and
257     * the first equation in the appendix of Thompson&Barnett is as follows:
258     *
259     *    b_0 = 0, b_n = 1 for all n > 0
260     *
261     *    a_1 = 1
262     *    a_n = (n/2-a)/x    for n even
263     *    a_n = (n-1)/(2x)   for n odd
264     *
265     */
266    static
267    int
268    gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
269    {
270      gsl_sf_result D;
271      gsl_sf_result F;
272      const int stat_D = gamma_inc_D(a, x, &D);
273      const int stat_F = gamma_inc_F_CF(a, x, &F);
274    
275      result->val  = D.val * (a/x) * F.val;
276      result->err  = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err);
277    
278      return GSL_ERROR_SELECT_2(stat_F, stat_D);
279  }  }
280    
281    
282  /* Useful for small a and x. Handles the subtraction analytically.  /* Useful for small a and x. Handles the subtraction analytically.
283   */   */
284  static  static
# Line 271  gamma_inc_Q_series(const double a, const Line 302  gamma_inc_Q_series(const double a, const
302      const double c4 = -0.04166666666666666667      const double c4 = -0.04166666666666666667
303                         * (-1.758243446661483480 + lnx)                         * (-1.758243446661483480 + lnx)
304                         * (-0.764428657272716373 + lnx)                         * (-0.764428657272716373 + lnx)
305                         * ( 0.723980571623507657 + lnx)                         * ( 0.723980571623507657 + lnx)
306                         * ( 4.107554191916823640 + lnx);                         * ( 4.107554191916823640 + lnx);
307      const double c5 = -0.0083333333333333333      const double c5 = -0.0083333333333333333
308                         * (-2.06563396085715900 + lnx)                         * (-2.06563396085715900 + lnx)
309                         * (-1.28459889470864700 + lnx)                         * (-1.28459889470864700 + lnx)
310                         * (-0.27583535756454143 + lnx)                         * (-0.27583535756454143 + lnx)
311                         * ( 1.33677371336239618 + lnx)                         * ( 1.33677371336239618 + lnx)
312                         * ( 5.17537282427561550 + lnx);                         * ( 5.17537282427561550 + lnx);
313      const double c6 = -0.0013888888888888889      const double c6 = -0.0013888888888888889
314                         * (-2.30814336454783200 + lnx)                         * (-2.30814336454783200 + lnx)
315                         * (-1.65846557706987300 + lnx)                         * (-1.65846557706987300 + lnx)
# Line 356  gamma_inc_Q_series(const double a, const Line 387  gamma_inc_Q_series(const double a, const
387  }  }
388    
389    
390    /* series for small a and x, but not defined for a == 0 */
391    static int
392    gamma_inc_series(double a, double x, gsl_sf_result * result)
393    {
394      gsl_sf_result Q;
395      gsl_sf_result G;
396      const int stat_Q = gamma_inc_Q_series(a, x, &Q);
397      const int stat_G = gsl_sf_gamma_e(a, &G);
398      result->val = Q.val * G.val;
399      result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val);
400      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
401    
402      return GSL_ERROR_SELECT_2(stat_Q, stat_G);
403    }
404    
405    
406    static int
407    gamma_inc_a_gt_0(double a, double x, gsl_sf_result * result)
408    {
409      /* x > 0 and a > 0; use result for Q */
410      gsl_sf_result Q;
411      gsl_sf_result G;
412      const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q);
413      const int stat_G = gsl_sf_gamma_e(a, &G);
414    
415      result->val = G.val * Q.val;
416      result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val);
417      result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val);
418    
419      return GSL_ERROR_SELECT_2(stat_G, stat_Q);
420    }
421    
422    
423    static int
424    gamma_inc_CF(double a, double x, gsl_sf_result * result)
425    {
426      gsl_sf_result F;
427      gsl_sf_result pre;
428      const int stat_F = gamma_inc_F_CF(a, x, &F);
429      const int stat_E = gsl_sf_exp_e((a-1.0)*log(x) - x, &pre);
430    
431      result->val = F.val * pre.val;
432      result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
433      result->err += (2.0 + fabs(a)) * GSL_DBL_EPSILON * fabs(result->val);
434    
435      return GSL_ERROR_SELECT_2(stat_F, stat_E);
436    }
437    
438    
439    /* evaluate Gamma(0,x), x > 0 */
440    #define GAMMA_INC_A_0(x, result) gsl_sf_expint_E1_e(x, result)
441    
442    
443  /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/  /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
444    
445  int  int
446  gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)  gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
447  {  {
448    if(a <= 0.0 || x < 0.0) {    if(a < 0.0 || x < 0.0) {
449      DOMAIN_ERROR(result);      DOMAIN_ERROR(result);
450    }    }
451    else if(x == 0.0) {    else if(x == 0.0) {
# Line 369  gsl_sf_gamma_inc_Q_e(const double a, con Line 453  gsl_sf_gamma_inc_Q_e(const double a, con
453      result->err = 0.0;      result->err = 0.0;
454      return GSL_SUCCESS;      return GSL_SUCCESS;
455    }    }
456      else if(a == 0.0)
457      {
458        result->val = 0.0;
459        result->err = 0.0;
460        return GSL_SUCCESS;
461      }
462    else if(x <= 0.5*a) {    else if(x <= 0.5*a) {
463      /* If the series is quick, do that. It is      /* If the series is quick, do that. It is
464       * robust and simple.       * robust and simple.
# Line 498  gsl_sf_gamma_inc_P_e(const double a, con Line 588  gsl_sf_gamma_inc_P_e(const double a, con
588  }  }
589    
590    
591    int
592    gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result)
593    {  
594      if(x < 0.0) {
595        DOMAIN_ERROR(result);
596      }
597      else if(x == 0.0) {
598        return gsl_sf_gamma_e(a, result);
599      }
600      else if(a == 0.0)
601      {
602        return GAMMA_INC_A_0(x, result);
603      }
604      else if(a > 0.0)
605      {
606        return gamma_inc_a_gt_0(a, x, result);
607      }
608      else if(x > 0.25)
609      {
610        /* continued fraction seems to fail for x too small; otherwise
611           it is ok, independent of the value of |x/a|, because of the
612           non-oscillation in the expansion, i.e. the CF is
613           un-conditionally convergent for a < 0 and x > 0
614         */
615        return gamma_inc_CF(a, x, result);
616      }
617      else if(fabs(a) < 0.5)
618      {
619        return gamma_inc_series(a, x, result);
620      }
621      else
622      {
623        /* a = fa + da; da >= 0 */
624        const double fa = floor(a);
625        const double da = a - fa;
626    
627        gsl_sf_result g_da;
628        const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da)
629                                         : GAMMA_INC_A_0(x, &g_da));
630    
631        double alpha = da;
632        double gax = g_da.val;
633    
634        /* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */
635        do
636        {
637          const double shift = exp(-x + (alpha-1.0)*log(x));
638          gax = (gax - shift) / (alpha - 1.0);
639          alpha -= 1.0;
640        } while(alpha > a);
641    
642        result->val = gax;
643        result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax);
644        return stat_g_da;
645      }
646    
647    }
648    
649    
650  /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/  /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
651    
652  #include "eval.h"  #include "eval.h"
# Line 511  double gsl_sf_gamma_inc_Q(const double a Line 660  double gsl_sf_gamma_inc_Q(const double a
660  {  {
661    EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));    EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
662  }  }
663    
664    double gsl_sf_gamma_inc(const double a, const double x)
665    {
666       EVAL_RESULT(gsl_sf_gamma_inc_e(a, x, &result));
667    }

Legend:
Removed from v.1.1.1.2  
changed lines
  Added in v.1.1.1.2.8.1

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26