/src/r-source/src/nmath/beta.c
Line | Count | Source |
1 | | /* |
2 | | * Mathlib : A C Library of Special Functions |
3 | | * Copyright (C) 1998 Ross Ihaka |
4 | | * Copyright (C) 2000-2014 The R 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 | | * https://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 | 0 | { |
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 | 0 | # define xmin -170.5674972726612 |
56 | 0 | # define xmax 171.61447887182298 |
57 | 0 | # define lnsml -708.39641853226412 |
58 | 0 | #endif |
59 | | |
60 | |
|
61 | 0 | #ifdef IEEE_754 |
62 | | /* NaNs propagated correctly */ |
63 | 0 | if(ISNAN(a) || ISNAN(b)) return a + b; |
64 | 0 | #endif |
65 | | |
66 | 0 | if (a < 0 || b < 0) |
67 | 0 | ML_WARN_return_NAN |
68 | 0 | else if (a == 0 || b == 0) |
69 | 0 | return ML_POSINF; |
70 | 0 | else if (!R_FINITE(a) || !R_FINITE(b)) |
71 | 0 | return 0; |
72 | | |
73 | 0 | 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. |
77 | | gammafn(x) can still overflow for x ~ 1e-308, |
78 | | but the result would too. |
79 | | */ |
80 | 0 | return (1 / gammafn(a+b)) * gammafn(a) * gammafn(b); |
81 | 0 | } else { |
82 | 0 | double val = lbeta(a, b); |
83 | | // underflow to 0 is not harmful per se; exp(-999) also gives no warning |
84 | | #ifndef IEEE_754 |
85 | | if (val < lnsml) { |
86 | | /* a and/or b so big that beta underflows */ |
87 | | ML_WARNING(ME_UNDERFLOW, "beta"); |
88 | | /* return ML_UNDERFLOW; pointless giving incorrect value */ |
89 | | } |
90 | | #endif |
91 | 0 | return exp(val); |
92 | 0 | } |
93 | 0 | } |