The R Project SVN R

Rev

Rev 42302 | Rev 42549 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 42302 Rev 42303
Line -... Line 1...
- 
 
 1 
/*
- 
 
 2 
 *  Mathlib : A C Library of Special Functions
- 
 
 3 
 *  Copyright (C) 1998 Ross Ihaka
- 
 
 4 
 *  Copyright (C) 2000-2001 The R Development Core Team
- 
 
 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 2 of the License, or
- 
 
 9 
 *  (at your option) any later version.
- 
 
 10 
 *
- 
 
 11 
 *  This program is distributed in the hope that it will be useful,
- 
 
 12 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- 
 
 13 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- 
 
 14 
 *  GNU 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, a copy is available at
- 
 
 18 
 *  http://www.r-project.org/Licenses/
- 
 
 19 
 *
- 
 
 20 
 *  SYNOPSIS
- 
 
 21 
 *
- 
 
 22 
 *    #include <Rmath.h>
- 
 
 23 
 *    double beta(double a, double b);
- 
 
 24 
 *
- 
 
 25 
 *  DESCRIPTION
- 
 
 26 
 *
- 
 
 27 
 *    This function returns the value of the beta function
- 
 
 28 
 *    evaluated with arguments a and b.
- 
 
 29 
 *
- 
 
 30 
 *  NOTES
- 
 
 31 
 *
- 
 
 32 
 *    This routine is a translation into C of a Fortran subroutine
- 
 
 33 
 *    by W. Fullerton of Los Alamos Scientific Laboratory.
- 
 
 34 
 *    Some modifications have been made so that the routines
- 
 
 35 
 *    conform to the IEEE 754 standard.
- 
 
 36 
 */
- 
 
 37 
 
- 
 
 38 
#include "nmath.h"
- 
 
 39 
 
- 
 
 40 
double beta(double a, double b)
- 
 
 41 
{
- 
 
 42 
    double val;
- 
 
 43 
 
- 
 
 44 
#ifdef NOMORE_FOR_THREADS
- 
 
 45 
    static double xmin, xmax = 0;/*-> typically = 171.61447887 for IEEE */
- 
 
 46 
    static double lnsml = 0;/*-> typically = -708.3964185 */
- 
 
 47 
 
- 
 
 48 
    if (xmax == 0) {
- 
 
 49 
	    gammalims(&xmin, &xmax);
- 
 
 50 
	    lnsml = log(d1mach(1));
- 
 
 51 
    }
- 
 
 52 
#else
- 
 
 53 
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
- 
 
 54 
 *   xmin, xmax : see ./gammalims.c
- 
 
 55 
 *   lnsml = log(DBL_MIN) = log(2 ^ -1022) = -1022 * log(2)
- 
 
 56 
*/
- 
 
 57 
# define xmin  -170.5674972726612
- 
 
 58 
# define xmax   171.61447887182298
- 
 
 59 
# define lnsml -708.39641853226412
- 
 
 60 
#endif
- 
 
 61 
 
- 
 
 62 
 
- 
 
 63 
#ifdef IEEE_754
- 
 
 64 
    /* NaNs propagated correctly */
- 
 
 65 
    if(ISNAN(a) || ISNAN(b)) return a + b;
- 
 
 66 
#endif
- 
 
 67 
 
- 
 
 68 
    if (a < 0 || b < 0)
- 
 
 69 
	ML_ERR_return_NAN
- 
 
 70 
    else if (a == 0 || b == 0) {
- 
 
 71 
	return ML_POSINF;
- 
 
 72 
    }
- 
 
 73 
    else if (!R_FINITE(a) || !R_FINITE(b)) {
- 
 
 74 
	return 0;
- 
 
 75 
    }
- 
 
 76 
 
- 
 
 77 
    if (a + b < xmax) /* ~= 171.61 for IEEE */
- 
 
 78 
	return gammafn(a) * gammafn(b) / gammafn(a+b);
- 
 
 79 
 
- 
 
 80 
    val = lbeta(a, b);
- 
 
 81 
    if (val < lnsml) {
- 
 
 82 
	/* a and/or b so big that beta underflows */
- 
 
 83 
	ML_ERROR(ME_UNDERFLOW, "beta");
- 
 
 84 
	return ML_UNDERFLOW;
- 
 
 85 
    }
- 
 
 86 
    return exp(val);
- 
 
 87 
}