The R Project SVN R

Rev

Rev 13433 | Rev 37461 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
574 ihaka 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
13433 maechler 4 
 *  Copyright (C) 2000-2001 The R Development Core Team
574 ihaka 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, write to the Free Software
5458 ripley 18 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
574 ihaka 19 
 *
20 
 *  SYNOPSIS
21 
 *
19500 hornik 22 
 *    #include <Rmath.h>
574 ihaka 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 
 
8431 ripley 38 
#include "nmath.h"
574 ihaka 39 
 
40 
double beta(double a, double b)
41 
{
13433 maechler 42 
    double val;
574 ihaka 43 
 
13433 maechler 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 
 
574 ihaka 48 
    if (xmax == 0) {
49 
	    gammalims(&xmin, &xmax);
13433 maechler 50 
	    lnsml = log(d1mach(1));
574 ihaka 51 
    }
13433 maechler 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
574 ihaka 61 
 
13433 maechler 62 
 
580 ihaka 63 
#ifdef IEEE_754
64 
    /* NaNs propagated correctly */
65 
    if(ISNAN(a) || ISNAN(b)) return a + b;
66 
#endif
67 
 
7671 maechler 68 
    if (a < 0 || b < 0)
69 
	ML_ERR_return_NAN
928 maechler 70 
    else if (a == 0 || b == 0) {
71 
	return ML_POSINF;
72 
    }
5108 maechler 73 
    else if (!R_FINITE(a) || !R_FINITE(b)) {
928 maechler 74 
	return 0;
75 
    }
574 ihaka 76 
 
10418 maechler 77 
    if (a + b < xmax) /* ~= 171.61 for IEEE */
2276 maechler 78 
	return gammafn(a) * gammafn(b) / gammafn(a+b);
574 ihaka 79 
 
80 
    val = lbeta(a, b);
13433 maechler 81 
    if (val < lnsml) {
574 ihaka 82 
	/* a and/or b so big that beta underflows */
928 maechler 83 
	ML_ERROR(ME_UNDERFLOW);
574 ihaka 84 
	return ML_UNDERFLOW;
85 
    }
86 
    return exp(val);
87 
}