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

   1 /* specfunc/beta.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
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_log.h>
  27 #include <gsl/gsl_sf_psi.h>
  28 #include <gsl/gsl_sf_gamma.h>
  29
  30 #include "error.h"
  31
  32 static double
  33 isnegint (const double x)
  34 {
  35   return (x < 0) && (x == floor(x));
  36 }
  37
  38 int
  39 gsl_sf_lnbeta_e(const double x, const double y, gsl_sf_result * result)
  40 {
  41   double sgn;
  42   int status = gsl_sf_lnbeta_sgn_e(x,y,result,&sgn);
  43   if (sgn == -1) {
  44     DOMAIN_ERROR(result);
  45   }
  46   return status;
  47 }
  48
  49 int
  50 gsl_sf_lnbeta_sgn_e(const double x, const double y, gsl_sf_result * result, double * sgn)
  51 {
  52   /* CHECK_POINTER(result) */
  53
  54   if(x == 0.0 || y == 0.0) {
  55     *sgn = 0.0;
  56     DOMAIN_ERROR(result);
  57   } else if (isnegint(x) || isnegint(y)) {
  58     *sgn = 0.0;
  59     DOMAIN_ERROR(result); /* not defined for negative integers */
  60   }
  61
  62   /* See if we can handle the postive case with min/max < 0.2 */
  63
  64   if (x > 0 && y > 0) {
  65     const double max = GSL_MAX(x,y);
  66     const double min = GSL_MIN(x,y);
  67     const double rat = min/max;
  68
  69     if(rat < 0.2) {
  70       /* min << max, so be careful
  71        * with the subtraction
  72        */
  73       double lnpre_val;
  74       double lnpre_err;
  75       double lnpow_val;
  76       double lnpow_err;
  77       double t1, t2, t3;
  78       gsl_sf_result lnopr;
  79       gsl_sf_result gsx, gsy, gsxy;
  80       gsl_sf_gammastar_e(x, &gsx);
  81       gsl_sf_gammastar_e(y, &gsy);
  82       gsl_sf_gammastar_e(x+y, &gsxy);
  83       gsl_sf_log_1plusx_e(rat, &lnopr);
  84       lnpre_val = log(gsx.val*gsy.val/gsxy.val * M_SQRT2*M_SQRTPI);
  85       lnpre_err = gsx.err/gsx.val + gsy.err/gsy.val + gsxy.err/gsxy.val;
  86       t1 = min*log(rat);
  87       t2 = 0.5*log(min);
  88       t3 = (x+y-0.5)*lnopr.val;
  89       lnpow_val  = t1 - t2 - t3;
  90       lnpow_err  = GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
  91       lnpow_err += fabs(x+y-0.5) * lnopr.err;
  92       result->val  = lnpre_val + lnpow_val;
  93       result->err  = lnpre_err + lnpow_err;
  94       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  95       *sgn = 1.0;
  96       return GSL_SUCCESS;
  97     }
  98   }
  99
 100   /* General case - Fallback */
 101   {
 102     gsl_sf_result lgx, lgy, lgxy;
 103     double sgx, sgy, sgxy, xy = x+y;
 104     int stat_gx  = gsl_sf_lngamma_sgn_e(x, &lgx, &sgx);
 105     int stat_gy  = gsl_sf_lngamma_sgn_e(y, &lgy, &sgy);
 106     int stat_gxy = gsl_sf_lngamma_sgn_e(xy, &lgxy, &sgxy);
 107     *sgn = sgx * sgy * sgxy;
 108     result->val  = lgx.val + lgy.val - lgxy.val;
 109     result->err  = lgx.err + lgy.err + lgxy.err;
 110     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(lgx.val) + fabs(lgy.val) + fabs(lgxy.val));
 111     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 112     return GSL_ERROR_SELECT_3(stat_gx, stat_gy, stat_gxy);
 113   }
 114 }
 115
 116
 117 int
 118 gsl_sf_beta_e(const double x, const double y, gsl_sf_result * result)
 119 {
 120   if((x > 0 && y > 0) && x < 50.0 && y < 50.0) {
 121     /* Handle the easy case */
 122     gsl_sf_result gx, gy, gxy;
 123     gsl_sf_gamma_e(x, &gx);
 124     gsl_sf_gamma_e(y, &gy);
 125     gsl_sf_gamma_e(x+y, &gxy);
 126     result->val  = (gx.val*gy.val)/gxy.val;
 127     result->err  = gx.err * fabs(gy.val/gxy.val);
 128     result->err += gy.err * fabs(gx.val/gxy.val);
 129     result->err += fabs((gx.val*gy.val)/(gxy.val*gxy.val)) * gxy.err;
 130     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 131     return GSL_SUCCESS;
 132   }
 133   else if (isnegint(x) || isnegint(y)) {
 134     DOMAIN_ERROR(result);
 135   } else if (isnegint(x+y)) {  /* infinity in the denominator */
 136     result->val = 0.0;
 137     result->err = 0.0;
 138     return GSL_SUCCESS;
 139   } else {
 140     gsl_sf_result lb;
 141     double sgn;
 142     int stat_lb = gsl_sf_lnbeta_sgn_e(x, y, &lb, &sgn);
 143     if(stat_lb == GSL_SUCCESS) {
 144       int status = gsl_sf_exp_err_e(lb.val, lb.err, result);
 145       result->val *= sgn;
 146       return status;
 147     }
 148     else {
 149       result->val = 0.0;
 150       result->err = 0.0;
 151       return stat_lb;
 152     }
 153   }
 154 }
 155
 156
 157 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 158
 159 #include "eval.h"
 160
 161 double gsl_sf_lnbeta(const double x, const double y)
 162 {
 163   EVAL_RESULT(gsl_sf_lnbeta_e(x, y, &result));
 164 }
 165
 166 double gsl_sf_beta(const double x, const double y)
 167 {
 168   EVAL_RESULT(gsl_sf_beta_e(x, y, &result));
 169 }