1 files changed, 122 insertions, 0 deletions
diff --git a/libm/float/betaf.c b/libm/float/betaf.c
new file mode 100644
index 000000000..7a1963191
--- /dev/null
+++ b/libm/float/betaf.c
@@ -0,0 +1,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);
+}