beta.c
Go to the documentation of this file.
1 /* beta.c
2  *
3  * Beta function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double a, b, y, beta();
10  *
11  * y = beta( a, b );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * - -
18  * | (a) | (b)
19  * beta( a, b ) = -----------.
20  * -
21  * | (a+b)
22  *
23  * For large arguments the logarithm of the function is
24  * evaluated using lgam(), then exponentiated.
25  *
26  *
27  *
28  * ACCURACY:
29  *
30  * Relative error:
31  * arithmetic domain # trials peak rms
32  * IEEE 0,30 30000 8.1e-14 1.1e-14
33  *
34  * ERROR MESSAGES:
35  *
36  * message condition value returned
37  * beta overflow log(beta) > MAXLOG 0.0
38  * a or b <0 integer 0.0
39  *
40  */
41 
42 
43 /*
44  * Cephes Math Library Release 2.0: April, 1987
45  * Copyright 1984, 1987 by Stephen L. Moshier
46  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
47  */
48 
49 #include "mconf.h"
50 
51 #define MAXGAM 171.624376956302725
52 
53 extern double MAXLOG;
54 
55 #define ASYMP_FACTOR 1e6
56 
57 static double lbeta_asymp(double a, double b, int *sgn);
58 static double lbeta_negint(int a, double b);
59 static double beta_negint(int a, double b);
60 
61 double beta(double a, double b)
62 {
63  double y;
64  int sign = 1;
65 
66  if (a <= 0.0) {
67  if (a == floor(a)) {
68  if (a == (int)a) {
69  return beta_negint((int)a, b);
70  }
71  else {
72  goto overflow;
73  }
74  }
75  }
76 
77  if (b <= 0.0) {
78  if (b == floor(b)) {
79  if (b == (int)b) {
80  return beta_negint((int)b, a);
81  }
82  else {
83  goto overflow;
84  }
85  }
86  }
87 
88  if (fabs(a) < fabs(b)) {
89  y = a; a = b; b = y;
90  }
91 
92  if (fabs(a) > ASYMP_FACTOR * fabs(b) && a > ASYMP_FACTOR) {
93  /* Avoid loss of precision in lgam(a + b) - lgam(a) */
94  y = lbeta_asymp(a, b, &sign);
95  return sign * exp(y);
96  }
97 
98  y = a + b;
99  if (fabs(y) > MAXGAM || fabs(a) > MAXGAM || fabs(b) > MAXGAM) {
100  int sgngam;
101  y = lgam_sgn(y, &sgngam);
102  sign *= sgngam; /* keep track of the sign */
103  y = lgam_sgn(b, &sgngam) - y;
104  sign *= sgngam;
105  y = lgam_sgn(a, &sgngam) + y;
106  sign *= sgngam;
107  if (y > MAXLOG) {
108  goto overflow;
109  }
110  return (sign * exp(y));
111  }
112 
113  y = Gamma(y);
114  a = Gamma(a);
115  b = Gamma(b);
116  if (y == 0.0)
117  goto overflow;
118 
119  if (fabs(fabs(a) - fabs(y)) > fabs(fabs(b) - fabs(y))) {
120  y = b / y;
121  y *= a;
122  }
123  else {
124  y = a / y;
125  y *= b;
126  }
127 
128  return (y);
129 
130 overflow:
131  sf_error("beta", SF_ERROR_OVERFLOW, NULL);
132  return (sign * INFINITY);
133 }
134 
135 
136 /* Natural log of |beta|. */
137 
138 double lbeta(double a, double b)
139 {
140  double y;
141  int sign;
142 
143  sign = 1;
144 
145  if (a <= 0.0) {
146  if (a == floor(a)) {
147  if (a == (int)a) {
148  return lbeta_negint((int)a, b);
149  }
150  else {
151  goto over;
152  }
153  }
154  }
155 
156  if (b <= 0.0) {
157  if (b == floor(b)) {
158  if (b == (int)b) {
159  return lbeta_negint((int)b, a);
160  }
161  else {
162  goto over;
163  }
164  }
165  }
166 
167  if (fabs(a) < fabs(b)) {
168  y = a; a = b; b = y;
169  }
170 
171  if (fabs(a) > ASYMP_FACTOR * fabs(b) && a > ASYMP_FACTOR) {
172  /* Avoid loss of precision in lgam(a + b) - lgam(a) */
173  y = lbeta_asymp(a, b, &sign);
174  return y;
175  }
176 
177  y = a + b;
178  if (fabs(y) > MAXGAM || fabs(a) > MAXGAM || fabs(b) > MAXGAM) {
179  int sgngam;
180  y = lgam_sgn(y, &sgngam);
181  sign *= sgngam; /* keep track of the sign */
182  y = lgam_sgn(b, &sgngam) - y;
183  sign *= sgngam;
184  y = lgam_sgn(a, &sgngam) + y;
185  sign *= sgngam;
186  return (y);
187  }
188 
189  y = Gamma(y);
190  a = Gamma(a);
191  b = Gamma(b);
192  if (y == 0.0) {
193  over:
194  sf_error("lbeta", SF_ERROR_OVERFLOW, NULL);
195  return (sign * INFINITY);
196  }
197 
198  if (fabs(fabs(a) - fabs(y)) > fabs(fabs(b) - fabs(y))) {
199  y = b / y;
200  y *= a;
201  }
202  else {
203  y = a / y;
204  y *= b;
205  }
206 
207  if (y < 0) {
208  y = -y;
209  }
210 
211  return (log(y));
212 }
213 
214 /*
215  * Asymptotic expansion for ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
216  */
217 static double lbeta_asymp(double a, double b, int *sgn)
218 {
219  double r = lgam_sgn(b, sgn);
220  r -= b * log(a);
221 
222  r += b*(1-b)/(2*a);
223  r += b*(1-b)*(1-2*b)/(12*a*a);
224  r += - b*b*(1-b)*(1-b)/(12*a*a*a);
225 
226  return r;
227 }
228 
229 
230 /*
231  * Special case for a negative integer argument
232  */
233 
234 static double beta_negint(int a, double b)
235 {
236  int sgn;
237  if (b == (int)b && 1 - a - b > 0) {
238  sgn = ((int)b % 2 == 0) ? 1 : -1;
239  return sgn * beta(1 - a - b, b);
240  }
241  else {
242  sf_error("lbeta", SF_ERROR_OVERFLOW, NULL);
243  return INFINITY;
244  }
245 }
246 
247 static double lbeta_negint(int a, double b)
248 {
249  double r;
250  if (b == (int)b && 1 - a - b > 0) {
251  r = lbeta(1 - a - b, b);
252  return r;
253  }
254  else {
255  sf_error("lbeta", SF_ERROR_OVERFLOW, NULL);
256  return INFINITY;
257  }
258 }
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
beta_negint
static double beta_negint(int a, double b)
Definition: beta.c:234
b
Scalar * b
Definition: benchVecAdd.cpp:17
Gamma
double Gamma(double x)
Definition: gamma.c:160
sign
const EIGEN_DEVICE_FUNC SignReturnType sign() const
Definition: ArrayCwiseUnaryOps.h:219
SF_ERROR_OVERFLOW
@ SF_ERROR_OVERFLOW
Definition: sf_error.h:12
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
beta
double beta(double a, double b)
Definition: beta.c:61
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
lgam_sgn
double lgam_sgn(double x, int *sign)
Definition: gamma.c:281
MAXGAM
#define MAXGAM
Definition: beta.c:51
y
Scalar * y
Definition: level1_cplx_impl.h:124
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
lbeta
double lbeta(double a, double b)
Definition: beta.c:138
MAXLOG
double MAXLOG
Definition: const.c:57
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
lbeta_asymp
static double lbeta_asymp(double a, double b, int *sgn)
Definition: beta.c:217
NULL
#define NULL
Definition: ccolamd.c:609
ASYMP_FACTOR
#define ASYMP_FACTOR
Definition: beta.c:55
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
lbeta_negint
static double lbeta_negint(int a, double b)
Definition: beta.c:247


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:01