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

Diff of /gsl/cdf/beta.c

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

revision 1.1.2.2 by bjg, Mon Jun 23 17:18:31 2003 UTC revision 1.1.2.3 by bjg, Wed Jun 25 15:46:15 2003 UTC
# Line 1  Line 1 
1  /* specfunc/beta.c  /* cdf/cdf_beta.c
2   *   *
3   * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman   * Copyright (C) 2003 Brian Gough.
4   *   *
5   * This program is free software; you can redistribute it and/or modify   * 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   * it under the terms of the GNU General Public License as published by
# Line 14  Line 14 
14   *   *
15   * You should have received a copy of the GNU General Public License   * You should have received a copy of the GNU General Public License
16   * along with this program; if not, write to the Free Software   * along with this program; if not, write to the Free Software
17   * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.   * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  */  
   
 /* Author:  G. Jungman  
  * Tweaked by Jason Stover for the cdf module.  
18   */   */
19    
20  #include <config.h>  #include <config.h>
21    #include <math.h>
22    #include <gsl/gsl_cdf.h>
23    #include <gsl/gsl_sf_gamma.h>
24  #include <gsl/gsl_math.h>  #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>  
25    
26  #include "exp.h"  #include "beta_inc.c"
 #include "log.h"  
 #include "psi.h"  
 #include "gamma.h"  
 #include "error.h"  
27    
28  int  double
29  gsl_cdf_lnbeta_e(const double x, const double y, gsl_cdf_result * result)  gsl_cdf_beta_P (double x, double a, double b)
30  {  {
31    /* CHECK_POINTER(result) */    double P;
32    
33    if(x <= 0.0 || y <= 0.0) {    if (x <= 0.0 || x > 1.0)
34      DOMAIN_ERROR(result);      {
35    }        return 0.0;
   else {  
     const double max = GSL_MAX(x,y);  
     const double min = GSL_MIN(x,y);  
     const double rat = min/max;  
   
     if(rat < 0.2) {  
       /* min << max, so be careful  
        * with the subtraction  
        */  
       double lnpre_val;  
       double lnpre_err;  
       double lnpow_val;  
       double lnpow_err;  
       double t1, t2, t3;  
       gsl_cdf_result lnopr;  
       gsl_cdf_result gsx, gsy, gsxy;  
       gsl_cdf_gammastar_e(x, &gsx);  
       gsl_cdf_gammastar_e(y, &gsy);  
       gsl_cdf_gammastar_e(x+y, &gsxy);  
       gsl_cdf_log_1plusx_e(rat, &lnopr);  
       lnpre_val = log(gsx.val*gsy.val/gsxy.val * M_SQRT2*M_SQRTPI);  
       lnpre_err = gsx.err/gsx.val + gsy.err/gsy.val + gsxy.err/gsxy.val;  
       t1 = min*log(rat);  
       t2 = 0.5*log(min);  
       t3 = (x+y-0.5)*lnopr.val;  
       lnpow_val  = t1 - t2 - t3;  
       lnpow_err  = GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));  
       lnpow_err += fabs(x+y-0.5) * lnopr.err;  
       result->val  = lnpre_val + lnpow_val;  
       result->err  = lnpre_err + lnpow_err;  
       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);  
       return GSL_SUCCESS;  
36      }      }
     else {  
       gsl_cdf_result lgx, lgy, lgxy;  
       int stat_gx  = gsl_cdf_lngamma_e(x, &lgx);  
       int stat_gy  = gsl_cdf_lngamma_e(y, &lgy);  
       int stat_gxy = gsl_cdf_lngamma_e(x+y, &lgxy);  
       result->val  = lgx.val + lgy.val - lgxy.val;  
       result->err  = lgx.err + lgy.err + lgxy.err;  
       result->err += GSL_DBL_EPSILON * (fabs(lgx.val) + fabs(lgy.val) + fabs(lgxy.val));  
       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);  
       return GSL_ERROR_SELECT_3(stat_gx, stat_gy, stat_gxy);  
     }  
   }  
 }  
37    
38      P = beta_inc_AXPY (1.0, 0.0, a, b, x);
39    
40  int    return P;
 gsl_cdf_beta_e(const double x, const double y, gsl_cdf_result * result)  
 {  
   if(x < 50.0 && y < 50.0) {  
     gsl_cdf_result gx, gy, gxy;  
     gsl_cdf_gamma_e(x, &gx);  
     gsl_cdf_gamma_e(y, &gy);  
     gsl_cdf_gamma_e(x+y, &gxy);  
     result->val  = (gx.val*gy.val)/gxy.val;  
     result->err  = gx.err * gy.val/gxy.val;  
     result->err += gy.err * gx.val/gxy.val;  
     result->err += (gx.val*gy.val)/(gxy.val*gxy.val) * gxy.err;  
     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);  
     return GSL_SUCCESS;  
   }  
   else {  
     gsl_cdf_result lb;  
     int stat_lb = gsl_cdf_lnbeta_e(x, y, &lb);  
     if(stat_lb == GSL_SUCCESS) {  
       return gsl_cdf_exp_err_e(lb.val, lb.err, result);  
     }  
     else {  
       result->val = 0.0;  
       result->err = 0.0;  
       return stat_lb;  
     }  
   }  
41  }  }
42    
43    double
44    gsl_cdf_beta_Q (double x, double a, double b)
45    {
46      double P;
47    
48  /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/    if (x < 0.0 || x >= 1.0)
49        {
50  #include "eval.h"        return 0.0;
51        }
52    
53  double gsl_cdf_lnbeta(const double x, const double y)    P = beta_inc_AXPY (-1.0, 1.0, a, b, x);
 {  
   EVAL_RESULT(gsl_cdf_lnbeta_e(x, y, &result));  
 }  
54    
55  double gsl_cdf_beta(const double x, const double y)    return P;
 {  
   EVAL_RESULT(gsl_cdf_beta_e(x, y, &result));  
56  }  }

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