* mingwex/math/powil.c: Rename powil to __powil.

	* mingwex/math/powl.c: Adjust declaration and call
	to __powil. Remove comment on powil.
	* mingwex/math/powi.c: New file.
	* mingwex/math/powif.c: New file.
	* mingwex/math/pow.c: New file.
	* mingwex/math/cephes_mconf.h. Add double and float
	versions of constants.
	(polevl): Add double precision function.
	(p1evl): Likewise.
	* mingwex/Makefile.in (MATH_DISTFILES): Add pow.c,
	powi.c, powif.c.
	(MATH_OBJS): Add pow.o,	powi.o, powif.o.

8 files changed, 1370 insertions, 67 deletions

diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog
index ad349dd72..bca53671c 100644
--- a/winsup/mingw/ChangeLog
+++ b/winsup/mingw/ChangeLog
@@ -1,3 +1,19 @@
+2002-10-07  Danny Smith  <dannysmith@users.sourceforge.net>
+
+	* mingwex/math/powil.c: Rename powil to __powil.
+	* mingwex/math/powl.c: Adjust declaration and call
+	to __powil. Remove comment on powil.
+	* mingwex/math/powi.c: New file.
+	* mingwex/math/powif.c: New file.
+	* mingwex/math/pow.c: New file.
+	* mingwex/math/cephes_mconf.h. Add double and float
+	versions of constants.
+	(polevl): Add double precision function.
+	(p1evl): Likewise.
+	* mingwex/Makefile.in (MATH_DISTFILES): Add pow.c,
+	powi.c, powif.c.
+	(MATH_OBJS): Add pow.o,	powi.o, powif.o.
+
 2002-10-03  Danny Smith  <dannysmith@users.sourceforge.net>
 
 	* include/cytpe.h (_imp____mbcur_max): Add missing ';'.
diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in
index 4c346d346..8ede9f9e9 100644
--- a/winsup/mingw/mingwex/Makefile.in
+++ b/winsup/mingw/mingwex/Makefile.in
@@ -50,8 +50,9 @@ MATH_DISTFILES = \
 	log10f.S log10l.S log1p.S log1pf.S log1pl.S log2.S log2f.S \
 	log2l.S logb.c logbf.c logbl.c logf.S logl.S lrint.c lrintf.c \
 	lrintl.c lround.c lroundf.c lroundl.c modff.c modfl.c \
-	nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c powf.c powil.c \
-	powl.c remainder.S remainderf.S remainderl.S remquo.S \
+	nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c \
+	pow.c powf.c powi.c powif.c powil.c powl.c \
+	remainder.S remainderf.S remainderl.S remquo.S \
 	remquof.S remquol.S rint.c rintf.c rintl.c round.c roundf.c \
 	roundl.c scalbn.S scalbnf.S scalbnl.S signbit.c signbitf.c \
 	signbitl.c sinf.S sinhf.c sinhl.c sinl.S  sqrtf.c sqrtl.c \
@@ -115,8 +116,9 @@ MATH_OBJS = \
 	log10f.o log10l.o log1p.o log1pf.o log1pl.o log2.o log2f.o \
 	log2l.o logb.o logbf.o logbl.o logf.o logl.o lrint.o lrintf.o \
 	lrintl.o lround.o lroundf.o lroundl.o modff.o modfl.o \
-	nearbyint.o nearbyintf.o nearbyintl.o nextafterf.o powf.o powil.o \
-	powl.o remainder.o remainderf.o remainderl.o remquo.o \
+	nearbyint.o nearbyintf.o nearbyintl.o nextafterf.o \
+	pow.o powf.o powi.o powif.o powil.o powl.o \
+	remainder.o remainderf.o remainderl.o remquo.o \
 	remquof.o remquol.o rint.o rintf.o rintl.o round.o roundf.o \
 	roundl.o scalbn.o scalbnf.o scalbnl.o signbit.o signbitf.o \
 	signbitl.o sinf.o sinhf.o sinhl.o sinl.o sqrtf.o sqrtl.o \
diff --git a/winsup/mingw/mingwex/math/cephes_mconf.h b/winsup/mingw/mingwex/math/cephes_mconf.h
index ba8400f34..1dda63d53 100644
--- a/winsup/mingw/mingwex/math/cephes_mconf.h
+++ b/winsup/mingw/mingwex/math/cephes_mconf.h
@@ -1,8 +1,45 @@
 #include <math.h>
 #include <errno.h>
 
+
+#define IBMPC 1
+#define ANSIPROT 1
+#define MINUSZERO 1
+#define INFINITIES 1
+#define NANS 1
+#define DENORMAL 1
+#define VOLATILE
+#define mtherr(fname, code) 
+#define XPD 0,
+
+#ifdef _CEPHES_USE_ERRNO
+#define _SET_ERRNO(x) errno = (x)
+#else
+#define _SET_ERRNO(x)
+#endif
+
 /* constants used by cephes functions */
 
+/* double */
+#define MAXNUM	1.7976931348623158E308
+#define MAXLOG	7.09782712893383996843E2
+#define MINLOG	-7.08396418532264106224E2
+#define LOGE2	6.93147180559945309417E-1
+#define LOG2E	1.44269504088896340736
+#define PI	3.14159265358979323846
+#define PIO2	1.57079632679489661923
+#define PIO4	7.85398163397448309616E-1
+
+#define NEGZERO (-0.0)
+extern double __INF;
+#undef INFINITY
+#define INFINITY (__INF)
+extern double __QNAN;
+#undef NAN
+#define NAN (__QNAN)
+
+
+/*long double*/
 #define MAXNUML 1.189731495357231765021263853E4932L
 #define MAXLOGL	1.1356523406294143949492E4L
 #define MINLOGL	-1.13994985314888605586758E4L
@@ -17,27 +54,133 @@
 #define isnanl isnan
 #define signbitl signbit
 
-
-#define IBMPC 1
-#define ANSIPROT 1
-#define MINUSZERO 1
-#define INFINITIES 1
-#define NANS 1
-#define DENORMAL 1
 #define NEGZEROL (-0.0L)
 extern long double __INFL;
 #define INFINITYL (__INFL)
 extern long double __QNANL;
 #define NANL (__QNANL)
-#define VOLATILE
-#define mtherr(fname, code) 
-#define XPD 0,
 
-#ifdef _CEPHES_USE_ERRNO
-#define _SET_ERRNO(x) errno = (x)
-#else
-#define _SET_ERRNO(x)
-#endif
+/* float */
+
+#define MAXNUMF	3.4028234663852885981170418348451692544e38F
+#define MAXLOGF	88.72283905206835F
+#define MINLOGF	-103.278929903431851103F /* log(2^-149) */
+#define LOG2EF	1.44269504088896341F
+#define LOGE2F	0.693147180559945309F
+#define PIF	3.141592653589793238F
+#define PIO2F	1.5707963267948966192F
+#define PIO4F	0.7853981633974483096F
+
+#define isfinitef isfinite
+#define isinff isinf
+#define isnanf isnan
+#define signbitf signbit
+
+#define NEGZEROF (-0.0F)
+extern float __INFF;
+#define INFINITYF (__INFF)
+extern float __QNANF;
+#define NANF (__QNANF)
+
+
+/* double */
+
+/*
+Cephes Math Library Release 2.2:  July, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+/*							polevl.c
+ *							p1evl.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * double x, y, coef[N+1], polevl[];
+ *
+ * y = polevl( x, coef, N );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Evaluates polynomial of degree N:
+ *
+ *                     2          N
+ * y  =  C  + C x + C x  +...+ C x
+ *        0    1     2          N
+ *
+ * Coefficients are stored in reverse order:
+ *
+ * coef[0] = C  , ..., coef[N] = C  .
+ *            N                   0
+ *
+ *  The function p1evl() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevl().
+ *
+ *
+ * SPEED:
+ *
+ * In the interest of speed, there are no checks for out
+ * of bounds arithmetic.  This routine is used by most of
+ * the functions in the library.  Depending on available
+ * equipment features, the user may wish to rewrite the
+ * program in microcode or assembly language.
+ *
+ */
+
+/* Polynomial evaluator:
+ *  P[0] x^n  +  P[1] x^(n-1)  +  ...  +  P[n]
+ */
+static  __inline__ double polevl( x, p, n )
+double x;
+const void *p;
+int n;
+{
+register double y;
+register double *P = (double *)p;
+
+y = *P++;
+do
+	{
+	y = y * x + *P++;
+	}
+while( --n );
+return(y);
+}
+
+
+
+/* Polynomial evaluator:
+ *  x^n  +  P[0] x^(n-1)  +  P[1] x^(n-2)  +  ...  +  P[n]
+ */
+static __inline__  double p1evl( x, p, n )
+double x;
+const void *p;
+int n;
+{
+register double y;
+register double *P = (double *)p;
+
+n -= 1;
+y = x + *P++;
+do
+	{
+	y = y * x + *P++;
+	}
+while( --n );
+return( y );
+}
+
+
+/* long double */
 /*
 Cephes Math Library Release 2.2:  July, 1992
 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
diff --git a/winsup/mingw/mingwex/math/pow.c b/winsup/mingw/mingwex/math/pow.c
new file mode 100644
index 000000000..1fa548e5e
--- /dev/null
+++ b/winsup/mingw/mingwex/math/pow.c
@@ -0,0 +1,781 @@
+/*							pow.c
+ *
+ *	Power function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, z, pow();
+ *
+ * z = pow( x, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes x raised to the yth power.  Analytically,
+ *
+ *      x**y  =  exp( y log(x) ).
+ *
+ * Following Cody and Waite, this program uses a lookup table
+ * of 2**-i/16 and pseudo extended precision arithmetic to
+ * obtain an extra three bits of accuracy in both the logarithm
+ * and the exponential.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -26,26       30000      4.2e-16      7.7e-17
+ *    DEC      -26,26       60000      4.8e-17      9.1e-18
+ * 1/26 < x < 26, with log(x) uniformly distributed.
+ * -26 < y < 26, y uniformly distributed.
+ *    IEEE     0,8700       30000      1.5e-14      2.1e-15
+ * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * pow overflow     x**y > MAXNUM      INFINITY
+ * pow underflow   x**y < 1/MAXNUM       0.0
+ * pow domain      x<0 and y noninteger  0.0
+ *
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1995, 2000 by Stephen L. Moshier
+*/
+
+/*
+Modified for mingw
+2002-09-27 Danny Smith <dannysmith@users.sourceforge.net>
+*/
+
+#ifdef __MINGW32__
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+static char fname[] = {"pow"};
+#endif
+
+#ifndef _SET_ERRNO
+#define _SET_ERRNO(x)
+#endif
+
+#define SQRTH 0.70710678118654752440
+
+#ifdef UNK
+static double P[] = {
+  4.97778295871696322025E-1,
+  3.73336776063286838734E0,
+  7.69994162726912503298E0,
+  4.66651806774358464979E0
+};
+static double Q[] = {
+/* 1.00000000000000000000E0, */
+  9.33340916416696166113E0,
+  2.79999886606328401649E1,
+  3.35994905342304405431E1,
+  1.39995542032307539578E1
+};
+/* 2^(-i/16), IEEE precision */
+static double A[] = {
+  1.00000000000000000000E0,
+  9.57603280698573700036E-1,
+  9.17004043204671215328E-1,
+  8.78126080186649726755E-1,
+  8.40896415253714502036E-1,
+  8.05245165974627141736E-1,
+  7.71105412703970372057E-1,
+  7.38413072969749673113E-1,
+  7.07106781186547572737E-1,
+  6.77127773468446325644E-1,
+  6.48419777325504820276E-1,
+  6.20928906036742001007E-1,
+  5.94603557501360513449E-1,
+  5.69394317378345782288E-1,
+  5.45253866332628844837E-1,
+  5.22136891213706877402E-1,
+  5.00000000000000000000E-1
+};
+static double B[] = {
+ 0.00000000000000000000E0,
+ 1.64155361212281360176E-17,
+ 4.09950501029074826006E-17,
+ 3.97491740484881042808E-17,
+-4.83364665672645672553E-17,
+ 1.26912513974441574796E-17,
+ 1.99100761573282305549E-17,
+-1.52339103990623557348E-17,
+ 0.00000000000000000000E0
+};
+static double R[] = {
+ 1.49664108433729301083E-5,
+ 1.54010762792771901396E-4,
+ 1.33335476964097721140E-3,
+ 9.61812908476554225149E-3,
+ 5.55041086645832347466E-2,
+ 2.40226506959099779976E-1,
+ 6.93147180559945308821E-1
+};
+
+#define douba(k) A[k]
+#define doubb(k) B[k]
+#define MEXP 16383.0
+#ifdef DENORMAL
+#define MNEXP -17183.0
+#else
+#define MNEXP -16383.0
+#endif
+#endif
+
+#ifdef DEC
+static unsigned short P[] = {
+0037776,0156313,0175332,0163602,
+0040556,0167577,0052366,0174245,
+0040766,0062753,0175707,0055564,
+0040625,0052035,0131344,0155636,
+};
+static unsigned short Q[] = {
+/*0040200,0000000,0000000,0000000,*/
+0041025,0052644,0154404,0105155,
+0041337,0177772,0007016,0047646,
+0041406,0062740,0154273,0020020,
+0041137,0177054,0106127,0044555,
+};
+static unsigned short A[] = {
+0040200,0000000,0000000,0000000,
+0040165,0022575,0012444,0103314,
+0040152,0140306,0163735,0022071,
+0040140,0146336,0166052,0112341,
+0040127,0042374,0145326,0116553,
+0040116,0022214,0012437,0102201,
+0040105,0063452,0010525,0003333,
+0040075,0004243,0117530,0006067,
+0040065,0002363,0031771,0157145,
+0040055,0054076,0165102,0120513,
+0040045,0177326,0124661,0050471,
+0040036,0172462,0060221,0120422,
+0040030,0033760,0050615,0134251,
+0040021,0141723,0071653,0010703,
+0040013,0112701,0161752,0105727,
+0040005,0125303,0063714,0044173,
+0040000,0000000,0000000,0000000
+};
+static unsigned short B[] = {
+0000000,0000000,0000000,0000000,
+0021473,0040265,0153315,0140671,
+0121074,0062627,0042146,0176454,
+0121413,0003524,0136332,0066212,
+0121767,0046404,0166231,0012553,
+0121257,0015024,0002357,0043574,
+0021736,0106532,0043060,0056206,
+0121310,0020334,0165705,0035326,
+0000000,0000000,0000000,0000000
+};
+
+static unsigned short R[] = {
+0034173,0014076,0137624,0115771,
+0035041,0076763,0003744,0111311,
+0035656,0141766,0041127,0074351,
+0036435,0112533,0073611,0116664,
+0037143,0054106,0134040,0152223,
+0037565,0176757,0176026,0025551,
+0040061,0071027,0173721,0147572
+};
+
+/*
+static double R[] = {
+0.14928852680595608186e-4,
+0.15400290440989764601e-3,
+0.13333541313585784703e-2,
+0.96181290595172416964e-2,
+0.55504108664085595326e-1,
+0.24022650695909537056e0,
+0.69314718055994529629e0
+};
+*/
+#define douba(k) (*(double *)&A[(k)<<2])
+#define doubb(k) (*(double *)&B[(k)<<2])
+#define MEXP 2031.0
+#define MNEXP -2031.0
+#endif
+
+#ifdef IBMPC
+static const unsigned short P[] = {
+0x5cf0,0x7f5b,0xdb99,0x3fdf,
+0xdf15,0xea9e,0xddef,0x400d,
+0xeb6f,0x7f78,0xccbd,0x401e,
+0x9b74,0xb65c,0xaa83,0x4012,
+};
+static const unsigned short Q[] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0x914e,0x9b20,0xaab4,0x4022,
+0xc9f5,0x41c1,0xffff,0x403b,
+0x6402,0x1b17,0xccbc,0x4040,
+0xe92e,0x918a,0xffc5,0x402b,
+};
+static const unsigned short A[] = {
+0x0000,0x0000,0x0000,0x3ff0,
+0x90da,0xa2a4,0xa4af,0x3fee,
+0xa487,0xdcfb,0x5818,0x3fed,
+0x529c,0xdd85,0x199b,0x3fec,
+0xd3ad,0x995a,0xe89f,0x3fea,
+0xf090,0x82a3,0xc491,0x3fe9,
+0xa0db,0x422a,0xace5,0x3fe8,
+0x0187,0x73eb,0xa114,0x3fe7,
+0x3bcd,0x667f,0xa09e,0x3fe6,
+0x5429,0xdd48,0xab07,0x3fe5,
+0x2a27,0xd536,0xbfda,0x3fe4,
+0x3422,0x4c12,0xdea6,0x3fe3,
+0xb715,0x0a31,0x06fe,0x3fe3,
+0x6238,0x6e75,0x387a,0x3fe2,
+0x517b,0x3c7d,0x72b8,0x3fe1,
+0x890f,0x6cf9,0xb558,0x3fe0,
+0x0000,0x0000,0x0000,0x3fe0
+};
+static const unsigned short B[] = {
+0x0000,0x0000,0x0000,0x0000,
+0x3707,0xd75b,0xed02,0x3c72,
+0xcc81,0x345d,0xa1cd,0x3c87,
+0x4b27,0x5686,0xe9f1,0x3c86,
+0x6456,0x13b2,0xdd34,0xbc8b,
+0x42e2,0xafec,0x4397,0x3c6d,
+0x82e4,0xd231,0xf46a,0x3c76,
+0x8a76,0xb9d7,0x9041,0xbc71,
+0x0000,0x0000,0x0000,0x0000
+};
+static const unsigned short R[] = {
+0x937f,0xd7f2,0x6307,0x3eef,
+0x9259,0x60fc,0x2fbe,0x3f24,
+0xef1d,0xc84a,0xd87e,0x3f55,
+0x33b7,0x6ef1,0xb2ab,0x3f83,
+0x1a92,0xd704,0x6b08,0x3fac,
+0xc56d,0xff82,0xbfbd,0x3fce,
+0x39ef,0xfefa,0x2e42,0x3fe6
+};
+
+#define douba(k) (*(double *)&A[(k)<<2])
+#define doubb(k) (*(double *)&B[(k)<<2])
+#define MEXP 16383.0
+#ifdef DENORMAL
+#define MNEXP -17183.0
+#else
+#define MNEXP -16383.0
+#endif
+#endif
+
+#ifdef MIEEE
+static unsigned short P[] = {
+0x3fdf,0xdb99,0x7f5b,0x5cf0,
+0x400d,0xddef,0xea9e,0xdf15,
+0x401e,0xccbd,0x7f78,0xeb6f,
+0x4012,0xaa83,0xb65c,0x9b74
+};
+static unsigned short Q[] = {
+0x4022,0xaab4,0x9b20,0x914e,
+0x403b,0xffff,0x41c1,0xc9f5,
+0x4040,0xccbc,0x1b17,0x6402,
+0x402b,0xffc5,0x918a,0xe92e
+};
+static unsigned short A[] = {
+0x3ff0,0x0000,0x0000,0x0000,
+0x3fee,0xa4af,0xa2a4,0x90da,
+0x3fed,0x5818,0xdcfb,0xa487,
+0x3fec,0x199b,0xdd85,0x529c,
+0x3fea,0xe89f,0x995a,0xd3ad,
+0x3fe9,0xc491,0x82a3,0xf090,
+0x3fe8,0xace5,0x422a,0xa0db,
+0x3fe7,0xa114,0x73eb,0x0187,
+0x3fe6,0xa09e,0x667f,0x3bcd,
+0x3fe5,0xab07,0xdd48,0x5429,
+0x3fe4,0xbfda,0xd536,0x2a27,
+0x3fe3,0xdea6,0x4c12,0x3422,
+0x3fe3,0x06fe,0x0a31,0xb715,
+0x3fe2,0x387a,0x6e75,0x6238,
+0x3fe1,0x72b8,0x3c7d,0x517b,
+0x3fe0,0xb558,0x6cf9,0x890f,
+0x3fe0,0x0000,0x0000,0x0000
+};
+static unsigned short B[] = {
+0x0000,0x0000,0x0000,0x0000,
+0x3c72,0xed02,0xd75b,0x3707,
+0x3c87,0xa1cd,0x345d,0xcc81,
+0x3c86,0xe9f1,0x5686,0x4b27,
+0xbc8b,0xdd34,0x13b2,0x6456,
+0x3c6d,0x4397,0xafec,0x42e2,
+0x3c76,0xf46a,0xd231,0x82e4,
+0xbc71,0x9041,0xb9d7,0x8a76,
+0x0000,0x0000,0x0000,0x0000
+};
+static unsigned short R[] = {
+0x3eef,0x6307,0xd7f2,0x937f,
+0x3f24,0x2fbe,0x60fc,0x9259,
+0x3f55,0xd87e,0xc84a,0xef1d,
+0x3f83,0xb2ab,0x6ef1,0x33b7,
+0x3fac,0x6b08,0xd704,0x1a92,
+0x3fce,0xbfbd,0xff82,0xc56d,
+0x3fe6,0x2e42,0xfefa,0x39ef
+};
+
+#define douba(k) (*(double *)&A[(k)<<2])
+#define doubb(k) (*(double *)&B[(k)<<2])
+#define MEXP 16383.0
+#ifdef DENORMAL
+#define MNEXP -17183.0
+#else
+#define MNEXP -16383.0
+#endif
+#endif
+
+/* log2(e) - 1 */
+#define LOG2EA 0.44269504088896340736
+
+#define F W
+#define Fa Wa
+#define Fb Wb
+#define G W
+#define Ga Wa
+#define Gb u
+#define H W
+#define Ha Wb
+#define Hb Wb
+
+#ifdef __MINGW32__
+static __inline__ double reduc( double );
+extern double __powi ( double, int );
+extern double pow ( double x,  double y);
+
+#else /* __MINGW32__ */
+
+#ifdef ANSIPROT
+extern double floor ( double );
+extern double fabs ( double );
+extern double frexp ( double, int * );
+extern double ldexp ( double, int );
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double __powi ( double, int );
+extern int signbit ( double );
+extern int isnan ( double );
+extern int isfinite ( double );
+static double reduc ( double );
+#else
+double floor(), fabs(), frexp(), ldexp();
+double polevl(), p1evl(), __powi();
+int signbit(), isnan(), isfinite();
+static double reduc();
+#endif
+extern double MAXNUM;
+#ifdef INFINITIES
+extern double INFINITY;
+#endif
+#ifdef NANS
+extern double NAN;
+#endif
+#ifdef MINUSZERO
+extern double NEGZERO;
+#endif
+
+#endif /* __MINGW32__ */
+
+double pow( x, y )
+double x, y;
+{
+double w, z, W, Wa, Wb, ya, yb, u;
+/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
+double aw, ay, wy;
+int e, i, nflg, iyflg, yoddint;
+
+if( y == 0.0 )
+	return( 1.0 );
+#ifdef NANS
+if( isnan(x) || isnan(y) )
+	{
+	_SET_ERRNO (EDOM);
+	return( x + y );
+	}
+#endif
+if( y == 1.0 )
+	return( x );
+
+
+#ifdef INFINITIES
+if( !isfinite(y) && (x == 1.0 || x == -1.0) )
+	{
+	mtherr( "pow", DOMAIN );
+#ifdef NANS
+	return( NAN );
+#else
+	return( INFINITY );
+#endif
+	}
+#endif
+
+if( x == 1.0 )
+	return( 1.0 );
+
+if( y >= MAXNUM )
+	{
+	_SET_ERRNO (ERANGE);
+#ifdef INFINITIES
+	if( x > 1.0 )
+		return( INFINITY );
+#else
+	if( x > 1.0 )
+		return( MAXNUM );
+#endif
+	if( x > 0.0 && x < 1.0 )
+		return( 0.0);
+	if( x < -1.0 )
+		{
+#ifdef INFINITIES
+		return( INFINITY );
+#else
+		return( MAXNUM );
+#endif
+		}
+	if( x > -1.0 && x < 0.0 )
+		return( 0.0 );
+	}
+if( y <= -MAXNUM )
+	{
+	_SET_ERRNO (ERANGE);
+	if( x > 1.0 )
+		return( 0.0 );
+#ifdef INFINITIES
+	if( x > 0.0 && x < 1.0 )
+		return( INFINITY );
+#else
+	if( x > 0.0 && x < 1.0 )
+		return( MAXNUM );
+#endif
+	if( x < -1.0 )
+		return( 0.0 );
+#ifdef INFINITIES
+	if( x > -1.0 && x < 0.0 )
+		return( INFINITY );
+#else
+	if( x > -1.0 && x < 0.0 )
+		return( MAXNUM );
+#endif
+	}
+if( x >= MAXNUM )
+	{
+#if INFINITIES
+	if( y > 0.0 )
+		return( INFINITY );
+#else
+	if( y > 0.0 )
+		return( MAXNUM );
+#endif
+	return(0.0);
+	}
+/* Set iyflg to 1 if y is an integer.  */
+iyflg = 0;
+w = floor(y);
+if( w == y )
+	iyflg = 1;
+
+/* Test for odd integer y.  */
+yoddint = 0;
+if( iyflg )
+	{
+	ya = fabs(y);
+	ya = floor(0.5 * ya);
+	yb = 0.5 * fabs(w);
+	if( ya != yb )
+		yoddint = 1;
+	}
+
+if( x <= -MAXNUM )
+	{
+	if( y > 0.0 )
+		{
+#ifdef INFINITIES
+		if( yoddint )
+			return( -INFINITY );
+		return( INFINITY );
+#else
+		if( yoddint )
+			return( -MAXNUM );
+		return( MAXNUM );
+#endif
+		}
+	if( y < 0.0 )
+		{
+#ifdef MINUSZERO
+		if( yoddint )
+			return( NEGZERO );
+#endif
+		return( 0.0 );
+		}
+ 	}
+
+nflg = 0;	/* flag = 1 if x<0 raised to integer power */
+if( x <= 0.0 )
+	{
+	if( x == 0.0 )
+		{
+		if( y < 0.0 )
+			{
+#ifdef MINUSZERO
+			if( signbit(x) && yoddint )
+				return( -INFINITY );
+#endif
+#ifdef INFINITIES
+			return( INFINITY );
+#else
+			return( MAXNUM );
+#endif
+			}
+		if( y > 0.0 )
+			{
+#ifdef MINUSZERO
+			if( signbit(x) && yoddint )
+				return( NEGZERO );
+#endif
+			return( 0.0 );
+			}
+		return( 1.0 );
+		}
+	else
+		{
+		if( iyflg == 0 )
+			{ /* noninteger power of negative number */
+			mtherr( fname, DOMAIN );
+			_SET_ERRNO (EDOM);
+#ifdef NANS
+			return(NAN);
+#else
+			return(0.0L);
+#endif
+			}
+		nflg = 1;
+		}
+	}
+
+/* Integer power of an integer.  */
+
+if( iyflg )
+	{
+	i = w;
+	w = floor(x);
+	if( (w == x) && (fabs(y) < 32768.0) )
+		{
+		w = __powi( x, (int) y );
+		return( w );
+		}
+	}
+
+if( nflg )
+	x = fabs(x);
+
+/* For results close to 1, use a series expansion.  */
+w = x - 1.0;
+aw = fabs(w);
+ay = fabs(y);
+wy = w * y;
+ya = fabs(wy);
+if((aw <= 1.0e-3 && ay <= 1.0)
+   || (ya <= 1.0e-3 && ay >= 1.0))
+	{
+	z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
+		+ 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
+	goto done;
+	}
+/* These are probably too much trouble.  */
+#if 0
+w = y * log(x);
+if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
+  {
+    z = ((((((
+    w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
+    goto done;
+  }
+
+if(ya <= 1.0e-3 && aw <= 1.0e-4)
+  {
+    z = (((((
+	     wy*1./720.
+	     + (-w*1./48. + 1./120.) )*wy
+	    + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
+	   + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
+	  + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
+	 + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
+	   + wy + 1.0;
+    goto done;
+  }
+#endif
+
+/* separate significand from exponent */
+x = frexp( x, &e );
+
+#if 0
+/* For debugging, check for gross overflow. */
+if( (e * y)  > (MEXP + 1024) )
+	goto overflow;
+#endif
+
+/* Find significand of x in antilog table A[]. */
+i = 1;
+if( x <= douba(9) )
+	i = 9;
+if( x <= douba(i+4) )
+	i += 4;
+if( x <= douba(i+2) )
+	i += 2;
+if( x >= douba(1) )
+	i = -1;
+i += 1;
+
+
+/* Find (x - A[i])/A[i]
+ * in order to compute log(x/A[i]):
+ *
+ * log(x) = log( a x/a ) = log(a) + log(x/a)
+ *
+ * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
+ */
+x -= douba(i);
+x -= doubb(i/2);
+x /= douba(i);
+
+
+/* rational approximation for log(1+v):
+ *
+ * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
+ */
+z = x*x;
+w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
+w = w - ldexp( z, -1 );   /*  w - 0.5 * z  */
+
+/* Convert to base 2 logarithm:
+ * multiply by log2(e)
+ */
+w = w + LOG2EA * w;
+/* Note x was not yet added in
+ * to above rational approximation,
+ * so do it now, while multiplying
+ * by log2(e).
+ */
+z = w + LOG2EA * x;
+z = z + x;
+
+/* Compute exponent term of the base 2 logarithm. */
+w = -i;
+w = ldexp( w, -4 );	/* divide by 16 */
+w += e;
+/* Now base 2 log of x is w + z. */
+
+/* Multiply base 2 log by y, in extended precision. */
+
+/* separate y into large part ya
+ * and small part yb less than 1/16
+ */
+ya = reduc(y);
+yb = y - ya;
+
+
+F = z * y  +  w * yb;
+Fa = reduc(F);
+Fb = F - Fa;
+
+G = Fa + w * ya;
+Ga = reduc(G);
+Gb = G - Ga;
+
+H = Fb + Gb;
+Ha = reduc(H);
+w = ldexp( Ga+Ha, 4 );
+
+/* Test the power of 2 for overflow */
+if( w > MEXP )
+	{
+#ifndef INFINITIES
+	mtherr( fname, OVERFLOW );
+#endif
+#ifdef INFINITIES
+	if( nflg && yoddint )
+	  return( -INFINITY );
+	return( INFINITY );
+#else
+	if( nflg && yoddint )
+	  return( -MAXNUM );
+	return( MAXNUM );
+#endif
+	}
+
+if( w < (MNEXP - 1) )
+	{
+#ifndef DENORMAL
+	mtherr( fname, UNDERFLOW );
+#endif
+#ifdef MINUSZERO
+	if( nflg && yoddint )
+	  return( NEGZERO );
+#endif
+	return( 0.0 );
+	}
+
+e = w;
+Hb = H - Ha;
+
+if( Hb > 0.0 )
+	{
+	e += 1;
+	Hb -= 0.0625;
+	}
+
+/* Now the product y * log2(x)  =  Hb + e/16.0.
+ *
+ * Compute base 2 exponential of Hb,
+ * where -0.0625 <= Hb <= 0.
+ */
+z = Hb * polevl( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
+
+/* Express e/16 as an integer plus a negative number of 16ths.
+ * Find lookup table entry for the fractional power of 2.
+ */
+if( e < 0 )
+	i = 0;
+else
+	i = 1;
+i = e/16 + i;
+e = 16*i - e;
+w = douba( e );
+z = w + w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
+z = ldexp( z, i );  /* multiply by integer power of 2 */
+
+done:
+
+/* Negate if odd integer power of negative number */
+if( nflg && yoddint )
+	{
+#ifdef MINUSZERO
+	if( z == 0.0 )
+		z = NEGZERO;
+	else
+#endif
+		z = -z;
+	}
+return( z );
+}
+
+
+/* Find a multiple of 1/16 that is within 1/16 of x. */
+static __inline__ double reduc(x)
+double x;
+{
+double t;
+
+t = ldexp( x, 4 );
+t = floor( t );
+t = ldexp( t, -4 );
+return(t);
+}
diff --git a/winsup/mingw/mingwex/math/powi.c b/winsup/mingw/mingwex/math/powi.c
new file mode 100644
index 000000000..9dd0c3d82
--- /dev/null
+++ b/winsup/mingw/mingwex/math/powi.c
@@ -0,0 +1,200 @@
+/*							powi.c
+ *
+ *	Real raised to integer power
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, __powi();
+ * int n;
+ *
+ * y = __powi( x, n );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns argument x raised to the nth power.
+ * The routine efficiently decomposes n as a sum of powers of
+ * two. The desired power is a product of two-to-the-kth
+ * powers of x.  Thus to compute the 32767 power of x requires
+ * 28 multiplications instead of 32767 multiplications.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *                      Relative error:
+ * arithmetic   x domain   n domain  # trials      peak         rms
+ *    DEC       .04,26     -26,26    100000       2.7e-16     4.3e-17
+ *    IEEE      .04,26     -26,26     50000       2.0e-15     3.8e-16
+ *    IEEE        1,2    -1022,1023   50000       8.6e-14     1.6e-14
+ *
+ * Returns MAXNUM on overflow, zero on underflow.
+ *
+ */
+
+/*							powi.c	*/
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1995, 2000 by Stephen L. Moshier
+*/
+
+/*
+Modified for mingw
+2002-07-22 Danny Smith <dannysmith@users.sourceforge.net>
+*/
+
+#ifdef __MINGW32__
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+#ifdef ANSIPROT
+extern double log ( double );
+extern double frexp ( double, int * );
+extern int signbit ( double );
+#else
+double log(), frexp();
+int signbit();
+#endif
+extern double NEGZERO, INFINITY, MAXNUM, MAXLOG, MINLOG, LOGE2;
+#endif /* __MINGW32__ */
+
+#ifndef _SET_ERRNO
+#define _SET_ERRNO(x)
+#endif
+
+double __powi( x, nn )
+double x;
+int nn;
+{
+int n, e, sign, asign, lx;
+double w, y, s;
+
+/* See pow.c for these tests.  */
+if( x == 0.0 )
+	{
+	if( nn == 0 )
+		return( 1.0 );
+	else if( nn < 0 )
+	    return( INFINITY );
+	else
+	  {
+	    if( nn & 1 )
+	      return( x );
+	    else
+	      return( 0.0 );
+	  }
+	}
+
+if( nn == 0 )
+	return( 1.0 );
+
+if( nn == -1 )
+	return( 1.0/x );
+
+if( x < 0.0 )
+	{
+	asign = -1;
+	x = -x;
+	}
+else
+	asign = 0;
+
+
+if( nn < 0 )
+	{
+	sign = -1;
+	n = -nn;
+	}
+else
+	{
+	sign = 1;
+	n = nn;
+	}
+
+/* Even power will be positive. */
+if( (n & 1) == 0 )
+	asign = 0;
+
+/* Overflow detection */
+
+/* Calculate approximate logarithm of answer */
+s = frexp( x, &lx );
+e = (lx - 1)*n;
+if( (e == 0) || (e > 64) || (e < -64) )
+	{
+	s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
+	s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
+	}
+else
+	{
+	s = LOGE2 * e;
+	}
+
+if( s > MAXLOG )
+	{
+	mtherr( "powi", OVERFLOW );
+	_SET_ERRNO(ERANGE);
+	y = INFINITY;
+	goto done;
+	}
+
+#if DENORMAL
+if( s < MINLOG )
+	{
+	y = 0.0;
+	goto done;
+	}
+
+/* Handle tiny denormal answer, but with less accuracy
+ * since roundoff error in 1.0/x will be amplified.
+ * The precise demarcation should be the gradual underflow threshold.
+ */
+if( (s < (-MAXLOG+2.0)) && (sign < 0) )
+	{
+	x = 1.0/x;
+	sign = -sign;
+	}
+#else
+/* do not produce denormal answer */
+if( s < -MAXLOG )
+	return(0.0);
+#endif
+
+
+/* First bit of the power */
+if( n & 1 )
+	y = x;
+		
+else
+	y = 1.0;
+
+w = x;
+n >>= 1;
+while( n )
+	{
+	w = w * w;	/* arg to the 2-to-the-kth power */
+	if( n & 1 )	/* if that bit is set, then include in product */
+		y *= w;
+	n >>= 1;
+	}
+
+if( sign < 0 )
+	y = 1.0/y;
+
+done:
+
+if( asign )
+	{
+	/* odd power of negative number */
+	if( y == 0.0 )
+		y = NEGZERO;
+	else
+		y = -y;
+	}
+return(y);
+}
diff --git a/winsup/mingw/mingwex/math/powif.c b/winsup/mingw/mingwex/math/powif.c
new file mode 100644
index 000000000..160fb5476
--- /dev/null
+++ b/winsup/mingw/mingwex/math/powif.c
@@ -0,0 +1,198 @@
+/*							powi.c
+ *
+ *	Real raised to integer power
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, __powif();
+ * int n;
+ *
+ * y = powi( x, n );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns argument x raised to the nth power.
+ * The routine efficiently decomposes n as a sum of powers of
+ * two. The desired power is a product of two-to-the-kth
+ * powers of x.  Thus to compute the 32767 power of x requires
+ * 28 multiplications instead of 32767 multiplications.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *                      Relative error:
+ * arithmetic   x domain   n domain  # trials      peak         rms
+ *    DEC       .04,26     -26,26    100000       2.7e-16     4.3e-17
+ *    IEEE      .04,26     -26,26     50000       2.0e-15     3.8e-16
+ *    IEEE        1,2    -1022,1023   50000       8.6e-14     1.6e-14
+ *
+ * Returns MAXNUM on overflow, zero on underflow.
+ *
+ */
+
+/*							powi.c	*/
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1995, 2000 by Stephen L. Moshier
+*/
+
+/*
+Modified for float from powi.c and adapted to mingw
+2002-10-01 Danny Smith <dannysmith@users.sourceforge.net>
+*/
+
+#ifdef __MINGW32__
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+#ifdef ANSIPROT
+extern float logf ( float );
+extern float frexpf ( float, int * );
+extern int signbitf ( float );
+#else
+float logf(), frexpf();
+int signbitf();
+#endif
+extern float NEGZEROF, INFINITYF, MAXNUMF, MAXLOGF, MINLOGF, LOGE2F;
+#endif /* __MINGW32__ */
+
+#ifndef _SET_ERRNO
+#define _SET_ERRNO(x)
+#endif
+
+float __powif( float x, int nn )
+{
+int n, e, sign, asign, lx;
+float w, y, s;
+
+/* See pow.c for these tests.  */
+if( x == 0.0F )
+	{
+	if( nn == 0 )
+		return( 1.0F );
+	else if( nn < 0 )
+	    return( INFINITYF );
+	else
+	  {
+	    if( nn & 1 )
+	      return( x );
+	    else
+	      return( 0.0 );
+	  }
+	}
+
+if( nn == 0 )
+	return( 1.0 );
+
+if( nn == -1 )
+	return( 1.0/x );
+
+if( x < 0.0 )
+	{
+	asign = -1;
+	x = -x;
+	}
+else
+	asign = 0;
+
+
+if( nn < 0 )
+	{
+	sign = -1;
+	n = -nn;
+	}
+else
+	{
+	sign = 1;
+	n = nn;
+	}
+
+/* Even power will be positive. */
+if( (n & 1) == 0 )
+	asign = 0;
+
+/* Overflow detection */
+
+/* Calculate approximate logarithm of answer */
+s = frexpf( x, &lx );
+e = (lx - 1)*n;
+if( (e == 0) || (e > 64) || (e < -64) )
+	{
+	s = (s - 7.0710678118654752e-1) / (s +  7.0710678118654752e-1);
+	s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
+	}
+else
+	{
+	s = LOGE2F * e;
+	}
+
+if( s > MAXLOGF )
+	{
+	mtherr( "__powif", OVERFLOW );
+	_SET_ERRNO(ERANGE);
+	y = INFINITYF;
+	goto done;
+	}
+
+#if DENORMAL
+if( s < MINLOGF )
+	{
+	y = 0.0;
+	goto done;
+	}
+
+/* Handle tiny denormal answer, but with less accuracy
+ * since roundoff error in 1.0/x will be amplified.
+ * The precise demarcation should be the gradual underflow threshold.
+ */
+if( (s < (-MAXLOGF+2.0)) && (sign < 0) )
+	{
+	x = 1.0/x;
+	sign = -sign;
+	}
+#else
+/* do not produce denormal answer */
+if( s < -MAXLOGF )
+	return(0.0);
+#endif
+
+
+/* First bit of the power */
+if( n & 1 )
+	y = x;
+		
+else
+	y = 1.0;
+
+w = x;
+n >>= 1;
+while( n )
+	{
+	w = w * w;	/* arg to the 2-to-the-kth power */
+	if( n & 1 )	/* if that bit is set, then include in product */
+		y *= w;
+	n >>= 1;
+	}
+
+if( sign < 0 )
+	y = 1.0/y;
+
+done:
+
+if( asign )
+	{
+	/* odd power of negative number */
+	if( y == 0.0 )
+		y = NEGZEROF;
+	else
+		y = -y;
+	}
+return(y);
+}
diff --git a/winsup/mingw/mingwex/math/powil.c b/winsup/mingw/mingwex/math/powil.c
index e900dd881..ec7a2866b 100644
--- a/winsup/mingw/mingwex/math/powil.c
+++ b/winsup/mingw/mingwex/math/powil.c
@@ -1,4 +1,4 @@
-/*							powil.c
+/*							__powil.c
  *
  *	Real raised to integer power, long double precision
  *
@@ -6,10 +6,10 @@
  *
  * SYNOPSIS:
  *
- * long double x, y, powil();
+ * long double x, y, __powil();
  * int n;
  *
- * y = powil( x, n );
+ * y = __powil( x, n );
  *
  *
  *
@@ -36,7 +36,7 @@
  *
  */
 
-/*							powil.c	*/
+/*							__powil.c	*/
 
 /*
 Cephes Math Library Release 2.2:  December, 1990
@@ -66,7 +66,7 @@ long double frexpl();
 #define _SET_ERRNO(x)
 #endif
 
-long double powil( x, nn )
+long double __powil( x, nn )
 long double x;
 int nn;
 {
@@ -126,7 +126,7 @@ else
 
 if( s > MAXLOGL )
 	{
-	mtherr( "powil", OVERFLOW );
+	mtherr( "__powil", OVERFLOW );
 	_SET_ERRNO(ERANGE);
 	y = INFINITYL;
 	goto done;
@@ -134,7 +134,7 @@ if( s > MAXLOGL )
 
 if( s < MINLOGL )
 	{
-	mtherr( "powil", UNDERFLOW );
+	mtherr( "__powil", UNDERFLOW );
 	_SET_ERRNO(ERANGE);
 	return(0.0L);
 	}
diff --git a/winsup/mingw/mingwex/math/powl.c b/winsup/mingw/mingwex/math/powl.c
index a94ede965..f066eeaee 100644
--- a/winsup/mingw/mingwex/math/powl.c
+++ b/winsup/mingw/mingwex/math/powl.c
@@ -382,7 +382,7 @@ static long double w, W, Wa, Wb, ya, yb, u;
 
 #ifdef __MINGW32__
 static __inline__ long double reducl( long double );
-extern long double powil ( long double, int );
+extern long double __powil ( long double, int );
 extern long double powl ( long double x, long double y);
 #else
 #ifdef ANSIPROT
@@ -392,14 +392,14 @@ extern long double frexpl ( long double, int * );
 extern long double ldexpl ( long double, int );
 extern long double polevll ( long double, void *, int );
 extern long double p1evll ( long double, void *, int );
-extern long double powil ( long double, int );
+extern long double __powil ( long double, int );
 extern int isnanl ( long double );
 extern int isfinitel ( long double );
 static long double reducl( long double );
 extern int signbitl ( long double );
 #else
 long double floorl(), fabsl(), frexpl(), ldexpl();
-long double polevll(), p1evll(), powil();
+long double polevll(), p1evll(), __powil();
 static long double reducl();
 int isnanl(), isfinitel(), signbitl();
 #endif
@@ -603,7 +603,7 @@ if( iyflg )
 	w = floorl(x);
 	if( (w == x) && (fabsl(y) < 32768.0) )
 		{
-		w = powil( x, (int) y );
+		w = __powil( x, (int) y );
 		return( w );
 		}
 	}
@@ -764,40 +764,3 @@ t = ldexpl( t, -LNXT );
 return(t);
 }
 
-/*							powil.c
- *
- *	Real raised to integer power, long double precision
- *
- *
- *
- * SYNOPSIS:
- *
- * long double x, y, powil();
- * int n;
- *
- * y = powil( x, n );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns argument x raised to the nth power.
- * The routine efficiently decomposes n as a sum of powers of
- * two. The desired power is a product of two-to-the-kth
- * powers of x.  Thus to compute the 32767 power of x requires
- * 28 multiplications instead of 32767 multiplications.
- *
- *
- *
- * ACCURACY:
- *
- *
- *                      Relative error:
- * arithmetic   x domain   n domain  # trials      peak         rms
- *    IEEE     .001,1000  -1022,1023  50000       4.3e-17     7.8e-18
- *    IEEE        1,2     -1022,1023  20000       3.9e-17     7.6e-18
- *    IEEE     .99,1.01     0,8700    10000       3.6e-16     7.2e-17
- *
- * Returns INFINITY on overflow, zero on underflow.
- *
- */