diff options
Diffstat (limited to 'gsl-1.9/cdf/beta_inc.c')
| -rw-r--r-- | gsl-1.9/cdf/beta_inc.c | 177 |
1 files changed, 177 insertions, 0 deletions
diff --git a/gsl-1.9/cdf/beta_inc.c b/gsl-1.9/cdf/beta_inc.c new file mode 100644 index 0000000..8a0165a --- /dev/null +++ b/gsl-1.9/cdf/beta_inc.c @@ -0,0 +1,177 @@ +/* specfunc/beta_inc.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ +/* Modified for cdfs by Brian Gough, June 2003 */ + +static double +beta_cont_frac (const double a, const double b, const double x, + const double epsabs) +{ + const unsigned int max_iter = 512; /* control iterations */ + const double cutoff = 2.0 * GSL_DBL_MIN; /* control the zero cutoff */ + unsigned int iter_count = 0; + double cf; + + /* standard initialization for continued fraction */ + double num_term = 1.0; + double den_term = 1.0 - (a + b) * x / (a + 1.0); + + if (fabs (den_term) < cutoff) + den_term = GSL_NAN; + + den_term = 1.0 / den_term; + cf = den_term; + + while (iter_count < max_iter) + { + const int k = iter_count + 1; + double coeff = k * (b - k) * x / (((a - 1.0) + 2 * k) * (a + 2 * k)); + double delta_frac; + + /* first step */ + den_term = 1.0 + coeff * den_term; + num_term = 1.0 + coeff / num_term; + + if (fabs (den_term) < cutoff) + den_term = GSL_NAN; + + if (fabs (num_term) < cutoff) + num_term = GSL_NAN; + + den_term = 1.0 / den_term; + + delta_frac = den_term * num_term; + cf *= delta_frac; + + coeff = -(a + k) * (a + b + k) * x / ((a + 2 * k) * (a + 2 * k + 1.0)); + + /* second step */ + den_term = 1.0 + coeff * den_term; + num_term = 1.0 + coeff / num_term; + + if (fabs (den_term) < cutoff) + den_term = GSL_NAN; + + if (fabs (num_term) < cutoff) + num_term = GSL_NAN; + + den_term = 1.0 / den_term; + + delta_frac = den_term * num_term; + cf *= delta_frac; + + if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON) + break; + + if (cf * fabs (delta_frac - 1.0) < epsabs) + break; + + ++iter_count; + } + + if (iter_count >= max_iter) + return GSL_NAN; + + return cf; +} + +/* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x) + + Y taking account of possible cancellations when using the + hypergeometric transformation beta_inc(a,b,x)=1-beta_inc(b,a,1-x). + + It also adjusts the accuracy of beta_inc() to fit the overall + absolute error when A*beta_inc is added to Y. (e.g. if Y >> + A*beta_inc then the accuracy of beta_inc can be reduced) */ + +static double +beta_inc_AXPY (const double A, const double Y, + const double a, const double b, const double x) +{ + if (x == 0.0) + { + return A * 0 + Y; + } + else if (x == 1.0) + { + return A * 1 + Y; + } + else + { + double ln_beta = gsl_sf_lnbeta (a, b); + double ln_pre = -ln_beta + a * log (x) + b * log1p (-x); + + double prefactor = exp (ln_pre); + + if (x < (a + 1.0) / (a + b + 2.0)) + { + /* Apply continued fraction directly. */ + double epsabs = fabs (Y / (A * prefactor / a)) * GSL_DBL_EPSILON; + + double cf = beta_cont_frac (a, b, x, epsabs); + + return A * (prefactor * cf / a) + Y; + } + else + { + /* Apply continued fraction after hypergeometric transformation. */ + double epsabs = + fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON; + double cf = beta_cont_frac (b, a, 1.0 - x, epsabs); + double term = prefactor * cf / b; + + if (A == -Y) + { + return -A * term; + } + else + { + return A * (1 - term) + Y; + } + } + } +} + +/* Direct series evaluation for testing purposes only */ + +#if 0 +static double +beta_series (const double a, const double b, const double x, + const double epsabs) +{ + double f = x / (1 - x); + double c = (b - 1) / (a + 1) * f; + double s = 1; + double n = 0; + + s += c; + + do + { + n++; + c *= -f * (2 + n - b) / (2 + n + a); + s += c; + } + while (n < 512 && fabs (c) > GSL_DBL_EPSILON * fabs (s) + epsabs); + + s /= (1 - x); + + return s; +} +#endif |