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

Diff of /gsl/cdf/beta_inc.c

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

revision 1.1.2.2 by bjg, Sun Jun 22 21:20:36 2003 UTC revision 1.1.2.3 by bjg, Mon Jun 23 14:05:03 2003 UTC
# Line 21  Line 21 
21  /* Modified for cdfs by Brian Gough, June 2003 */  /* Modified for cdfs by Brian Gough, June 2003 */
22    
23  static double  static double
24  beta_cont_frac (const double a, const double b, const double x)  beta_cont_frac (const double a, const double b, const double x,
25                    const double epsabs)
26  {  {
27    const unsigned int max_iter = 512;    /* control iterations      */    const unsigned int max_iter = 512;    /* control iterations      */
28    const double cutoff = 2.0 * GSL_DBL_MIN;      /* control the zero cutoff */    const double cutoff = 2.0 * GSL_DBL_MIN;      /* control the zero cutoff */
# Line 33  beta_cont_frac (const double a, const do Line 34  beta_cont_frac (const double a, const do
34    double den_term = 1.0 - (a + b) * x / (a + 1.0);    double den_term = 1.0 - (a + b) * x / (a + 1.0);
35    
36    if (fabs (den_term) < cutoff)    if (fabs (den_term) < cutoff)
37      den_term = cutoff;      den_term = GSL_NAN;
38    
39    den_term = 1.0 / den_term;    den_term = 1.0 / den_term;
40    cf = den_term;    cf = den_term;
# Line 47  beta_cont_frac (const double a, const do Line 48  beta_cont_frac (const double a, const do
48        /* first step */        /* first step */
49        den_term = 1.0 + coeff * den_term;        den_term = 1.0 + coeff * den_term;
50        num_term = 1.0 + coeff / num_term;        num_term = 1.0 + coeff / num_term;
51    
52        if (fabs (den_term) < cutoff)        if (fabs (den_term) < cutoff)
53          den_term = cutoff;          den_term = GSL_NAN;
54    
55        if (fabs (num_term) < cutoff)        if (fabs (num_term) < cutoff)
56          num_term = cutoff;          num_term = GSL_NAN;
57    
58        den_term = 1.0 / den_term;        den_term = 1.0 / den_term;
59    
60        delta_frac = den_term * num_term;        delta_frac = den_term * num_term;
# Line 61  beta_cont_frac (const double a, const do Line 65  beta_cont_frac (const double a, const do
65        /* second step */        /* second step */
66        den_term = 1.0 + coeff * den_term;        den_term = 1.0 + coeff * den_term;
67        num_term = 1.0 + coeff / num_term;        num_term = 1.0 + coeff / num_term;
68    
69        if (fabs (den_term) < cutoff)        if (fabs (den_term) < cutoff)
70          den_term = cutoff;          den_term = GSL_NAN;
71    
72        if (fabs (num_term) < cutoff)        if (fabs (num_term) < cutoff)
73          num_term = cutoff;          num_term = GSL_NAN;
74    
75        den_term = 1.0 / den_term;        den_term = 1.0 / den_term;
76    
# Line 74  beta_cont_frac (const double a, const do Line 80  beta_cont_frac (const double a, const do
80        if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON)        if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON)
81          break;          break;
82    
83          if (cf * fabs (delta_frac - 1.0) < epsabs)
84            break;
85    
86        ++iter_count;        ++iter_count;
87      }      }
88    
# Line 83  beta_cont_frac (const double a, const do Line 92  beta_cont_frac (const double a, const do
92    return cf;    return cf;
93  }  }
94    
95    /* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x)
96       + Y taking account of possible cancellations when using the
97  /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/     hypergeometric transformation beta_inc(a,b,x)=1-beta(b,a,1-x).
98    
99       It also adjusts the accuracy of beta_inc() to fit the overall
100       absolute error when A*beta_inc is added to Y. (e.g. if Y >>
101       A*beta_inc then the accuracy of beta_inc can be reduced) */
102    
103  static double  static double
104  beta_inc (const double a, const double b, const double x)  beta_inc_AXPY (const double A, const double Y,
105                   const double a, const double b, const double x)
106  {  {
107    if (x == 0.0)    if (x == 0.0)
108      {      {
109        return 0;        return A * 0 + Y;
110      }      }
111    else if (x == 1.0)    else if (x == 1.0)
112      {      {
113        return 1;        return A * 1 + Y;
114      }      }
115    else    else
116      {      {
117        double ln_beta = gsl_sf_lnbeta (a, b);        double ln_beta = gsl_sf_lnbeta (a, b);
118        double ln_pre = -ln_beta + a * log (x) + b * log1p (-x);        double ln_pre = -ln_beta + a * log (x) + b * log1p (-x);
119    
120        double prefactor = exp(ln_pre);        double prefactor = exp (ln_pre);
121    
122        if (x < (a + 1.0) / (a + b + 2.0))        if (x < (a + 1.0) / (a + b + 2.0))
123          {          {
124            /* Apply continued fraction directly. */            /* Apply continued fraction directly. */
125            double cf = beta_cont_frac (a, b, x);            double epsabs = fabs (Y / (A * prefactor / a)) * GSL_DBL_EPSILON;
126    
127            return prefactor * cf / a;            double cf = beta_cont_frac (a, b, x, epsabs);
128    
129              return A * (prefactor * cf / a) + Y;
130          }          }
131        else        else
132          {          {
133            /* Apply continued fraction after hypergeometric transformation. */            /* Apply continued fraction after hypergeometric transformation. */
134            double cf = beta_cont_frac (b, a, 1.0 - x);            double epsabs =
135                fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON;
136              double cf = beta_cont_frac (b, a, 1.0 - x, epsabs);
137            double term = prefactor * cf / b;            double term = prefactor * cf / b;
138    
139            return 1 - term;            if (A == -Y)
140                {
141                  return -A * term;
142                }
143              else
144                {
145                  return A * (1 - term) + Y;
146                }
147          }          }
148      }      }
149  }  }
150    
151    /* Direct series evaluation for testing purposes only */
152    
153    #if 0
154    static double
155    beta_series (const double a, const double b, const double x,
156                 const double epsabs)
157    {
158      double f = x / (1 - x);
159      double c = (b - 1) / (a + 1) * f;
160      double s = 1;
161      double n = 0;
162    
163      s += c;
164    
165      do
166        {
167          n++;
168          c *= -f * (2 + n - b) / (2 + n + a);
169          s += c;
170        }
171      while (n < 512 && fabs (c) > GSL_DBL_EPSILON * fabs (s) + epsabs);
172    
173      s /= (1 - x);
174    
175      return s;
176    }
177    #endif
178    

Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.1.2.3

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