1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
/* betaf.c
*
* Beta function
*
*
*
* SYNOPSIS:
*
* float a, b, y, betaf();
*
* y = betaf( a, b );
*
*
*
* DESCRIPTION:
*
* - -
* | (a) | (b)
* beta( a, b ) = -----------.
* -
* | (a+b)
*
* For large arguments the logarithm of the function is
* evaluated using lgam(), then exponentiated.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0,30 10000 4.0e-5 6.0e-6
* IEEE -20,0 10000 4.9e-3 5.4e-5
*
* ERROR MESSAGES:
*
* message condition value returned
* betaf overflow log(beta) > MAXLOG 0.0
* a or b <0 integer 0.0
*
*/
/* beta.c */
/*
Cephes Math Library Release 2.2: July, 1992
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <math.h>
#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
#define MAXGAM 34.84425627277176174
extern float MAXLOGF, MAXNUMF;
extern int sgngamf;
#ifdef ANSIC
float gammaf(float), lgamf(float), expf(float), floorf(float);
#else
float gammaf(), lgamf(), expf(), floorf();
#endif
float betaf( float aa, float bb )
{
float a, b, y;
int sign;
sign = 1;
a = aa;
b = bb;
if( a <= 0.0 )
{
if( a == floorf(a) )
goto over;
}
if( b <= 0.0 )
{
if( b == floorf(b) )
goto over;
}
y = a + b;
if( fabsf(y) > MAXGAM )
{
y = lgamf(y);
sign *= sgngamf; /* keep track of the sign */
y = lgamf(b) - y;
sign *= sgngamf;
y = lgamf(a) + y;
sign *= sgngamf;
if( y > MAXLOGF )
{
over:
mtherr( "betaf", OVERFLOW );
return( sign * MAXNUMF );
}
return( sign * expf(y) );
}
y = gammaf(y);
if( y == 0.0 )
goto over;
if( a > b )
{
y = gammaf(a)/y;
y *= gammaf(b);
}
else
{
y = gammaf(b)/y;
y *= gammaf(a);
}
return(y);
}
|
