diff options
Diffstat (limited to 'libm/double/beta.c')
| -rw-r--r-- | libm/double/beta.c | 201 |
1 files changed, 201 insertions, 0 deletions
diff --git a/libm/double/beta.c b/libm/double/beta.c new file mode 100644 index 000000000..410760f32 --- /dev/null +++ b/libm/double/beta.c @@ -0,0 +1,201 @@ +/* beta.c + * + * Beta function + * + * + * + * SYNOPSIS: + * + * double a, b, y, beta(); + * + * y = beta( 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 + * DEC 0,30 1700 7.7e-15 1.5e-15 + * IEEE 0,30 30000 8.1e-14 1.1e-14 + * + * ERROR MESSAGES: + * + * message condition value returned + * beta overflow log(beta) > MAXLOG 0.0 + * a or b <0 integer 0.0 + * + */ + +/* beta.c */ + + +/* +Cephes Math Library Release 2.0: April, 1987 +Copyright 1984, 1987 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include <math.h> + +#ifdef UNK +#define MAXGAM 34.84425627277176174 +#endif +#ifdef DEC +#define MAXGAM 34.84425627277176174 +#endif +#ifdef IBMPC +#define MAXGAM 171.624376956302725 +#endif +#ifdef MIEEE +#define MAXGAM 171.624376956302725 +#endif + +#ifdef ANSIPROT +extern double fabs ( double ); +extern double gamma ( double ); +extern double lgam ( double ); +extern double exp ( double ); +extern double log ( double ); +extern double floor ( double ); +#else +double fabs(), gamma(), lgam(), exp(), log(), floor(); +#endif +extern double MAXLOG, MAXNUM; +extern int sgngam; + +double beta( a, b ) +double a, b; +{ +double y; +int sign; + +sign = 1; + +if( a <= 0.0 ) + { + if( a == floor(a) ) + goto over; + } +if( b <= 0.0 ) + { + if( b == floor(b) ) + goto over; + } + + +y = a + b; +if( fabs(y) > MAXGAM ) + { + y = lgam(y); + sign *= sgngam; /* keep track of the sign */ + y = lgam(b) - y; + sign *= sgngam; + y = lgam(a) + y; + sign *= sgngam; + if( y > MAXLOG ) + { +over: + mtherr( "beta", OVERFLOW ); + return( sign * MAXNUM ); + } + return( sign * exp(y) ); + } + +y = gamma(y); +if( y == 0.0 ) + goto over; + +if( a > b ) + { + y = gamma(a)/y; + y *= gamma(b); + } +else + { + y = gamma(b)/y; + y *= gamma(a); + } + +return(y); +} + + + +/* Natural log of |beta|. Return the sign of beta in sgngam. */ + +double lbeta( a, b ) +double a, b; +{ +double y; +int sign; + +sign = 1; + +if( a <= 0.0 ) + { + if( a == floor(a) ) + goto over; + } +if( b <= 0.0 ) + { + if( b == floor(b) ) + goto over; + } + + +y = a + b; +if( fabs(y) > MAXGAM ) + { + y = lgam(y); + sign *= sgngam; /* keep track of the sign */ + y = lgam(b) - y; + sign *= sgngam; + y = lgam(a) + y; + sign *= sgngam; + sgngam = sign; + return( y ); + } + +y = gamma(y); +if( y == 0.0 ) + { +over: + mtherr( "lbeta", OVERFLOW ); + return( sign * MAXNUM ); + } + +if( a > b ) + { + y = gamma(a)/y; + y *= gamma(b); + } +else + { + y = gamma(b)/y; + y *= gamma(a); + } + +if( y < 0 ) + { + sgngam = -1; + y = -y; + } +else + sgngam = 1; + +return( log(y) ); +} |