uClibc now has a math library. muahahahaha!

 -Erik

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) );
+}