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

   1 /* specfunc/beta_inc.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 /* Author:  G. Jungman */
  22
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_log.h>
  27 #include <gsl/gsl_sf_exp.h>
  28 #include <gsl/gsl_sf_gamma.h>
  29 #include <gsl/gsl_sf_hyperg.h>
  30
  31 #include "error.h"
  32 #include "check.h"
  33
  34 static double
  35 isnegint (const double x)
  36 {
  37   return (x < 0) && (x == floor(x));
  38 }
  39
  40 static
  41 int
  42 beta_cont_frac(
  43   const double a,
  44   const double b,
  45   const double x,
  46   gsl_sf_result * result
  47   )
  48 {
  49   const unsigned int max_iter = 512;        /* control iterations      */
  50   const double cutoff = 2.0 * GSL_DBL_MIN;  /* control the zero cutoff */
  51   unsigned int iter_count = 0;
  52   double cf;
  53
  54   /* standard initialization for continued fraction */
  55   double num_term = 1.0;
  56   double den_term = 1.0 - (a+b)*x/(a+1.0);
  57   if (fabs(den_term) < cutoff) den_term = cutoff;
  58   den_term = 1.0/den_term;
  59   cf = den_term;
  60
  61   while(iter_count < max_iter) {
  62     const int k  = iter_count + 1;
  63     double coeff = k*(b-k)*x/(((a-1.0)+2*k)*(a+2*k));
  64     double delta_frac;
  65
  66     /* first step */
  67     den_term = 1.0 + coeff*den_term;
  68     num_term = 1.0 + coeff/num_term;
  69     if(fabs(den_term) < cutoff) den_term = cutoff;
  70     if(fabs(num_term) < cutoff) num_term = cutoff;
  71     den_term  = 1.0/den_term;
  72
  73     delta_frac = den_term * num_term;
  74     cf *= delta_frac;
  75
  76     coeff = -(a+k)*(a+b+k)*x/((a+2*k)*(a+2*k+1.0));
  77
  78     /* second step */
  79     den_term = 1.0 + coeff*den_term;
  80     num_term = 1.0 + coeff/num_term;
  81     if(fabs(den_term) < cutoff) den_term = cutoff;
  82     if(fabs(num_term) < cutoff) num_term = cutoff;
  83     den_term = 1.0/den_term;
  84
  85     delta_frac = den_term*num_term;
  86     cf *= delta_frac;
  87
  88     if(fabs(delta_frac-1.0) < 2.0*GSL_DBL_EPSILON) break;
  89
  90     ++iter_count;
  91   }
  92
  93   result->val = cf;
  94   result->err = iter_count * 4.0 * GSL_DBL_EPSILON * fabs(cf);
  95
  96   if(iter_count >= max_iter)
  97     GSL_ERROR ("error", GSL_EMAXITER);
  98   else
  99     return GSL_SUCCESS;
 100 }
 101
 102
 103
 104 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 105
 106 int
 107 gsl_sf_beta_inc_e(
 108   const double a,
 109   const double b,
 110   const double x,
 111   gsl_sf_result * result
 112   )
 113 {
 114   if(x < 0.0 || x > 1.0) {
 115     DOMAIN_ERROR(result);
 116   } else if (isnegint(a) || isnegint(b)) {
 117     DOMAIN_ERROR(result);
 118   } else if (isnegint(a+b)) {
 119     DOMAIN_ERROR(result);
 120   } else if(x == 0.0) {
 121     result->val = 0.0;
 122     result->err = 0.0;
 123     return GSL_SUCCESS;
 124   }
 125   else if(x == 1.0) {
 126     result->val = 1.0;
 127     result->err = 0.0;
 128     return GSL_SUCCESS;
 129   } else if (a <= 0 || b <= 0) {
 130     gsl_sf_result f, beta;
 131     int stat;
 132     const int stat_f = gsl_sf_hyperg_2F1_e(a, 1-b, a+1, x, &f);
 133     const int stat_beta = gsl_sf_beta_e(a, b, &beta);
 134     double prefactor = (pow(x, a) / a);
 135     result->val = prefactor * f.val / beta.val;
 136     result->err = fabs(prefactor) * f.err/ fabs(beta.val) + fabs(result->val/beta.val) * beta.err;
 137
 138     stat = GSL_ERROR_SELECT_2(stat_f, stat_beta);
 139     if(stat == GSL_SUCCESS) {
 140       CHECK_UNDERFLOW(result);
 141     }
 142     return stat;
 143   } else {
 144     gsl_sf_result ln_beta;
 145     gsl_sf_result ln_x;
 146     gsl_sf_result ln_1mx;
 147     gsl_sf_result prefactor;
 148     const int stat_ln_beta = gsl_sf_lnbeta_e(a, b, &ln_beta);
 149     const int stat_ln_1mx = gsl_sf_log_1plusx_e(-x, &ln_1mx);
 150     const int stat_ln_x = gsl_sf_log_e(x, &ln_x);
 151     const int stat_ln = GSL_ERROR_SELECT_3(stat_ln_beta, stat_ln_1mx, stat_ln_x);
 152
 153     const double ln_pre_val = -ln_beta.val + a * ln_x.val + b * ln_1mx.val;
 154     const double ln_pre_err =  ln_beta.err + fabs(a*ln_x.err) + fabs(b*ln_1mx.err);
 155     const int stat_exp = gsl_sf_exp_err_e(ln_pre_val, ln_pre_err, &prefactor);
 156
 157     if(stat_ln != GSL_SUCCESS) {
 158       result->val = 0.0;
 159       result->err = 0.0;
 160       GSL_ERROR ("error", GSL_ESANITY);
 161     }
 162
 163     if(x < (a + 1.0)/(a+b+2.0)) {
 164       /* Apply continued fraction directly. */
 165       gsl_sf_result cf;
 166       const int stat_cf = beta_cont_frac(a, b, x, &cf);
 167       int stat;
 168       result->val = prefactor.val * cf.val / a;
 169       result->err = (fabs(prefactor.err * cf.val) + fabs(prefactor.val * cf.err))/a;
 170
 171       stat = GSL_ERROR_SELECT_2(stat_exp, stat_cf);
 172       if(stat == GSL_SUCCESS) {
 173         CHECK_UNDERFLOW(result);
 174       }
 175       return stat;
 176     }
 177     else {
 178       /* Apply continued fraction after hypergeometric transformation. */
 179       gsl_sf_result cf;
 180       const int stat_cf = beta_cont_frac(b, a, 1.0-x, &cf);
 181       int stat;
 182       const double term = prefactor.val * cf.val / b;
 183       result->val  = 1.0 - term;
 184       result->err  = fabs(prefactor.err * cf.val)/b;
 185       result->err += fabs(prefactor.val * cf.err)/b;
 186       result->err += 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(term));
 187       stat = GSL_ERROR_SELECT_2(stat_exp, stat_cf);
 188       if(stat == GSL_SUCCESS) {
 189         CHECK_UNDERFLOW(result);
 190       }
 191       return stat;
 192     }
 193   }
 194 }
 195
 196
 197 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 198
 199 #include "eval.h"
 200
 201 double gsl_sf_beta_inc(const double a, const double b, const double x)
 202 {
 203   EVAL_RESULT(gsl_sf_beta_inc_e(a, b, x, &result));
 204 }