summaryrefslogtreecommitdiff
path: root/libm/double/beta.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/beta.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/beta.c')
-rw-r--r--libm/double/beta.c201
1 files changed, 0 insertions, 201 deletions
diff --git a/libm/double/beta.c b/libm/double/beta.c
deleted file mode 100644
index 410760f32..000000000
--- a/libm/double/beta.c
+++ /dev/null
@@ -1,201 +0,0 @@
-/* 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) );
-}