htswanalysis/MACS/lib/gsl/gsl-1.11/cdf/beta_inc.c

   1 /* specfunc/beta_inc.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21 /* Modified for cdfs by Brian Gough, June 2003 */
  22
  23 #include <gsl/gsl_sf_gamma.h>
  24
  25 static double
  26 beta_cont_frac (const double a, const double b, const double x,
  27                 const double epsabs)
  28 {
  29   const unsigned int max_iter = 512;    /* control iterations      */
  30   const double cutoff = 2.0 * GSL_DBL_MIN;      /* control the zero cutoff */
  31   unsigned int iter_count = 0;
  32   double cf;
  33
  34   /* standard initialization for continued fraction */
  35   double num_term = 1.0;
  36   double den_term = 1.0 - (a + b) * x / (a + 1.0);
  37
  38   if (fabs (den_term) < cutoff)
  39     den_term = GSL_NAN;
  40
  41   den_term = 1.0 / den_term;
  42   cf = den_term;
  43
  44   while (iter_count < max_iter)
  45     {
  46       const int k = iter_count + 1;
  47       double coeff = k * (b - k) * x / (((a - 1.0) + 2 * k) * (a + 2 * k));
  48       double delta_frac;
  49
  50       /* first step */
  51       den_term = 1.0 + coeff * den_term;
  52       num_term = 1.0 + coeff / num_term;
  53
  54       if (fabs (den_term) < cutoff)
  55         den_term = GSL_NAN;
  56
  57       if (fabs (num_term) < cutoff)
  58         num_term = GSL_NAN;
  59
  60       den_term = 1.0 / den_term;
  61
  62       delta_frac = den_term * num_term;
  63       cf *= delta_frac;
  64
  65       coeff = -(a + k) * (a + b + k) * x / ((a + 2 * k) * (a + 2 * k + 1.0));
  66
  67       /* second step */
  68       den_term = 1.0 + coeff * den_term;
  69       num_term = 1.0 + coeff / num_term;
  70
  71       if (fabs (den_term) < cutoff)
  72         den_term = GSL_NAN;
  73
  74       if (fabs (num_term) < cutoff)
  75         num_term = GSL_NAN;
  76
  77       den_term = 1.0 / den_term;
  78
  79       delta_frac = den_term * num_term;
  80       cf *= delta_frac;
  81
  82       if (fabs (delta_frac - 1.0) < 2.0 * GSL_DBL_EPSILON)
  83         break;
  84
  85       if (cf * fabs (delta_frac - 1.0) < epsabs)
  86         break;
  87
  88       ++iter_count;
  89     }
  90
  91   if (iter_count >= max_iter)
  92     return GSL_NAN;
  93
  94   return cf;
  95 }
  96
  97 /* The function beta_inc_AXPY(A,Y,a,b,x) computes A * beta_inc(a,b,x)
  98    + Y taking account of possible cancellations when using the
  99    hypergeometric transformation beta_inc(a,b,x)=1-beta_inc(b,a,1-x).
 100
 101    It also adjusts the accuracy of beta_inc() to fit the overall
 102    absolute error when A*beta_inc is added to Y. (e.g. if Y >>
 103    A*beta_inc then the accuracy of beta_inc can be reduced) */
 104
 105
 106
 107 static double
 108 beta_inc_AXPY (const double A, const double Y,
 109                const double a, const double b, const double x)
 110 {
 111   if (x == 0.0)
 112     {
 113       return A * 0 + Y;
 114     }
 115   else if (x == 1.0)
 116     {
 117       return A * 1 + Y;
 118     }
 119   else if (a > 1e5 && b < 10 && x > a / (a + b))
 120     {
 121       /* Handle asymptotic regime, large a, small b, x > peak [AS 26.5.17] */
 122       double N = a + (b - 1.0) / 2.0;
 123       return A * gsl_sf_gamma_inc_Q (b, -N * log (x)) + Y;
 124     }
 125   else if (b > 1e5 && a < 10 && x < b / (a + b))
 126     {
 127       /* Handle asymptotic regime, small a, large b, x < peak [AS 26.5.17] */
 128       double N = b + (a - 1.0) / 2.0;
 129       return A * gsl_sf_gamma_inc_P (a, -N * log1p (-x)) + Y;
 130     }
 131   else
 132     {
 133       double ln_beta = gsl_sf_lnbeta (a, b);
 134       double ln_pre = -ln_beta + a * log (x) + b * log1p (-x);
 135
 136       double prefactor = exp (ln_pre);
 137
 138       if (x < (a + 1.0) / (a + b + 2.0))
 139         {
 140           /* Apply continued fraction directly. */
 141           double epsabs = fabs (Y / (A * prefactor / a)) * GSL_DBL_EPSILON;
 142
 143           double cf = beta_cont_frac (a, b, x, epsabs);
 144
 145           return A * (prefactor * cf / a) + Y;
 146         }
 147       else
 148         {
 149           /* Apply continued fraction after hypergeometric transformation. */
 150           double epsabs =
 151             fabs ((A + Y) / (A * prefactor / b)) * GSL_DBL_EPSILON;
 152           double cf = beta_cont_frac (b, a, 1.0 - x, epsabs);
 153           double term = prefactor * cf / b;
 154
 155           if (A == -Y)
 156             {
 157               return -A * term;
 158             }
 159           else
 160             {
 161               return A * (1 - term) + Y;
 162             }
 163         }
 164     }
 165 }
 166
 167 /* Direct series evaluation for testing purposes only */
 168
 169 #if 0
 170 static double
 171 beta_series (const double a, const double b, const double x,
 172              const double epsabs)
 173 {
 174   double f = x / (1 - x);
 175   double c = (b - 1) / (a + 1) * f;
 176   double s = 1;
 177   double n = 0;
 178
 179   s += c;
 180
 181   do
 182     {
 183       n++;
 184       c *= -f * (2 + n - b) / (2 + n + a);
 185       s += c;
 186     }
 187   while (n < 512 && fabs (c) > GSL_DBL_EPSILON * fabs (s) + epsabs);
 188
 189   s /= (1 - x);
 190
 191   return s;
 192 }
 193 #endif