summaryrefslogtreecommitdiff
path: root/libm/double/beta.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/beta.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/beta.c')
-rw-r--r--libm/double/beta.c201
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) );
+}