* mingwex/math/lgamma.c: New file.

	* mingwex/math/lgammaf.c: New file.
	* mingwex/math/lgammal.c: New file.
	* mingwex/math/tgamma.c: New file.
	* mingwex/math/tgammaf.c: New file.
	* mingwex/math/tgammal.c: New file.
	* mingwex/math/cephes_mconf (polevlf): Add float version.
	(p1evlf): Likewise.
	Define _CEPHES_USE_ERRNO.
	* mingwex/Makefile.in (MATH_DISTFILES): Add new files.
	(MATH_OBJS): Add new objects.
	* include/math.h (lgamma[fl]): Add prototypes.
	(tgamma[fl]): Add prototypes.

1 files changed, 253 insertions, 0 deletions

diff --git a/winsup/mingw/mingwex/math/lgammaf.c b/winsup/mingw/mingwex/math/lgammaf.c
new file mode 100644
index 000000000..20982f999
--- /dev/null
+++ b/winsup/mingw/mingwex/math/lgammaf.c
@@ -0,0 +1,253 @@
+/*							lgamf()
+ *
+ *	Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, __lgammaf_r();
+ * int* sgngamf;
+ * y = __lgammaf_r( x, sgngamf );
+ * 
+ * float x, y, lgammaf();
+ * y = lgammaf( x);
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the base e (2.718...) logarithm of the absolute
+ * value of the gamma function of the argument. In the reentrant
+ * version the sign (+1 or -1) of the gamma function is returned in
+ * variable referenced by sgngamf.
+ *
+ * For arguments greater than 6.5, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula.  Arguments between 0 and +6.5 are reduced by
+ * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational
+ * approximation.  The cosecant reflection formula is employed for
+ * arguments less than zero.
+ *
+ * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an
+ * error message.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *
+ * arithmetic      domain        # trials     peak         rms
+ *    IEEE        -100,+100       500,000    7.4e-7       6.8e-8
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ * The routine has low relative error for positive arguments.
+ *
+ * The following test used the relative error criterion.
+ *    IEEE    -2, +3              100000     4.0e-7      5.6e-8
+ *
+ */
+
+
+/*
+  Cephes Math Library Release 2.7:  July, 1998
+  Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier
+*/
+
+/*
+  26-11-2002 Modified for mingw.
+  Danny Smith <dannysmith@users.sourceforge.net>
+*/
+
+
+/* log gamma(x+2), -.5 < x < .5 */
+static const float B[] = {
+ 6.055172732649237E-004,
+-1.311620815545743E-003,
+ 2.863437556468661E-003,
+-7.366775108654962E-003,
+ 2.058355474821512E-002,
+-6.735323259371034E-002,
+ 3.224669577325661E-001,
+ 4.227843421859038E-001
+};
+
+/* log gamma(x+1), -.25 < x < .25 */
+static const float C[] = {
+ 1.369488127325832E-001,
+-1.590086327657347E-001,
+ 1.692415923504637E-001,
+-2.067882815621965E-001,
+ 2.705806208275915E-001,
+-4.006931650563372E-001,
+ 8.224670749082976E-001,
+-5.772156501719101E-001
+};
+
+/* log( sqrt( 2*pi ) ) */
+static const float LS2PI  =  0.91893853320467274178;
+#define MAXLGM 2.035093e36
+static const float PIINV =  0.318309886183790671538;
+
+#ifndef __MINGW32__
+#include "mconf.h"
+float floorf(float);
+float polevlf( float, float *, int );
+float p1evlf( float, float *, int );
+#else
+#include "cephes_mconf.h"
+#endif
+
+/* Reentrant version */ 
+/* Logarithm of gamma function */
+
+float __lgammaf_r( float x, int* sgngamf )
+{
+float p, q, w, z;
+float nx, tx;
+int i, direction;
+
+*sgngamf = 1;
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+#endif
+
+#ifdef INFINITIES
+if( !isfinite(x) )
+	return(x);
+#endif
+
+
+if( x < 0.0 )
+	{
+	q = -x;
+	w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */
+	p = floorf(q);
+	if( p == q )
+		{
+lgsing:
+		_SET_ERRNO(EDOM);
+		mtherr( "lgamf", SING );
+#ifdef INFINITIES
+		return (INFINITYF);
+#else
+	return( *sgngamf * MAXNUMF );
+#endif
+		}
+	i = p;
+	if( (i & 1) == 0 )
+		*sgngamf = -1;
+	else
+		*sgngamf = 1;
+	z = q - p;
+	if( z > 0.5 )
+		{
+		p += 1.0;
+		z = p - q;
+		}
+	z = q * sinf( PIF * z );
+	if( z == 0.0 )
+		goto lgsing;
+	z = -logf( PIINV*z ) - w;
+	return( z );
+	}
+
+if( x < 6.5 )
+	{
+	direction = 0;
+	z = 1.0;
+	tx = x;
+	nx = 0.0;
+	if( x >= 1.5 )
+		{
+		while( tx > 2.5 )
+			{
+			nx -= 1.0;
+			tx = x + nx;
+			z *=tx;
+			}
+		x += nx - 2.0;
+iv1r5:
+		p = x * polevlf( x, B, 7 );
+		goto cont;
+		}
+	if( x >= 1.25 )
+		{
+		z *= x;
+		x -= 1.0; /* x + 1 - 2 */
+		direction = 1;
+		goto iv1r5;
+		}
+	if( x >= 0.75 )
+		{
+		x -= 1.0;
+		p = x * polevlf( x, C, 7 );
+		q = 0.0;
+		goto contz;
+		}
+	while( tx < 1.5 )
+		{
+		if( tx == 0.0 )
+			goto lgsing;
+		z *=tx;
+		nx += 1.0;
+		tx = x + nx;
+		}
+	direction = 1;
+	x += nx - 2.0;
+	p = x * polevlf( x, B, 7 );
+
+cont:
+	if( z < 0.0 )
+		{
+		*sgngamf = -1;
+		z = -z;
+		}
+	else
+		{
+		*sgngamf = 1;
+		}
+	q = logf(z);
+	if( direction )
+		q = -q;
+contz:
+	return( p + q );
+	}
+
+if( x > MAXLGM )
+	{
+	_SET_ERRNO(ERANGE);
+	mtherr( "lgamf", OVERFLOW );
+#ifdef INFINITIES
+	return( *sgngamf * INFINITYF );
+#else
+	return( *sgngamf * MAXNUMF );
+#endif
+
+	}
+
+/* Note, though an asymptotic formula could be used for x >= 3,
+ * there is cancellation error in the following if x < 6.5.  */
+q = LS2PI - x;
+q += ( x - 0.5 ) * logf(x);
+
+if( x <= 1.0e4 )
+	{
+	z = 1.0/x;
+	p = z * z;
+	q += ((    6.789774945028216E-004 * p
+		 - 2.769887652139868E-003 ) * p
+		+  8.333316229807355E-002 ) * z;
+	}
+return( q );
+}
+
+/* This is the C99 version */
+
+float lgammaf(float x)
+{
+  int local_sgngamf=0;
+  return (__lgammaf_r(x, &local_sgngamf));
+}