| 42303 |
ripley |
1 |
/*
|
|
|
2 |
* Mathlib : A C Library of Special Functions
|
|
|
3 |
* Copyright (C) 1998 Ross Ihaka
|
| 60970 |
ripley |
4 |
* Copyright (C) 2000-2012 The R Core Team
|
| 42303 |
ripley |
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 |
#ifdef NOMORE_FOR_THREADS
|
|
|
43 |
static double xmin, xmax = 0;/*-> typically = 171.61447887 for IEEE */
|
|
|
44 |
static double lnsml = 0;/*-> typically = -708.3964185 */
|
|
|
45 |
|
|
|
46 |
if (xmax == 0) {
|
|
|
47 |
gammalims(&xmin, &xmax);
|
|
|
48 |
lnsml = log(d1mach(1));
|
|
|
49 |
}
|
|
|
50 |
#else
|
|
|
51 |
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
|
|
|
52 |
* xmin, xmax : see ./gammalims.c
|
|
|
53 |
* lnsml = log(DBL_MIN) = log(2 ^ -1022) = -1022 * log(2)
|
|
|
54 |
*/
|
|
|
55 |
# define xmin -170.5674972726612
|
|
|
56 |
# define xmax 171.61447887182298
|
|
|
57 |
# define lnsml -708.39641853226412
|
|
|
58 |
#endif
|
|
|
59 |
|
|
|
60 |
|
|
|
61 |
#ifdef IEEE_754
|
|
|
62 |
/* NaNs propagated correctly */
|
|
|
63 |
if(ISNAN(a) || ISNAN(b)) return a + b;
|
|
|
64 |
#endif
|
|
|
65 |
|
|
|
66 |
if (a < 0 || b < 0)
|
|
|
67 |
ML_ERR_return_NAN
|
| 60970 |
ripley |
68 |
else if (a == 0 || b == 0)
|
| 42303 |
ripley |
69 |
return ML_POSINF;
|
| 60970 |
ripley |
70 |
else if (!R_FINITE(a) || !R_FINITE(b))
|
| 42303 |
ripley |
71 |
return 0;
|
|
|
72 |
|
| 60970 |
ripley |
73 |
if (a + b < xmax) {/* ~= 171.61 for IEEE */
|
|
|
74 |
// return gammafn(a) * gammafn(b) / gammafn(a+b);
|
|
|
75 |
/* All the terms are positive, and all can be large for large
|
|
|
76 |
or small arguments. They are never much less than one.
|
| 60971 |
ripley |
77 |
gammafn(x) can still overflow for x ~ 1e-308,
|
|
|
78 |
but the result would too.
|
| 60970 |
ripley |
79 |
*/
|
|
|
80 |
return (1 / gammafn(a+b)) * gammafn(a) * gammafn(b);
|
|
|
81 |
} else {
|
| 50406 |
maechler |
82 |
double val = lbeta(a, b);
|
|
|
83 |
if (val < lnsml) {
|
|
|
84 |
/* a and/or b so big that beta underflows */
|
|
|
85 |
ML_ERROR(ME_UNDERFLOW, "beta");
|
|
|
86 |
/* return ML_UNDERFLOW; pointless giving incorrect value */
|
|
|
87 |
}
|
|
|
88 |
return exp(val);
|
| 42303 |
ripley |
89 |
}
|
|
|
90 |
}
|
The R Project SVN
R