/[gsl]/gsl/cdf/gamma.c
ViewVC logotype

Diff of /gsl/cdf/gamma.c

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

revision 1.1.2.4 by bjg, Wed Jun 25 10:36:26 2003 UTC revision 1.1.2.5 by bjg, Mon Jul 14 12:07:33 2003 UTC
# Line 27  Line 27 
27  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
28  #include <gsl/gsl_sf_gamma.h>  #include <gsl/gsl_sf_gamma.h>
29    
30  #define LARGE_B 85  #define LARGE_A 85
31    
32  /*  /*
33   * Use the normal approximation as defined in   * Use the normal approximation as defined in
# Line 42  Line 42 
42   */   */
43    
44  static double  static double
45  norm_arg (double x, double b)  norm_arg (double x, double a)
46  {  {
47    double t;    double t;
48    double arg = x + (1.0 / 3.0) - b - (0.02 / b);    double arg = x + (1.0 / 3.0) - a - (0.02 / a);
49    double u = (b - 0.5) / x;    double u = (a - 0.5) / x;
50    
51    if (fabs (u - 1.0) < GSL_DBL_EPSILON)    if (fabs (u - 1.0) < GSL_DBL_EPSILON)
52      {      {
# Line 78  double Line 78  double
78  gsl_cdf_gamma_P (double x, double a, double b)  gsl_cdf_gamma_P (double x, double a, double b)
79  {  {
80    double P;    double P;
81    double y = x / a;    double y = x / b;
82    
83    if (x <= 0.0)    if (x <= 0.0)
84      {      {
85        return 0.0;        return 0.0;
86      }      }
87    
88    if (b < LARGE_B)    if (b < LARGE_A)
89      {      {
90        P = gsl_sf_gamma_inc_P (b, y);        P = gsl_sf_gamma_inc_P (a, y);
91      }      }
92    else    else
93      {      {
94        /* Use Peizer and Pratt's normal approximation. */        /* Use Peizer and Pratt's normal approximation. */
95    
96        double z = norm_arg (y, b);        double z = norm_arg (y, a);
97        P = gsl_cdf_ugaussian_P (z);        P = gsl_cdf_ugaussian_P (z);
98      }      }
99    
# Line 104  double Line 104  double
104  gsl_cdf_gamma_Q (double x, double a, double b)  gsl_cdf_gamma_Q (double x, double a, double b)
105  {  {
106    double P;    double P;
107    double y = x / a;    double y = x / b;
108    
109    if (x <= 0.0)    if (x <= 0.0)
110      {      {
111        return 1.0;        return 1.0;
112      }      }
113    
114    if (b < LARGE_B)    if (a < LARGE_A)
115      {      {
116        P = gsl_sf_gamma_inc_Q (b, y);        P = gsl_sf_gamma_inc_Q (a, y);
117      }      }
118    else    else
119      {      {
120        /*        /*
121         * Peizer and Pratt's approximation mentioned above.         * Peizer and Pratt's approximation mentioned above.
122         */         */
123        double z = norm_arg (y, b);        double z = norm_arg (y, a);
124        P = gsl_cdf_ugaussian_Q (z);        P = gsl_cdf_ugaussian_Q (z);
125      }      }
126    

Legend:
Removed from v.1.1.2.4  
changed lines
  Added in v.1.1.2.5

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