newlib: update

git-svn-id: svn://kolibrios.org@1906 a494cfbc-eb01-0410-851d-a64ba20cac60
2011-03-11 18:52:24 +00:00
parent a316fa7c9d
commit 2336060a0c
297 changed files with 26930 additions and 2094 deletions
@@ -0,0 +1,40 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+#include <math.h>
+
+float
+acosf (float x)
+{
+  float res;
+
+  /* acosl = atanl (sqrtl(1 - x^2) / x) */
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"		/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"			/* 1 - x^2 */
+	"fsqrt\n\t"			/* sqrtl (1 - x^2) */
+	"fxch	%%st(1)\n\t"
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
+
+double
+acos (double x)
+{
+  double res;
+
+  /* acosl = atanl (sqrtl(1 - x^2) / x) */
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"		/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"			/* 1 - x^2 */
+	"fsqrt\n\t"			/* sqrtl (1 - x^2) */
+	"fxch	%%st(1)\n\t"
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,26 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* acosh(x) = log (x + sqrt(x * x - 1)) */
+double acosh (double x)
+{
+  if (isnan (x)) 
+    return x;
+
+  if (x < 1.0)
+    {
+      errno = EDOM;
+      return nan("");
+    }
+
+  if (x > 0x1p32)
+    /*  Avoid overflow (and unnecessary calculation when
+        sqrt (x * x - 1) == x). GCC optimizes by replacing
+        the long double M_LN2 const with a fldln2 insn.  */ 
+    return __fast_log (x) + 6.9314718055994530941723E-1L;
+
+  /* Since  x >= 1, the arg to log will always be greater than
+     the fyl2xp1 limit (approx 0.29) so just use logl. */ 
+  return __fast_log (x + __fast_sqrt((x + 1.0) * (x - 1.0)));
+}
@@ -0,0 +1,25 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* acosh(x) = log (x + sqrt(x * x - 1)) */
+float acoshf (float x)
+{
+  if (isnan (x)) 
+    return x;
+  if (x < 1.0f)
+    {
+      errno = EDOM;
+      return nan("");
+    }
+
+ if (x > 0x1p32f)
+    /*  Avoid overflow (and unnecessary calculation when
+        sqrt (x * x - 1) == x). GCC optimizes by replacing
+        the long double M_LN2 const with a fldln2 insn.  */ 
+    return __fast_log (x) + 6.9314718055994530941723E-1L;
+
+  /* Since  x >= 1, the arg to log will always be greater than
+     the fyl2xp1 limit (approx 0.29) so just use logl. */ 
+  return __fast_log (x + __fast_sqrt((x + 1.0) * (x - 1.0)));
+}
@@ -0,0 +1,27 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* acosh(x) = log (x + sqrt(x * x - 1)) */
+long double acoshl (long double x)
+{
+  if (isnan (x)) 
+    return x;
+
+  if (x < 1.0L)
+    {
+      errno = EDOM;
+      return nanl("");
+    }
+  if (x > 0x1p32L)
+    /* Avoid overflow (and unnecessary calculation when
+       sqrt (x * x - 1) == x).
+       The M_LN2 define doesn't have enough precison for
+       long double so use this one. GCC optimizes by replacing
+       the const with a fldln2 insn. */
+    return __fast_logl (x) + 6.9314718055994530941723E-1L;
+
+   /* Since  x >= 1, the arg to log will always be greater than
+      the fyl2xp1 limit (approx 0.29) so just use logl. */ 
+   return __fast_logl (x + __fast_sqrtl((x + 1.0L) * (x - 1.0L)));
+}
@@ -0,0 +1,25 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+#include <math.h>
+
+long double
+acosl (long double x)
+{
+  long double res;
+
+  /* acosl = atanl (sqrtl(1 - x^2) / x) */
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"		/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"			/* 1 - x^2 */
+	"fsqrt\n\t"			/* sqrtl (1 - x^2) */
+	"fxch	%%st(1)\n\t"
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,34 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+/* asin = atan (x / sqrt(1 - x^2)) */
+
+float asinf (float x)
+{
+  float res;
+
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"			/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"				/* 1 - x^2 */
+	"fsqrt\n\t"				/* sqrt (1 - x^2) */
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
+
+double asin (double x)
+{
+  double res;
+
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"			/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"				/* 1 - x^2 */
+	"fsqrt\n\t"				/* sqrt (1 - x^2) */
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,28 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+ /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
+double asinh(double x)
+{
+  double z;
+  if (!isfinite (x))
+    return x;
+  z = fabs (x);
+
+  /* Avoid setting FPU underflow exception flag in x * x. */
+#if 0
+  if ( z < 0x1p-32)
+    return x;
+#endif
+
+  /* Use log1p to avoid cancellation with small x. Put
+     x * x in denom, so overflow is harmless. 
+     asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
+              = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
+
+  z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
+
+  return ( x > 0.0 ? z : -z);
+}
+
@@ -0,0 +1,28 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+ /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
+float asinhf(float x)
+{
+  float z;
+  if (!isfinite (x))
+    return x;
+  z = fabsf (x);
+
+  /* Avoid setting FPU underflow exception flag in x * x. */
+#if 0
+  if ( z < 0x1p-32)
+    return x;
+#endif
+
+
+  /* Use log1p to avoid cancellation with small x. Put
+     x * x in denom, so overflow is harmless. 
+     asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
+              = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
+
+  z = __fast_log1p (z + z * z / (__fast_sqrt (z * z + 1.0) + 1.0));
+
+  return ( x > 0.0 ? z : -z);
+}
@@ -0,0 +1,28 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+ /* asinh(x) = copysign(log(fabs(x) + sqrt(x * x + 1.0)), x) */
+long double asinhl(long double x)
+{
+  long double z;
+  if (!isfinite (x))
+    return x;
+
+  z = fabsl (x);
+
+  /* Avoid setting FPU underflow exception flag in x * x. */
+#if 0
+  if ( z < 0x1p-32)
+    return x;
+#endif
+
+  /* Use log1p to avoid cancellation with small x. Put
+     x * x in denom, so overflow is harmless. 
+     asinh(x) = log1p (x + sqrt (x * x + 1.0) - 1.0)
+              = log1p (x + x * x / (sqrt (x * x + 1.0) + 1.0))  */
+
+  z = __fast_log1pl (z + z * z / (__fast_sqrtl (z * z + 1.0L) + 1.0L));
+
+  return ( x > 0.0 ? z : -z);
+}
@@ -0,0 +1,21 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ * Adapted for long double type by Danny Smith <dannysmith@users.sourceforge.net>.
+ */
+
+/* asin = atan (x / sqrt(1 - x^2)) */
+
+long double asinl (long double x)
+{
+  long double res;
+
+  asm (	"fld	%%st\n\t"
+	"fmul	%%st(0)\n\t"			/* x^2 */
+	"fld1\n\t"
+	"fsubp\n\t"				/* 1 - x^2 */
+	"fsqrt\n\t"				/* sqrt (1 - x^2) */
+	"fpatan"
+	: "=t" (res) : "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,15 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ */
+
+#include <math.h>
+
+float
+atan2f (float y, float x)
+{
+  float res;
+  asm ("fpatan" : "=t" (res) : "u" (y), "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,16 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+#include <math.h>
+
+long double
+atan2l (long double y, long double x)
+{
+  long double res;
+  asm ("fpatan" : "=t" (res) : "u" (y), "0" (x) : "st(1)");
+  return res;
+}
@@ -0,0 +1,28 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ */
+
+#include <math.h>
+
+float
+atanf (float x)
+{
+  float res;
+
+  asm ("fld1\n\t"
+       "fpatan" : "=t" (res) : "0" (x));
+  return res;
+}
+
+double
+atan (double x)
+{
+  double res;
+
+  asm ("fld1 \n\t"
+       "fpatan" : "=t" (res) : "0" (x));
+  return res;
+}
+
@@ -0,0 +1,31 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */
+
+double atanh(double x)
+{
+  double z;
+  if isnan (x)
+    return x;
+  z = fabs (x);
+  if (z == 1.0)
+    {
+      errno  = ERANGE;
+      return (x > 0 ? INFINITY : -INFINITY);
+    }
+  if (z > 1.0)
+    {
+      errno = EDOM;
+      return nan("");
+    }
+  /* Rearrange formula to avoid precision loss for small x.
+
+  atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
+	   = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
+           = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) 
+           = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
+  z = 0.5 * __fast_log1p ((z + z) / (1.0 - z));
+  return x >= 0 ? z : -z;
+}
@@ -0,0 +1,30 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */
+float atanhf (float x)
+{
+  float z;
+  if isnan (x)
+    return x;
+  z = fabsf (x);
+  if (z == 1.0)
+    {
+      errno  = ERANGE;
+      return (x > 0 ? INFINITY : -INFINITY);
+    }
+  if ( z > 1.0)
+    {
+      errno = EDOM;
+      return nanf("");
+    }
+  /* Rearrange formula to avoid precision loss for small x.
+
+  atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
+	   = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
+           = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) 
+           = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
+  z = 0.5 * __fast_log1p ((z + z) / (1.0 - z));
+  return x >= 0 ? z : -z;
+}
@@ -0,0 +1,29 @@
+#include <math.h>
+#include <errno.h>
+#include "fastmath.h"
+
+/* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */
+long double atanhl (long double x)
+{
+  long double z;
+  if isnan (x)
+    return x;
+  z = fabsl (x);
+  if (z == 1.0L)
+    {
+      errno  = ERANGE;
+      return (x > 0 ? INFINITY : -INFINITY);
+    }
+  if ( z > 1.0L)
+    {
+      errno = EDOM;
+      return nanl("");
+    }
+  /* Rearrange formula to avoid precision loss for small x.
+  atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x))
+ 	   = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0)
+           = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) 
+           = 0.5 * log1p ((2.0 * x ) / (1.0 - x))  */
+  z = 0.5L * __fast_log1pl ((z + z) / (1.0L - z));
+  return x >= 0 ? z : -z;
+}
@@ -0,0 +1,19 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+#include <math.h>
+
+long double
+atanl (long double x)
+{
+  long double res;
+
+  asm ("fld1\n\t"
+       "fpatan"
+       : "=t" (res) : "0" (x));
+  return res;
+}
@@ -0,0 +1,162 @@
+/*							cbrt.c
+ *
+ *	Cube root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, cbrt();
+ *
+ * y = cbrt( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the cube root of the argument, which may be negative.
+ *
+ * Range reduction involves determining the power of 2 of
+ * the argument.  A polynomial of degree 2 applied to the
+ * mantissa, and multiplication by the cube root of 1, 2, or 4
+ * approximates the root to within about 0.1%.  Then Newton's
+ * iteration is used three times to converge to an accurate
+ * result.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC        -10,10     200000      1.8e-17     6.2e-18
+ *    IEEE       0,1e308     30000      1.5e-16     5.0e-17
+ *
+ */
+/*							cbrt.c  */
+
+/*
+Cephes Math Library Release 2.2:  January, 1991
+Copyright 1984, 1991 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/*
+  Modified for mingwex.a 
+  2002-07-01  Danny Smith  <dannysmith@users.sourceforge.net>
+ */
+#ifdef __MINGW32__
+#include <math.h>
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+#endif
+
+
+static const double CBRT2  = 1.2599210498948731647672;
+static const double CBRT4  = 1.5874010519681994747517;
+static const double CBRT2I = 0.79370052598409973737585;
+static const double CBRT4I = 0.62996052494743658238361;
+
+#ifndef __MINGW32__
+#ifdef ANSIPROT
+extern double frexp ( double, int * );
+extern double ldexp ( double, int );
+extern int isnan ( double );
+extern int isfinite ( double );
+#else
+double frexp(), ldexp();
+int isnan(), isfinite();
+#endif
+#endif
+
+double cbrt(x)
+double x;
+{
+int e, rem, sign;
+double z;
+
+#ifdef __MINGW32__
+if (!isfinite (x) || x == 0 )
+  return x;
+#else
+
+#ifdef NANS
+if( isnan(x) )
+  return x;
+#endif
+#ifdef INFINITIES
+if( !isfinite(x) )
+  return x;
+#endif
+if( x == 0 )
+	return( x );
+
+#endif /* __MINGW32__ */
+
+if( x > 0 )
+	sign = 1;
+else
+	{
+	sign = -1;
+	x = -x;
+	}
+
+z = x;
+/* extract power of 2, leaving
+ * mantissa between 0.5 and 1
+ */
+x = frexp( x, &e );
+
+/* Approximate cube root of number between .5 and 1,
+ * peak relative error = 9.2e-6
+ */
+x = (((-1.3466110473359520655053e-1  * x
+      + 5.4664601366395524503440e-1) * x
+      - 9.5438224771509446525043e-1) * x
+      + 1.1399983354717293273738e0 ) * x
+      + 4.0238979564544752126924e-1;
+
+/* exponent divided by 3 */
+if( e >= 0 )
+	{
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x *= CBRT2;
+	else if( rem == 2 )
+		x *= CBRT4;
+	}
+
+
+/* argument less than 1 */
+
+else
+	{
+	e = -e;
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x *= CBRT2I;
+	else if( rem == 2 )
+		x *= CBRT4I;
+	e = -e;
+	}
+
+/* multiply by power of 2 */
+x = ldexp( x, e );
+
+/* Newton iteration */
+x -= ( x - (z/(x*x)) )*0.33333333333333333333;
+#ifdef DEC
+x -= ( x - (z/(x*x)) )/3.0;
+#else
+x -= ( x - (z/(x*x)) )*0.33333333333333333333;
+#endif
+
+if( sign < 0 )
+	x = -x;
+return(x);
+}
@@ -0,0 +1,147 @@
+/*							cbrtf.c
+ *
+ *	Cube root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, cbrtf();
+ *
+ * y = cbrtf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the cube root of the argument, which may be negative.
+ *
+ * Range reduction involves determining the power of 2 of
+ * the argument.  A polynomial of degree 2 applied to the
+ * mantissa, and multiplication by the cube root of 1, 2, or 4
+ * approximates the root to within about 0.1%.  Then Newton's
+ * iteration is used to converge to an accurate result.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,1e38      100000      7.6e-8      2.7e-8
+ *
+ */
+/*							cbrt.c  */
+
+/*
+Cephes Math Library Release 2.2:  June, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/*
+  Modified for mingwex.a 
+  2002-07-01  Danny Smith  <dannysmith@users.sourceforge.net>
+ */
+#ifdef __MINGW32__
+#include <math.h>
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+#endif
+
+static const float CBRT2 = 1.25992104989487316477;
+static const float CBRT4 = 1.58740105196819947475;
+
+#ifndef __MINGW32__
+#ifdef ANSIC
+float frexpf(float, int *), ldexpf(float, int);
+
+float cbrtf( float xx )
+#else
+float frexpf(), ldexpf();
+
+float cbrtf(xx)
+double xx;
+#endif
+{
+int e, rem, sign;
+float x, z;
+
+x = xx;
+
+#else /* __MINGW32__ */
+float cbrtf (float x)
+{
+int e, rem, sign;
+float z;
+#endif /* __MINGW32__ */
+
+#ifdef __MINGW32__
+if (!isfinite (x) || x == 0.0F )
+  return x;
+#else
+if( x == 0 )
+	return( 0.0 );
+#endif
+if( x > 0 )
+	sign = 1;
+else
+	{
+	sign = -1;
+	x = -x;
+	}
+
+z = x;
+/* extract power of 2, leaving
+ * mantissa between 0.5 and 1
+ */
+x = frexpf( x, &e );
+
+/* Approximate cube root of number between .5 and 1,
+ * peak relative error = 9.2e-6
+ */
+x = (((-0.13466110473359520655053  * x
+      + 0.54664601366395524503440 ) * x
+      - 0.95438224771509446525043 ) * x
+      + 1.1399983354717293273738  ) * x
+      + 0.40238979564544752126924;
+
+/* exponent divided by 3 */
+if( e >= 0 )
+	{
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x *= CBRT2;
+	else if( rem == 2 )
+		x *= CBRT4;
+	}
+
+
+/* argument less than 1 */
+
+else
+	{
+	e = -e;
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x /= CBRT2;
+	else if( rem == 2 )
+		x /= CBRT4;
+	e = -e;
+	}
+
+/* multiply by power of 2 */
+x = ldexpf( x, e );
+
+/* Newton iteration */
+x -= ( x - (z/(x*x)) ) * 0.333333333333;
+
+if( sign < 0 )
+	x = -x;
+return(x);
+}
@@ -0,0 +1,161 @@
+/*							cbrtl.c
+ *
+ *	Cube root, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, cbrtl();
+ *
+ * y = cbrtl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the cube root of the argument, which may be negative.
+ *
+ * Range reduction involves determining the power of 2 of
+ * the argument.  A polynomial of degree 2 applied to the
+ * mantissa, and multiplication by the cube root of 1, 2, or 4
+ * approximates the root to within about 0.1%.  Then Newton's
+ * iteration is used three times to converge to an accurate
+ * result.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     .125,8        80000      7.0e-20     2.2e-20
+ *    IEEE    exp(+-707)    100000      7.0e-20     2.4e-20
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.2: January, 1991
+Copyright 1984, 1991 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/*
+  Modified for mingwex.a 
+  2002-07-01  Danny Smith  <dannysmith@users.sourceforge.net>
+ */
+#ifdef __MINGW32__
+#include "cephes_mconf.h"
+#else
+#include "mconf.h"
+#endif
+
+static const long double CBRT2  = 1.2599210498948731647672L;
+static const long double CBRT4  = 1.5874010519681994747517L;
+static const long double CBRT2I = 0.79370052598409973737585L;
+static const long double CBRT4I = 0.62996052494743658238361L;
+
+#ifndef __MINGW32__
+
+#ifdef ANSIPROT
+extern long double frexpl ( long double, int * );
+extern long double ldexpl ( long double, int );
+extern int isnanl ( long double );
+#else
+long double frexpl(), ldexpl();
+extern int isnanl();
+#endif
+
+#ifdef INFINITIES
+extern long double INFINITYL;
+#endif
+
+#endif /* __MINGW32__ */
+
+long double cbrtl(x)
+long double x;
+{
+int e, rem, sign;
+long double z;
+
+#ifdef __MINGW32__
+if (!isfinite (x) || x == 0.0L)
+	return(x);
+#else     
+
+#ifdef NANS
+if(isnanl(x))
+	return(x);
+#endif
+#ifdef INFINITIES
+if( x == INFINITYL)
+	return(x);
+if( x == -INFINITYL)
+	return(x);
+#endif
+if( x == 0 )
+	return( x );
+
+#endif /* __MINGW32__ */
+
+if( x > 0 )
+	sign = 1;
+else
+	{
+	sign = -1;
+	x = -x;
+	}
+
+z = x;
+/* extract power of 2, leaving
+ * mantissa between 0.5 and 1
+ */
+x = frexpl( x, &e );
+
+/* Approximate cube root of number between .5 and 1,
+ * peak relative error = 1.2e-6
+ */
+x = (((( 1.3584464340920900529734e-1L * x
+       - 6.3986917220457538402318e-1L) * x
+       + 1.2875551670318751538055e0L) * x
+       - 1.4897083391357284957891e0L) * x
+       + 1.3304961236013647092521e0L) * x
+       + 3.7568280825958912391243e-1L;
+
+/* exponent divided by 3 */
+if( e >= 0 )
+	{
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x *= CBRT2;
+	else if( rem == 2 )
+		x *= CBRT4;
+	}
+else
+	{ /* argument less than 1 */
+	e = -e;
+	rem = e;
+	e /= 3;
+	rem -= 3*e;
+	if( rem == 1 )
+		x *= CBRT2I;
+	else if( rem == 2 )
+		x *= CBRT4I;
+	e = -e;
+	}
+
+/* multiply by power of 2 */
+x = ldexpl( x, e );
+
+/* Newton iteration */
+
+x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
+x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
+
+if( sign < 0 )
+	x = -x;
+return(x);
+}
@@ -0,0 +1,31 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+	.file	"ceil.s"
+	.text
+	.align 4
+.globl _ceil
+	.def	_ceil;	.scl	2;	.type	32;	.endef
+_ceil:
+	fldl	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x0800,%edx		/* round towards +oo */
+	orl	4(%esp),%edx
+	andl	$0xfbff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,31 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+	.file	"ceilf.S"
+	.text
+	.align 4
+.globl _ceilf
+	.def	_ceilf;	.scl	2;	.type	32;	.endef
+_ceilf:
+	flds	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x0800,%edx		/* round towards +oo */
+	orl	4(%esp),%edx
+	andl	$0xfbff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,33 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ */
+
+
+	.file	"ceill.S"
+	.text
+	.align 4
+.globl _ceill
+	.def	_ceill;	.scl	2;	.type	32;	.endef
+_ceill:
+	fldt	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x0800,%edx		/* round towards +oo */
+	orl	4(%esp),%edx
+	andl	$0xfbff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,395 @@
+#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,
+
+//#define _CEPHES_USE_ERRNO
+
+#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)
+#undef NAN
+#undef INFINITY
+#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
+#define INFINITY __builtin_huge_val()
+#define NAN __builtin_nan("")
+#else
+extern double __INF;
+#define INFINITY (__INF)
+extern double __QNAN;
+#define NAN (__QNAN)
+#endif
+
+/*long double*/
+#define MAXNUML 1.189731495357231765021263853E4932L
+#define MAXLOGL	1.1356523406294143949492E4L
+#define MINLOGL	-1.13994985314888605586758E4L
+#define LOGE2L	6.9314718055994530941723E-1L
+#define LOG2EL	1.4426950408889634073599E0L
+#define PIL	3.1415926535897932384626L
+#define PIO2L	1.5707963267948966192313L
+#define PIO4L	7.8539816339744830961566E-1L
+
+#define isfinitel isfinite
+#define isinfl isinf
+#define isnanl isnan
+#define signbitl signbit
+
+#define NEGZEROL (-0.0L)
+
+#undef NANL
+#undef INFINITYL
+#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
+#define INFINITYL __builtin_huge_vall()
+#define NANL __builtin_nanl("")
+#else
+extern long double __INFL;
+#define INFINITYL (__INFL)
+extern long double __QNANL;
+#define NANL (__QNANL)
+#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)
+
+#undef NANF
+#undef INFINITYF
+#if (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ > 2))
+#define INFINITYF __builtin_huge_valf()
+#define NANF __builtin_nanf("")
+#else
+extern float __INFF;
+#define INFINITYF (__INFF)
+extern float __QNANF;
+#define NANF (__QNANF)
+#endif
+
+
+/* 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
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+/*							polevll.c
+ *							p1evll.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * long double x, y, coef[N+1], polevl[];
+ *
+ * y = polevll( 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 p1evll() assumes that coef[N] = 1.0 and is
+ * omitted from the array.  Its calling arguments are
+ * otherwise the same as polevll().
+ *
+ *
+ * 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__ long double polevll( x, p, n )
+long double x;
+const void *p;
+int n;
+{
+register long double y;
+register long double *P = (long 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__ long double p1evll( x, p, n )
+long double x;
+const void *p;
+int n;
+{
+register long double y;
+register long double *P = (long double *)p;
+
+n -= 1;
+y = x + *P++;
+do
+	{
+	y = y * x + *P++;
+	}
+while( --n );
+return( y );
+}
+
+/* Float version */
+
+/*							polevlf.c
+ *							p1evlf.c
+ *
+ *	Evaluate polynomial
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int N;
+ * float x, y, coef[N+1], polevlf[];
+ *
+ * y = polevlf( 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.
+ *
+ */
+
+/*
+Cephes Math Library Release 2.1:  December, 1988
+Copyright 1984, 1987, 1988 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+static __inline__ float polevlf(float x, const float* coef, int N )
+{
+float ans;
+float *p;
+int i;
+
+p = (float*)coef;
+ans = *p++;
+
+/*
+for( i=0; i<N; i++ )
+	ans = ans * x  +  *p++;
+*/
+
+i = N;
+do
+	ans = ans * x  +  *p++;
+while( --i );
+
+return( ans );
+}
+
+/*							p1evl()	*/
+/*                                          N
+ * Evaluate polynomial when coefficient of x  is 1.0.
+ * Otherwise same as polevl.
+ */
+
+static __inline__ float p1evlf( float x, const float *coef, int N )
+{
+float ans;
+float *p;
+int i;
+
+p = (float*)coef;
+ans = x + *p++;
+i = N-1;
+
+do
+	ans = ans * x  + *p++;
+while( --i );
+
+return( ans );
+}
@@ -0,0 +1,19 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+	.file	"copysign.S"
+	.text
+	.align 4
+.globl _copysign
+	.def	_copysign;	.scl	2;	.type	32;	.endef
+_copysign:
+	movl	16(%esp),%edx
+	movl	8(%esp),%eax
+	andl	$0x80000000,%edx
+	andl	$0x7fffffff,%eax
+	orl	%edx,%eax
+	movl	%eax,8(%esp)
+	fldl	4(%esp)
+	ret
@@ -0,0 +1,19 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+	.file	"copysignf.S"
+	.text
+	.align 4
+.globl _copysignf
+	.def	_copysignf;	.scl	2;	.type	32;	.endef
+_copysignf:
+	movl	8(%esp),%edx
+	movl	4(%esp),%eax
+	andl	$0x80000000,%edx
+	andl	$0x7fffffff,%eax
+	orl	%edx,%eax
+	movl	%eax,4(%esp)
+	flds	4(%esp)
+	ret
@@ -0,0 +1,20 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ * Public domain.
+ */
+
+	.file	"copysignl.S"
+	.text
+	.align 4
+.globl _copysignl
+	.def	_copysignl;	.scl	2;	.type	32;	.endef
+_copysignl:
+	movl	24(%esp),%edx
+	movl	12(%esp),%eax
+	andl	$0x8000,%edx
+	andl	$0x7fff,%eax
+	orl	%edx,%eax
+	movl	%eax,12(%esp)
+	fldt	4(%esp)
+	ret
@@ -0,0 +1,29 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Removed glibc header dependancy by Danny Smith
+ *  <dannysmith@users.sourceforge.net>
+ */
+	.file	"cos.s"
+	.text
+	.align 4
+.globl _cos
+	.def	_cos;	.scl	2;	.type	32;	.endef
+_cos:
+	fldl	4(%esp)
+	fcos
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	1f
+	ret
+1:	fldpi
+	fadd	%st(0)
+	fxch	%st(1)
+2:	fprem1
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	2b
+	fstp	%st(1)
+	fcos
+	ret
@@ -0,0 +1,29 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Removed glibc header dependancy by Danny Smith
+ *  <dannysmith@users.sourceforge.net>
+ */
+	.file	"cosf.S"
+	.text
+	.align 4
+.globl _cosl
+	.def	_cosf;	.scl	2;	.type	32;	.endef
+_cosf:
+	flds	4(%esp)
+	fcos
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	1f
+	ret
+1:	fldpi
+	fadd	%st(0)
+	fxch	%st(1)
+2:	fprem1
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	2b
+	fstp	%st(1)
+	fcos
+	ret
@@ -0,0 +1,3 @@
+#include <math.h>
+float coshf (float x)
+  {return (float) cosh (x);}
@@ -0,0 +1,110 @@
+/*							coshl.c
+ *
+ *	Hyperbolic cosine, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, coshl();
+ *
+ * y = coshl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns hyperbolic cosine of argument in the range MINLOGL to
+ * MAXLOGL.
+ *
+ * cosh(x)  =  ( exp(x) + exp(-x) )/2.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     +-10000      30000       1.1e-19     2.8e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition              value returned
+ * cosh overflow    |x| > MAXLOGL+LOGE2L      INFINITYL
+ *
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.7:  May, 1998
+Copyright 1985, 1991, 1998 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"
+#endif
+
+#ifndef _SET_ERRNO
+#define _SET_ERRNO(x)
+#endif
+
+
+#ifndef __MINGW32__
+extern long double MAXLOGL, MAXNUML, LOGE2L;
+#ifdef ANSIPROT
+extern long double expl ( long double );
+extern int isnanl ( long double );
+#else
+long double expl(), isnanl();
+#endif
+#ifdef INFINITIES
+extern long double INFINITYL;
+#endif
+#ifdef NANS
+extern long double NANL;
+#endif
+#endif /* __MINGW32__ */
+
+long double coshl(x)
+long double x;
+{
+long double y;
+
+#ifdef NANS
+if( isnanl(x) )
+	{
+	_SET_ERRNO(EDOM);
+	return(x);
+  	}
+#endif
+if( x < 0 )
+	x = -x;
+if( x > (MAXLOGL + LOGE2L) )
+	{
+	mtherr( "coshl", OVERFLOW );
+	_SET_ERRNO(ERANGE);
+#ifdef INFINITIES
+	return( INFINITYL );
+#else
+	return( MAXNUML );
+#endif
+	}	
+if( x >= (MAXLOGL - LOGE2L) )
+	{
+	y = expl(0.5L * x);
+	y = (0.5L * y) * y;
+	return(y);
+	}
+y = expl(x);
+y = 0.5L * (y + 1.0L / y);
+return( y );
+}
@@ -0,0 +1,30 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ * Removed glibc header dependancy by Danny Smith
+ *  <dannysmith@users.sourceforge.net>
+ */
+	.file	"cosl.S"
+	.text
+	.align 4
+.globl _cosl
+	.def	_cosl;	.scl	2;	.type	32;	.endef
+_cosl:
+	fldt	4(%esp)
+	fcos
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	1f
+	ret
+1:	fldpi
+	fadd	%st(0)
+	fxch	%st(1)
+2:	fprem1
+	fnstsw	%ax
+	testl	$0x400,%eax
+	jnz	2b
+	fstp	%st(1)
+	fcos
+	ret
@@ -0,0 +1,265 @@
+/* Software floating-point emulation.
+   Definitions for IEEE Double Precision
+   Copyright (C) 1997, 1998, 1999, 2006, 2007, 2008, 2009
+   Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   In addition to the permissions in the GNU Lesser General Public
+   License, the Free Software Foundation gives you unlimited
+   permission to link the compiled version of this file into
+   combinations with other programs, and to distribute those
+   combinations without any restriction coming from the use of this
+   file.  (The Lesser General Public License restrictions do apply in
+   other respects; for example, they cover modification of the file,
+   and distribution when not linked into a combine executable.)
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, 51 Franklin Street, Fifth Floor, Boston,
+   MA 02110-1301, USA.  */
+
+#if _FP_W_TYPE_SIZE < 32
+#error "Here's a nickel kid.  Go buy yourself a real computer."
+#endif
+
+#if _FP_W_TYPE_SIZE < 64
+#define _FP_FRACTBITS_D		(2 * _FP_W_TYPE_SIZE)
+#else
+#define _FP_FRACTBITS_D		_FP_W_TYPE_SIZE
+#endif
+
+#define _FP_FRACBITS_D		53
+#define _FP_FRACXBITS_D		(_FP_FRACTBITS_D - _FP_FRACBITS_D)
+#define _FP_WFRACBITS_D		(_FP_WORKBITS + _FP_FRACBITS_D)
+#define _FP_WFRACXBITS_D	(_FP_FRACTBITS_D - _FP_WFRACBITS_D)
+#define _FP_EXPBITS_D		11
+#define _FP_EXPBIAS_D		1023
+#define _FP_EXPMAX_D		2047
+
+#define _FP_QNANBIT_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-2) % _FP_W_TYPE_SIZE)
+#define _FP_QNANBIT_SH_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-2+_FP_WORKBITS) % _FP_W_TYPE_SIZE)
+#define _FP_IMPLBIT_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-1) % _FP_W_TYPE_SIZE)
+#define _FP_IMPLBIT_SH_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-1+_FP_WORKBITS) % _FP_W_TYPE_SIZE)
+#define _FP_OVERFLOW_D		\
+	((_FP_W_TYPE)1 << _FP_WFRACBITS_D % _FP_W_TYPE_SIZE)
+
+typedef float DFtype __attribute__((mode(DF)));
+
+#if _FP_W_TYPE_SIZE < 64
+
+union _FP_UNION_D
+{
+  DFtype flt;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign  : 1;
+    unsigned exp   : _FP_EXPBITS_D;
+    unsigned frac1 : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0) - _FP_W_TYPE_SIZE;
+    unsigned frac0 : _FP_W_TYPE_SIZE;
+#else
+    unsigned frac0 : _FP_W_TYPE_SIZE;
+    unsigned frac1 : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0) - _FP_W_TYPE_SIZE;
+    unsigned exp   : _FP_EXPBITS_D;
+    unsigned sign  : 1;
+#endif
+  } bits __attribute__((packed));
+};
+
+#define FP_DECL_D(X)		_FP_DECL(2,X)
+#define FP_UNPACK_RAW_D(X,val)	_FP_UNPACK_RAW_2(D,X,val)
+#define FP_UNPACK_RAW_DP(X,val)	_FP_UNPACK_RAW_2_P(D,X,val)
+#define FP_PACK_RAW_D(val,X)	_FP_PACK_RAW_2(D,val,X)
+#define FP_PACK_RAW_DP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(D,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_D(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2(D,X,val);		\
+    _FP_UNPACK_CANONICAL(D,2,X);	\
+  } while (0)
+
+#define FP_UNPACK_DP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2_P(D,X,val);	\
+    _FP_UNPACK_CANONICAL(D,2,X);	\
+  } while (0)
+
+#define FP_UNPACK_SEMIRAW_D(X,val)	\
+  do {					\
+    _FP_UNPACK_RAW_2(D,X,val);		\
+    _FP_UNPACK_SEMIRAW(D,2,X);		\
+  } while (0)
+
+#define FP_UNPACK_SEMIRAW_DP(X,val)	\
+  do {					\
+    _FP_UNPACK_RAW_2_P(D,X,val);	\
+    _FP_UNPACK_SEMIRAW(D,2,X);		\
+  } while (0)
+
+#define FP_PACK_D(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,2,X);		\
+    _FP_PACK_RAW_2(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_DP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,2,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(D,val,X);	\
+  } while (0)
+
+#define FP_PACK_SEMIRAW_D(val,X)	\
+  do {					\
+    _FP_PACK_SEMIRAW(D,2,X);		\
+    _FP_PACK_RAW_2(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_SEMIRAW_DP(val,X)	\
+  do {					\
+    _FP_PACK_SEMIRAW(D,2,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(D,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_D(X)		_FP_ISSIGNAN(D,2,X)
+#define FP_NEG_D(R,X)			_FP_NEG(D,2,R,X)
+#define FP_ADD_D(R,X,Y)			_FP_ADD(D,2,R,X,Y)
+#define FP_SUB_D(R,X,Y)			_FP_SUB(D,2,R,X,Y)
+#define FP_MUL_D(R,X,Y)			_FP_MUL(D,2,R,X,Y)
+#define FP_DIV_D(R,X,Y)			_FP_DIV(D,2,R,X,Y)
+#define FP_SQRT_D(R,X)			_FP_SQRT(D,2,R,X)
+#define _FP_SQRT_MEAT_D(R,S,T,X,Q)	_FP_SQRT_MEAT_2(R,S,T,X,Q)
+
+#define FP_CMP_D(r,X,Y,un)	_FP_CMP(D,2,r,X,Y,un)
+#define FP_CMP_EQ_D(r,X,Y)	_FP_CMP_EQ(D,2,r,X,Y)
+#define FP_CMP_UNORD_D(r,X,Y)	_FP_CMP_UNORD(D,2,r,X,Y)
+
+#define FP_TO_INT_D(r,X,rsz,rsg)	_FP_TO_INT(D,2,r,X,rsz,rsg)
+#define FP_FROM_INT_D(X,r,rs,rt)	_FP_FROM_INT(D,2,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_D(X)	_FP_FRAC_HIGH_2(X)
+#define _FP_FRAC_HIGH_RAW_D(X)	_FP_FRAC_HIGH_2(X)
+
+#else
+
+union _FP_UNION_D
+{
+  DFtype flt;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign   : 1;
+    unsigned exp    : _FP_EXPBITS_D;
+    _FP_W_TYPE frac : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0);
+#else
+    _FP_W_TYPE frac : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0);
+    unsigned exp    : _FP_EXPBITS_D;
+    unsigned sign   : 1;
+#endif
+  } bits __attribute__((packed));
+};
+
+#define FP_DECL_D(X)		_FP_DECL(1,X)
+#define FP_UNPACK_RAW_D(X,val)	_FP_UNPACK_RAW_1(D,X,val)
+#define FP_UNPACK_RAW_DP(X,val)	_FP_UNPACK_RAW_1_P(D,X,val)
+#define FP_PACK_RAW_D(val,X)	_FP_PACK_RAW_1(D,val,X)
+#define FP_PACK_RAW_DP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(D,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_D(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1(D,X,val);		\
+    _FP_UNPACK_CANONICAL(D,1,X);	\
+  } while (0)
+
+#define FP_UNPACK_DP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1_P(D,X,val);	\
+    _FP_UNPACK_CANONICAL(D,1,X);	\
+  } while (0)
+
+#define FP_UNPACK_SEMIRAW_D(X,val)	\
+  do {					\
+    _FP_UNPACK_RAW_1(D,X,val);		\
+    _FP_UNPACK_SEMIRAW(D,1,X);		\
+  } while (0)
+
+#define FP_UNPACK_SEMIRAW_DP(X,val)	\
+  do {					\
+    _FP_UNPACK_RAW_1_P(D,X,val);	\
+    _FP_UNPACK_SEMIRAW(D,1,X);		\
+  } while (0)
+
+#define FP_PACK_D(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,1,X);		\
+    _FP_PACK_RAW_1(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_DP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,1,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(D,val,X);	\
+  } while (0)
+
+#define FP_PACK_SEMIRAW_D(val,X)	\
+  do {					\
+    _FP_PACK_SEMIRAW(D,1,X);		\
+    _FP_PACK_RAW_1(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_SEMIRAW_DP(val,X)	\
+  do {					\
+    _FP_PACK_SEMIRAW(D,1,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(D,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_D(X)		_FP_ISSIGNAN(D,1,X)
+#define FP_NEG_D(R,X)			_FP_NEG(D,1,R,X)
+#define FP_ADD_D(R,X,Y)			_FP_ADD(D,1,R,X,Y)
+#define FP_SUB_D(R,X,Y)			_FP_SUB(D,1,R,X,Y)
+#define FP_MUL_D(R,X,Y)			_FP_MUL(D,1,R,X,Y)
+#define FP_DIV_D(R,X,Y)			_FP_DIV(D,1,R,X,Y)
+#define FP_SQRT_D(R,X)			_FP_SQRT(D,1,R,X)
+#define _FP_SQRT_MEAT_D(R,S,T,X,Q)	_FP_SQRT_MEAT_1(R,S,T,X,Q)
+
+/* The implementation of _FP_MUL_D and _FP_DIV_D should be chosen by
+   the target machine.  */
+
+#define FP_CMP_D(r,X,Y,un)	_FP_CMP(D,1,r,X,Y,un)
+#define FP_CMP_EQ_D(r,X,Y)	_FP_CMP_EQ(D,1,r,X,Y)
+#define FP_CMP_UNORD_D(r,X,Y)	_FP_CMP_UNORD(D,1,r,X,Y)
+
+#define FP_TO_INT_D(r,X,rsz,rsg)	_FP_TO_INT(D,1,r,X,rsz,rsg)
+#define FP_FROM_INT_D(X,r,rs,rt)	_FP_FROM_INT(D,1,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_D(X)	_FP_FRAC_HIGH_1(X)
+#define _FP_FRAC_HIGH_RAW_D(X)	_FP_FRAC_HIGH_1(X)
+
+#endif /* W_TYPE_SIZE < 64 */
@@ -0,0 +1,131 @@
+
+/* @(#)e_atan2.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+/* __ieee754_atan2(y,x)
+ * Method :
+ *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+ *	2. Reduce x to positive by (if x and y are unexceptional): 
+ *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
+ *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+ *
+ * Special cases:
+ *
+ *	ATAN2((anything), NaN ) is NaN;
+ *	ATAN2(NAN , (anything) ) is NaN;
+ *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
+ *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ *	ATAN2(+-INF,+INF ) is +-pi/4 ;
+ *	ATAN2(+-INF,-INF ) is +-3pi/4;
+ *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following 
+ * constants. The decimal values may be used, provided that the 
+ * compiler will convert from decimal to binary accurately enough 
+ * to produce the hexadecimal values shown.
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static const double 
+#else
+static double 
+#endif
+tiny  = 1.0e-300,
+zero  = 0.0,
+pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
+pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
+pi      = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
+pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
+
+#ifdef __STDC__
+	double __ieee754_atan2(double y, double x)
+#else
+	double __ieee754_atan2(y,x)
+	double  y,x;
+#endif
+{  
+	double z;
+	__int32_t k,m,hx,hy,ix,iy;
+	__uint32_t lx,ly;
+
+	EXTRACT_WORDS(hx,lx,x);
+	ix = hx&0x7fffffff;
+	EXTRACT_WORDS(hy,ly,y);
+	iy = hy&0x7fffffff;
+	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
+	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
+	   return x+y;
+	if((hx-0x3ff00000|lx)==0) return atan(y);   /* x=1.0 */
+	m = ((hy>>31)&1)|((hx>>30)&2);	/* 2*sign(x)+sign(y) */
+
+    /* when y = 0 */
+	if((iy|ly)==0) {
+	    switch(m) {
+		case 0: 
+		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+	    }
+	}
+    /* when x = 0 */
+	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
+	    
+    /* when x is INF */
+	if(ix==0x7ff00000) {
+	    if(iy==0x7ff00000) {
+		switch(m) {
+		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
+		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
+		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
+	    }
+	}
+    /* when y is INF */
+	if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* compute y/x */
+	k = (iy-ix)>>20;
+	if(k > 60) z=pi_o_2+0.5*pi_lo; 	/* |y/x| >  2**60 */
+	else if(hx<0&&k<-60) z=0.0; 	/* |y|/x < -2**60 */
+	else z=atan(fabs(y/x));		/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: {
+	    	      __uint32_t zh;
+		      GET_HIGH_WORD(zh,z);
+		      SET_HIGH_WORD(z,zh ^ 0x80000000);
+		    }
+		    return       z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
@@ -0,0 +1,90 @@
+
+/* @(#)e_cosh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_cosh(x)
+ * Method : 
+ * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
+ *	2. 
+ *		                                        [ exp(x) - 1 ]^2 
+ *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
+ *			       			           2*exp(x)
+ *
+ *		                                  exp(x) +  1/exp(x)
+ *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
+ *			       			          2
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
+ *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
+ *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ *	cosh(x) is |x| if x is +INF, -INF, or NaN.
+ *	only cosh(0)=1 is exact for finite x.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const double one = 1.0, half=0.5, huge = 1.0e300;
+#else
+static double one = 1.0, half=0.5, huge = 1.0e300;
+#endif
+
+#ifdef __STDC__
+	double cosh(double x)
+#else
+	double cosh(x)
+	double x;
+#endif
+{	
+	double t,w;
+	__int32_t ix;
+	__uint32_t lx;
+
+    /* High word of |x|. */
+	GET_HIGH_WORD(ix,x);
+	ix &= 0x7fffffff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x*x;	
+
+    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+	if(ix<0x3fd62e43) {
+	    t = expm1(fabs(x));
+	    w = one+t;
+	    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
+	    return one+(t*t)/(w+w);
+	}
+
+    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+	if (ix < 0x40360000) {
+		t = exp(fabs(x));
+		return half*t+half/t;
+	}
+
+    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+	if (ix < 0x40862E42)  return half*exp(fabs(x));
+
+    /* |x| in [log(maxdouble), overflowthresold] */
+	GET_LOW_WORD(lx,x);
+	if (ix<0x408633CE || 
+             (ix==0x408633ce && lx<=(__uint32_t)0x8fb9f87d)) {
+	    w = exp(half*fabs(x));
+	    t = half*w;
+	    return t*w;
+	}
+
+    /* |x| > overflowthresold, cosh(x) overflow */
+	return huge*huge;
+}
+
@@ -0,0 +1,128 @@
+
+/* @(#)e_hypot.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_hypot(x,y)
+ *
+ * Method :                  
+ *	If (assume round-to-nearest) z=x*x+y*y 
+ *	has error less than sqrt(2)/2 ulp, than 
+ *	sqrt(z) has error less than 1 ulp (exercise).
+ *
+ *	So, compute sqrt(x*x+y*y) with some care as 
+ *	follows to get the error below 1 ulp:
+ *
+ *	Assume x>y>0;
+ *	(if possible, set rounding to round-to-nearest)
+ *	1. if x > 2y  use
+ *		x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
+ *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
+ *	2. if x <= 2y use
+ *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
+ *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
+ *	y1= y with lower 32 bits chopped, y2 = y-y1.
+ *		
+ *	NOTE: scaling may be necessary if some argument is too 
+ *	      large or too tiny
+ *
+ * Special cases:
+ *	hypot(x,y) is INF if x or y is +INF or -INF; else
+ *	hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypot(x,y) returns sqrt(x^2+y^2) with error less 
+ * 	than 1 ulps (units in the last place) 
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+	double __ieee754_hypot(double x, double y)
+#else
+	double __ieee754_hypot(x,y)
+	double x, y;
+#endif
+{
+	double a=x,b=y,t1,t2,y1,y2,w;
+	__int32_t j,k,ha,hb;
+
+	GET_HIGH_WORD(ha,x);
+	ha &= 0x7fffffff;
+	GET_HIGH_WORD(hb,y);
+	hb &= 0x7fffffff;
+	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
+	SET_HIGH_WORD(a,ha);	/* a <- |a| */
+	SET_HIGH_WORD(b,hb);	/* b <- |b| */
+	if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
+	k=0;
+	if(ha > 0x5f300000) {	/* a>2**500 */
+	   if(ha >= 0x7ff00000) {	/* Inf or NaN */
+	       __uint32_t low;
+	       w = a+b;			/* for sNaN */
+	       GET_LOW_WORD(low,a);
+	       if(((ha&0xfffff)|low)==0) w = a;
+	       GET_LOW_WORD(low,b);
+	       if(((hb^0x7ff00000)|low)==0) w = b;
+	       return w;
+	   }
+	   /* scale a and b by 2**-600 */
+	   ha -= 0x25800000; hb -= 0x25800000;	k += 600;
+	   SET_HIGH_WORD(a,ha);
+	   SET_HIGH_WORD(b,hb);
+	}
+	if(hb < 0x20b00000) {	/* b < 2**-500 */
+	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */	
+	        __uint32_t low;
+		GET_LOW_WORD(low,b);
+		if((hb|low)==0) return a;
+		t1=0;
+		SET_HIGH_WORD(t1,0x7fd00000);	/* t1=2^1022 */
+		b *= t1;
+		a *= t1;
+		k -= 1022;
+	    } else {		/* scale a and b by 2^600 */
+	        ha += 0x25800000; 	/* a *= 2^600 */
+		hb += 0x25800000;	/* b *= 2^600 */
+		k -= 600;
+		SET_HIGH_WORD(a,ha);
+		SET_HIGH_WORD(b,hb);
+	    }
+	}
+    /* medium size a and b */
+	w = a-b;
+	if (w>b) {
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha);
+	    t2 = a-t1;
+	    w  = __ieee754_sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
+	} else {
+	    a  = a+a;
+	    y1 = 0;
+	    SET_HIGH_WORD(y1,hb);
+	    y2 = b - y1;
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha+0x00100000);
+	    t2 = a - t1;
+	    w  = __ieee754_sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+	}
+	if(k!=0) {
+	    __uint32_t high;
+	    t1 = 1.0;
+	    GET_HIGH_WORD(high,t1);
+	    SET_HIGH_WORD(t1,high+(k<<20));
+	    return t1*w;
+	} else return w;
+}
+
+#endif /* defined(_DOUBLE_IS_32BITS) */
@@ -0,0 +1,83 @@
+
+/* @(#)e_sinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_sinh(x)
+ * Method : 
+ * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
+ *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
+ *	2. 
+ *		                                    E + E/(E+1)
+ *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
+ *			       			        2
+ *
+ *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2 
+ *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
+ *	    ln2ovft  <  x	    :  sinh(x) := x*shuge (overflow)
+ *
+ * Special cases:
+ *	sinh(x) is |x| if x is +INF, -INF, or NaN.
+ *	only sinh(0)=0 is exact for finite x.
+ */
+
+#include "fdlibm.h"
+
+#ifdef __STDC__
+static const double one = 1.0, shuge = 1.0e307;
+#else
+static double one = 1.0, shuge = 1.0e307;
+#endif
+
+#ifdef __STDC__
+	double sinh(double x)
+#else
+	double sinh(x)
+	double x;
+#endif
+{	
+	double t,w,h;
+	__int32_t ix,jx;
+	__uint32_t lx;
+
+    /* High word of |x|. */
+	GET_HIGH_WORD(jx,x);
+	ix = jx&0x7fffffff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x+x;	
+
+	h = 0.5;
+	if (jx<0) h = -h;
+    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
+	if (ix < 0x40360000) {		/* |x|<22 */
+	    if (ix<0x3e300000) 		/* |x|<2**-28 */
+		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    t = expm1(fabs(x));
+	    if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
+	    return h*(t+t/(t+one));
+	}
+
+    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
+	if (ix < 0x40862E42)  return h * exp(fabs(x));
+
+    /* |x| in [log(maxdouble), overflowthresold] */
+	GET_LOW_WORD(lx,x);
+       if (ix<0x408633CE || (ix==0x408633ce && lx<=(__uint32_t)0x8fb9f87d)) {
+	    w = exp(0.5*fabs(x));
+	    t = h*w;
+	    return t*w;
+	}
+
+    /* |x| > overflowthresold, sinh(x) overflow */
+	return x*shuge;
+}
+
@@ -0,0 +1,452 @@
+
+/* @(#)e_sqrt.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_sqrt(x)
+ * Return correctly rounded sqrt.
+ *           ------------------------------------------
+ *	     |  Use the hardware sqrt if you have one |
+ *           ------------------------------------------
+ * Method: 
+ *   Bit by bit method using integer arithmetic. (Slow, but portable) 
+ *   1. Normalization
+ *	Scale x to y in [1,4) with even powers of 2: 
+ *	find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+ *		sqrt(x) = 2^k * sqrt(y)
+ *   2. Bit by bit computation
+ *	Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+ *	     i							 0
+ *                                     i+1         2
+ *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
+ *	     i      i            i                 i
+ *                                                        
+ *	To compute q    from q , one checks whether 
+ *		    i+1       i                       
+ *
+ *			      -(i+1) 2
+ *			(q + 2      ) <= y.			(2)
+ *     			  i
+ *							      -(i+1)
+ *	If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+ *		 	       i+1   i             i+1   i
+ *
+ *	With some algebric manipulation, it is not difficult to see
+ *	that (2) is equivalent to 
+ *                             -(i+1)
+ *			s  +  2       <= y			(3)
+ *			 i                i
+ *
+ *	The advantage of (3) is that s  and y  can be computed by 
+ *				      i      i
+ *	the following recurrence formula:
+ *	    if (3) is false
+ *
+ *	    s     =  s  ,	y    = y   ;			(4)
+ *	     i+1      i		 i+1    i
+ *
+ *	    otherwise,
+ *                         -i                     -(i+1)
+ *	    s	  =  s  + 2  ,  y    = y  -  s  - 2  		(5)
+ *           i+1      i          i+1    i     i
+ *				
+ *	One may easily use induction to prove (4) and (5). 
+ *	Note. Since the left hand side of (3) contain only i+2 bits,
+ *	      it does not necessary to do a full (53-bit) comparison 
+ *	      in (3).
+ *   3. Final rounding
+ *	After generating the 53 bits result, we compute one more bit.
+ *	Together with the remainder, we can decide whether the
+ *	result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ *	(it will never equal to 1/2ulp).
+ *	The rounding mode can be detected by checking whether
+ *	huge + tiny is equal to huge, and whether huge - tiny is
+ *	equal to huge for some floating point number "huge" and "tiny".
+ *		
+ * Special cases:
+ *	sqrt(+-0) = +-0 	... exact
+ *	sqrt(inf) = inf
+ *	sqrt(-ve) = NaN		... with invalid signal
+ *	sqrt(NaN) = NaN		... with invalid signal for signaling NaN
+ *
+ * Other methods : see the appended file at the end of the program below.
+ *---------------
+ */
+
+#include "fdlibm.h"
+
+#ifndef _DOUBLE_IS_32BITS
+
+#ifdef __STDC__
+static	const double	one	= 1.0, tiny=1.0e-300;
+#else
+static	double	one	= 1.0, tiny=1.0e-300;
+#endif
+
+#ifdef __STDC__
+	double __ieee754_sqrt(double x)
+#else
+	double __ieee754_sqrt(x)
+	double x;
+#endif
+{
+	double z;
+	__int32_t sign = (int)0x80000000; 
+	__uint32_t r,t1,s1,ix1,q1;
+	__int32_t ix0,s0,q,m,t,i;
+
+	EXTRACT_WORDS(ix0,ix1,x);
+
+    /* take care of Inf and NaN */
+	if((ix0&0x7ff00000)==0x7ff00000) {			
+	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+					   sqrt(-inf)=sNaN */
+	} 
+    /* take care of zero */
+	if(ix0<=0) {
+	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
+	    else if(ix0<0)
+		return (x-x)/(x-x);		/* sqrt(-ve) = sNaN */
+	}
+    /* normalize x */
+	m = (ix0>>20);
+	if(m==0) {				/* subnormal x */
+	    while(ix0==0) {
+		m -= 21;
+		ix0 |= (ix1>>11); ix1 <<= 21;
+	    }
+	    for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
+	    m -= i-1;
+	    ix0 |= (ix1>>(32-i));
+	    ix1 <<= i;
+	}
+	m -= 1023;	/* unbias exponent */
+	ix0 = (ix0&0x000fffff)|0x00100000;
+	if(m&1){	/* odd m, double x to make it even */
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	}
+	m >>= 1;	/* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+	ix0 += ix0 + ((ix1&sign)>>31);
+	ix1 += ix1;
+	q = q1 = s0 = s1 = 0;	/* [q,q1] = sqrt(x) */
+	r = 0x00200000;		/* r = moving bit from right to left */
+
+	while(r!=0) {
+	    t = s0+r; 
+	    if(t<=ix0) { 
+		s0   = t+r; 
+		ix0 -= t; 
+		q   += r; 
+	    } 
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	    r>>=1;
+	}
+
+	r = sign;
+	while(r!=0) {
+	    t1 = s1+r; 
+	    t  = s0;
+	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
+		s1  = t1+r;
+		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
+		ix0 -= t;
+		if (ix1 < t1) ix0 -= 1;
+		ix1 -= t1;
+		q1  += r;
+	    }
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	    r>>=1;
+	}
+
+    /* use floating add to find out rounding direction */
+	if((ix0|ix1)!=0) {
+	    z = one-tiny; /* trigger inexact flag */
+	    if (z>=one) {
+	        z = one+tiny;
+	        if (q1==(__uint32_t)0xffffffff) { q1=0; q += 1;}
+		else if (z>one) {
+		    if (q1==(__uint32_t)0xfffffffe) q+=1;
+		    q1+=2; 
+		} else
+	            q1 += (q1&1);
+	    }
+	}
+	ix0 = (q>>1)+0x3fe00000;
+	ix1 =  q1>>1;
+	if ((q&1)==1) ix1 |= sign;
+	ix0 += (m <<20);
+	INSERT_WORDS(z,ix0,ix1);
+	return z;
+}
+ 
+#endif /* defined(_DOUBLE_IS_32BITS) */
+
+/*
+Other methods  (use floating-point arithmetic)
+-------------
+(This is a copy of a drafted paper by Prof W. Kahan 
+and K.C. Ng, written in May, 1986)
+
+	Two algorithms are given here to implement sqrt(x) 
+	(IEEE double precision arithmetic) in software.
+	Both supply sqrt(x) correctly rounded. The first algorithm (in
+	Section A) uses newton iterations and involves four divisions.
+	The second one uses reciproot iterations to avoid division, but
+	requires more multiplications. Both algorithms need the ability
+	to chop results of arithmetic operations instead of round them, 
+	and the INEXACT flag to indicate when an arithmetic operation
+	is executed exactly with no roundoff error, all part of the 
+	standard (IEEE 754-1985). The ability to perform shift, add,
+	subtract and logical AND operations upon 32-bit words is needed
+	too, though not part of the standard.
+
+A.  sqrt(x) by Newton Iteration
+
+   (1)	Initial approximation
+
+	Let x0 and x1 be the leading and the trailing 32-bit words of
+	a floating point number x (in IEEE double format) respectively 
+
+	    1    11		     52				  ...widths
+	   ------------------------------------------------------
+	x: |s|	  e     |	      f				|
+	   ------------------------------------------------------
+	      msb    lsb  msb				      lsb ...order
+
+ 
+	     ------------------------  	     ------------------------
+	x0:  |s|   e    |    f1     |	 x1: |          f2           |
+	     ------------------------  	     ------------------------
+
+	By performing shifts and subtracts on x0 and x1 (both regarded
+	as integers), we obtain an 8-bit approximation of sqrt(x) as
+	follows.
+
+		k  := (x0>>1) + 0x1ff80000;
+		y0 := k - T1[31&(k>>15)].	... y ~ sqrt(x) to 8 bits
+	Here k is a 32-bit integer and T1[] is an integer array containing
+	correction terms. Now magically the floating value of y (y's
+	leading 32-bit word is y0, the value of its trailing word is 0)
+	approximates sqrt(x) to almost 8-bit.
+
+	Value of T1:
+	static int T1[32]= {
+	0,	1024,	3062,	5746,	9193,	13348,	18162,	23592,
+	29598,	36145,	43202,	50740,	58733,	67158,	75992,	85215,
+	83599,	71378,	60428,	50647,	41945,	34246,	27478,	21581,
+	16499,	12183,	8588,	5674,	3403,	1742,	661,	130,};
+
+    (2)	Iterative refinement
+
+	Apply Heron's rule three times to y, we have y approximates 
+	sqrt(x) to within 1 ulp (Unit in the Last Place):
+
+		y := (y+x/y)/2		... almost 17 sig. bits
+		y := (y+x/y)/2		... almost 35 sig. bits
+		y := y-(y-x/y)/2	... within 1 ulp
+
+
+	Remark 1.
+	    Another way to improve y to within 1 ulp is:
+
+		y := (y+x/y)		... almost 17 sig. bits to 2*sqrt(x)
+		y := y - 0x00100006	... almost 18 sig. bits to sqrt(x)
+
+				2
+			    (x-y )*y
+		y := y + 2* ----------	...within 1 ulp
+			       2
+			     3y  + x
+
+
+	This formula has one division fewer than the one above; however,
+	it requires more multiplications and additions. Also x must be
+	scaled in advance to avoid spurious overflow in evaluating the
+	expression 3y*y+x. Hence it is not recommended uless division
+	is slow. If division is very slow, then one should use the 
+	reciproot algorithm given in section B.
+
+    (3) Final adjustment
+
+	By twiddling y's last bit it is possible to force y to be 
+	correctly rounded according to the prevailing rounding mode
+	as follows. Let r and i be copies of the rounding mode and
+	inexact flag before entering the square root program. Also we
+	use the expression y+-ulp for the next representable floating
+	numbers (up and down) of y. Note that y+-ulp = either fixed
+	point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+	mode.
+
+		I := FALSE;	... reset INEXACT flag I
+		R := RZ;	... set rounding mode to round-toward-zero
+		z := x/y;	... chopped quotient, possibly inexact
+		If(not I) then {	... if the quotient is exact
+		    if(z=y) {
+		        I := i;	 ... restore inexact flag
+		        R := r;  ... restore rounded mode
+		        return sqrt(x):=y.
+		    } else {
+			z := z - ulp;	... special rounding
+		    }
+		}
+		i := TRUE;		... sqrt(x) is inexact
+		If (r=RN) then z=z+ulp	... rounded-to-nearest
+		If (r=RP) then {	... round-toward-+inf
+		    y = y+ulp; z=z+ulp;
+		}
+		y := y+z;		... chopped sum
+		y0:=y0-0x00100000;	... y := y/2 is correctly rounded.
+	        I := i;	 		... restore inexact flag
+	        R := r;  		... restore rounded mode
+	        return sqrt(x):=y.
+		    
+    (4)	Special cases
+
+	Square root of +inf, +-0, or NaN is itself;
+	Square root of a negative number is NaN with invalid signal.
+
+
+B.  sqrt(x) by Reciproot Iteration
+
+   (1)	Initial approximation
+
+	Let x0 and x1 be the leading and the trailing 32-bit words of
+	a floating point number x (in IEEE double format) respectively
+	(see section A). By performing shifs and subtracts on x0 and y0,
+	we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
+
+	    k := 0x5fe80000 - (x0>>1);
+	    y0:= k - T2[63&(k>>14)].	... y ~ 1/sqrt(x) to 7.8 bits
+
+	Here k is a 32-bit integer and T2[] is an integer array 
+	containing correction terms. Now magically the floating
+	value of y (y's leading 32-bit word is y0, the value of
+	its trailing word y1 is set to zero) approximates 1/sqrt(x)
+	to almost 7.8-bit.
+
+	Value of T2:
+	static int T2[64]= {
+	0x1500,	0x2ef8,	0x4d67,	0x6b02,	0x87be,	0xa395,	0xbe7a,	0xd866,
+	0xf14a,	0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
+	0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
+	0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
+	0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
+	0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
+	0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
+	0x1527f,0x1334a,0x11051,0xe951,	0xbe01,	0x8e0d,	0x5924,	0x1edd,};
+
+    (2)	Iterative refinement
+
+	Apply Reciproot iteration three times to y and multiply the
+	result by x to get an approximation z that matches sqrt(x)
+	to about 1 ulp. To be exact, we will have 
+		-1ulp < sqrt(x)-z<1.0625ulp.
+	
+	... set rounding mode to Round-to-nearest
+	   y := y*(1.5-0.5*x*y*y)	... almost 15 sig. bits to 1/sqrt(x)
+	   y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
+	... special arrangement for better accuracy
+	   z := x*y			... 29 bits to sqrt(x), with z*y<1
+	   z := z + 0.5*z*(1-z*y)	... about 1 ulp to sqrt(x)
+
+	Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
+	(a) the term z*y in the final iteration is always less than 1; 
+	(b) the error in the final result is biased upward so that
+		-1 ulp < sqrt(x) - z < 1.0625 ulp
+	    instead of |sqrt(x)-z|<1.03125ulp.
+
+    (3)	Final adjustment
+
+	By twiddling y's last bit it is possible to force y to be 
+	correctly rounded according to the prevailing rounding mode
+	as follows. Let r and i be copies of the rounding mode and
+	inexact flag before entering the square root program. Also we
+	use the expression y+-ulp for the next representable floating
+	numbers (up and down) of y. Note that y+-ulp = either fixed
+	point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+	mode.
+
+	R := RZ;		... set rounding mode to round-toward-zero
+	switch(r) {
+	    case RN:		... round-to-nearest
+	       if(x<= z*(z-ulp)...chopped) z = z - ulp; else
+	       if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
+	       break;
+	    case RZ:case RM:	... round-to-zero or round-to--inf
+	       R:=RP;		... reset rounding mod to round-to-+inf
+	       if(x<z*z ... rounded up) z = z - ulp; else
+	       if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
+	       break;
+	    case RP:		... round-to-+inf
+	       if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
+	       if(x>z*z ...chopped) z = z+ulp;
+	       break;
+	}
+
+	Remark 3. The above comparisons can be done in fixed point. For
+	example, to compare x and w=z*z chopped, it suffices to compare
+	x1 and w1 (the trailing parts of x and w), regarding them as
+	two's complement integers.
+
+	...Is z an exact square root?
+	To determine whether z is an exact square root of x, let z1 be the
+	trailing part of z, and also let x0 and x1 be the leading and
+	trailing parts of x.
+
+	If ((z1&0x03ffffff)!=0)	... not exact if trailing 26 bits of z!=0
+	    I := 1;		... Raise Inexact flag: z is not exact
+	else {
+	    j := 1 - [(x0>>20)&1]	... j = logb(x) mod 2
+	    k := z1 >> 26;		... get z's 25-th and 26-th 
+					    fraction bits
+	    I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
+	}
+	R:= r		... restore rounded mode
+	return sqrt(x):=z.
+
+	If multiplication is cheaper then the foregoing red tape, the 
+	Inexact flag can be evaluated by
+
+	    I := i;
+	    I := (z*z!=x) or I.
+
+	Note that z*z can overwrite I; this value must be sensed if it is 
+	True.
+
+	Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
+	zero.
+
+		    --------------------
+		z1: |        f2        | 
+		    --------------------
+		bit 31		   bit 0
+
+	Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
+	or even of logb(x) have the following relations:
+
+	-------------------------------------------------
+	bit 27,26 of z1		bit 1,0 of x1	logb(x)
+	-------------------------------------------------
+	00			00		odd and even
+	01			01		even
+	10			10		odd
+	10			00		even
+	11			01		even
+	-------------------------------------------------
+
+    (4)	Special cases (see (4) of Section A).	
+ 
+ */
@@ -0,0 +1,299 @@
+/*							erfl.c
+ *
+ *	Error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, erfl();
+ *
+ * y = erfl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * The integral is
+ *
+ *                           x 
+ *                            -
+ *                 2         | |          2
+ *   erf(x)  =  --------     |    exp( - t  ) dt.
+ *              sqrt(pi)   | |
+ *                          -
+ *                           0
+ *
+ * The magnitude of x is limited to about 106.56 for IEEE
+ * arithmetic; 1 or -1 is returned outside this range.
+ *
+ * For 0 <= |x| < 1, erf(x) = x * P6(x^2)/Q6(x^2);
+ * Otherwise: erf(x) = 1 - erfc(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,1         50000       2.0e-19     5.7e-20
+ *
+ */
+
+/*							erfcl.c
+ *
+ *	Complementary error function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, erfcl();
+ *
+ * y = erfcl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *
+ *  1 - erf(x) =
+ *
+ *                           inf. 
+ *                             -
+ *                  2         | |          2
+ *   erfc(x)  =  --------     |    exp( - t  ) dt
+ *               sqrt(pi)   | |
+ *                           -
+ *                            x
+ *
+ *
+ * For small x, erfc(x) = 1 - erf(x); otherwise rational
+ * approximations are computed.
+ *
+ * A special function expx2l.c is used to suppress error amplification
+ * in computing exp(-x^2).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,13        50000      8.4e-19      9.7e-20
+ *    IEEE      6,106.56    20000      2.9e-19      7.1e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message          condition              value returned
+ * erfcl underflow    x^2 > MAXLOGL              0.0
+ *
+ *
+ */
+
+
+/*
+Modified from file ndtrl.c
+Cephes Math Library Release 2.3:  January, 1995
+Copyright 1984, 1995 by Stephen L. Moshier
+*/
+
+#include <math.h>
+#include "cephes_mconf.h"
+
+/* erfc(x) = exp(-x^2) P(1/x)/Q(1/x)
+   1/8 <= 1/x <= 1
+   Peak relative error 5.8e-21  */
+
+static const unsigned short P[] = {
+0x4bf0,0x9ad8,0x7a03,0x86c7,0x401d, XPD
+0xdf23,0xd843,0x4032,0x8881,0x401e, XPD
+0xd025,0xcfd5,0x8494,0x88d3,0x401e, XPD
+0xb6d0,0xc92b,0x5417,0xacb1,0x401d, XPD
+0xada8,0x356a,0x4982,0x94a6,0x401c, XPD
+0x4e13,0xcaee,0x9e31,0xb258,0x401a, XPD
+0x5840,0x554d,0x37a3,0x9239,0x4018, XPD
+0x3b58,0x3da2,0xaf02,0x9780,0x4015, XPD
+0x0144,0x489e,0xbe68,0x9c31,0x4011, XPD
+0x333b,0xd9e6,0xd404,0x986f,0xbfee, XPD
+};
+static const unsigned short Q[] = {
+/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */
+0x0e43,0x302d,0x79ed,0x86c7,0x401d, XPD
+0xf817,0x9128,0xc0f8,0xd48b,0x401e, XPD
+0x8eae,0x8dad,0x6eb4,0x9aa2,0x401f, XPD
+0x00e7,0x7595,0xcd06,0x88bb,0x401f, XPD
+0x4991,0xcfda,0x52f1,0xa2a9,0x401e, XPD
+0xc39d,0xe415,0xc43d,0x87c0,0x401d, XPD
+0xa75d,0x436f,0x30dd,0xa027,0x401b, XPD
+0xc4cb,0x305a,0xbf78,0x8220,0x4019, XPD
+0x3708,0x33b1,0x07fa,0x8644,0x4016, XPD
+0x24fa,0x96f6,0x7153,0x8a6c,0x4012, XPD
+};
+
+/* erfc(x) = exp(-x^2) 1/x R(1/x^2) / S(1/x^2)
+   1/128 <= 1/x < 1/8
+   Peak relative error 1.9e-21  */
+
+static const unsigned short R[] = {
+0x260a,0xab95,0x2fc7,0xe7c4,0x4000, XPD
+0x4761,0x613e,0xdf6d,0xe58e,0x4001, XPD
+0x0615,0x4b00,0x575f,0xdc7b,0x4000, XPD
+0x521d,0x8527,0x3435,0x8dc2,0x3ffe, XPD
+0x22cf,0xc711,0x6c5b,0xdcfb,0x3ff9, XPD
+};
+static const unsigned short S[] = {
+/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */
+0x5de6,0x17d7,0x54d6,0xaba9,0x4002, XPD
+0x55d5,0xd300,0xe71e,0xf564,0x4002, XPD
+0xb611,0x8f76,0xf020,0xd255,0x4001, XPD
+0x3684,0x3798,0xb793,0x80b0,0x3fff, XPD
+0xf5af,0x2fb2,0x1e57,0xc3d7,0x3ffa, XPD
+};
+
+/* erf(x)  = x T(x^2)/U(x^2)
+   0 <= x <= 1
+   Peak relative error 7.6e-23  */
+
+static const unsigned short T[] = {
+0xfd7a,0x3a1a,0x705b,0xe0c4,0x3ffb, XPD
+0x3128,0xc337,0x3716,0xace5,0x4001, XPD
+0x9517,0x4e93,0x540e,0x8f97,0x4007, XPD
+0x6118,0x6059,0x9093,0xa757,0x400a, XPD
+0xb954,0xa987,0xc60c,0xbc83,0x400e, XPD
+0x7a56,0xe45a,0xa4bd,0x975b,0x4010, XPD
+0xc446,0x6bab,0x0b2a,0x86d0,0x4013, XPD
+};
+
+static const unsigned short U[] = {
+/* 0x0000,0x0000,0x0000,0x8000,0x3fff, XPD */
+0x3453,0x1f8e,0xf688,0xb507,0x4004, XPD
+0x71ac,0xb12f,0x21ca,0xf2e2,0x4008, XPD
+0xffe8,0x9cac,0x3b84,0xc2ac,0x400c, XPD
+0x481d,0x445b,0xc807,0xc232,0x400f, XPD
+0x9ad5,0x1aef,0x45b1,0xe25e,0x4011, XPD
+0x71a7,0x1cad,0x012e,0xeef3,0x4012, XPD
+};
+
+/*							expx2l.c
+ *
+ *	Exponential of squared argument
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, expmx2l();
+ * int sign;
+ *
+ * y = expx2l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes y = exp(x*x) while suppressing error amplification
+ * that would ordinarily arise from the inexactness of the
+ * exponential argument x*x.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic      domain        # trials      peak         rms
+ *   IEEE     -106.566, 106.566    10^5       1.6e-19     4.4e-20
+ *
+ */
+
+#define M 32768.0L
+#define MINV 3.0517578125e-5L
+
+static long double expx2l (long double x)
+{
+  long double u, u1, m, f;
+
+  x = fabsl (x);
+  /* Represent x as an exact multiple of M plus a residual.
+     M is a power of 2 chosen so that exp(m * m) does not overflow
+     or underflow and so that |x - m| is small.  */
+  m = MINV * floorl(M * x + 0.5L);
+  f = x - m;
+
+  /* x^2 = m^2 + 2mf + f^2 */
+  u = m * m;
+  u1 = 2 * m * f  +  f * f;
+
+  if ((u+u1) > MAXLOGL)
+    return (INFINITYL);
+
+  /* u is exact, u1 is small.  */
+  u = expl(u) * expl(u1);
+  return(u);
+}
+
+long double erfcl(long double a)
+{
+long double p,q,x,y,z;
+
+if (isinf (a))
+	return (signbit (a) ? 2.0 : 0.0);
+
+x = fabsl (a);
+
+if (x < 1.0L)
+	return (1.0L - erfl(a));
+
+z = a * a;
+
+if( z  > MAXLOGL )
+	{
+under:
+	mtherr( "erfcl", UNDERFLOW );
+        errno = ERANGE;
+	return (signbit (a) ? 2.0 : 0.0);
+	}
+
+/* Compute z = expl(a * a).  */
+z = expx2l (a);
+y = 1.0L/x;
+
+if (x < 8.0L)
+	{
+	p = polevll (y, P, 9);
+	q = p1evll (y, Q, 10);
+	}
+else
+	{
+	q = y * y;
+	p = y * polevll (q, R, 4);
+	q = p1evll (q, S, 5);
+	}
+y = p/(q * z);
+
+if (a < 0.0L)
+	y = 2.0L - y;
+
+if (y == 0.0L)
+	goto under;
+
+return (y);
+}
+
+long double erfl(long double x)
+{
+long double y, z;
+
+if( x == 0.0L )
+	return (x);
+
+if (isinf (x))
+	return (signbit (x) ?  -1.0L : 1.0L);
+
+if (fabsl(x) > 1.0L)
+	return (1.0L - erfcl (x));
+
+z = x * x;
+y = x * polevll( z, T, 6 ) / p1evll( z, U, 6 );
+return( y );
+}
@@ -0,0 +1,49 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+/* e^x = 2^(x * log2(e)) */
+
+  .file "exp.s"
+	.text
+	.p2align 4,,15
+
+.globl _exp
+.def   _exp;  .scl  2;  .type 32; .endef
+
+_exp:
+
+      fldl  4(%esp)
+/* I added the following ugly construct because exp(+-Inf) resulted
+   in NaN.  The ugliness results from the bright minds at Intel.
+   For the i686 the code can be written better.
+   -- drepper@cygnus.com.  */
+    fxam        /* Is NaN or +-Inf?  */
+    fstsw %ax
+    movb  $0x45, %dh
+    andb  %ah, %dh
+    cmpb  $0x05, %dh
+    je  1f      /* Is +-Inf, jump.  */
+
+    fldl2e
+    fmulp       /* x * log2(e) */
+    fld %st
+    frndint       /* int(x * log2(e)) */
+    fsubr %st,%st(1)    /* fract(x * log2(e)) */
+    fxch
+    f2xm1       /* 2^(fract(x * log2(e))) - 1 */
+    fld1
+    faddp       /* 2^(fract(x * log2(e))) */
+    fscale        /* e^x */
+    fstp  %st(1)
+    ret
+
+1:
+    testl $0x200, %eax    /* Test sign.  */
+    jz  2f      /* If positive, jump.  */
+
+    fstp  %st
+    fldz        /* Set result to 0.  */
+2:
+    ret
@@ -0,0 +1,39 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ */
+
+	.file	"exp2.S"
+	.text
+	.align 4
+.globl _exp2
+	.def	_exp2;	.scl	2;	.type	32;	.endef
+_exp2:
+	fldl	4(%esp)
+/* I added the following ugly construct because exp(+-Inf) resulted
+   in NaN.  The ugliness results from the bright minds at Intel.
+   For the i686 the code can be written better.
+   -- drepper@cygnus.com.  */
+	fxam				/* Is NaN or +-Inf?  */
+	fstsw	%ax
+	movb	$0x45, %dh
+	andb	%ah, %dh
+	cmpb	$0x05, %dh
+	je	1f			/* Is +-Inf, jump.  */
+	fld	%st
+	frndint				/* int(x) */
+	fsubr	%st,%st(1)		/* fract(x) */
+	fxch
+	f2xm1				/* 2^(fract(x)) - 1 */
+	fld1
+	faddp				/* 2^(fract(x)) */
+	fscale				/* e^x */
+	fstp	%st(1)
+	ret
+
+1:	testl	$0x200, %eax		/* Test sign.  */
+	jz	2f			/* If positive, jump.  */
+	fstp	%st
+	fldz				/* Set result to 0.  */
+2:	ret
@@ -0,0 +1,39 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ */
+
+	.file	"exp2f.S"
+	.text
+	.align 4
+.globl _exp2f
+	.def	_exp2f;	.scl	2;	.type	32;	.endef
+_exp2f:
+	flds	4(%esp)
+/* I added the following ugly construct because exp(+-Inf) resulted
+   in NaN.  The ugliness results from the bright minds at Intel.
+   For the i686 the code can be written better.
+   -- drepper@cygnus.com.  */
+	fxam				/* Is NaN or +-Inf?  */
+	fstsw	%ax
+	movb	$0x45, %dh
+	andb	%ah, %dh
+	cmpb	$0x05, %dh
+	je	1f			/* Is +-Inf, jump.  */
+	fld	%st
+	frndint				/* int(x) */
+	fsubr	%st,%st(1)		/* fract(x) */
+	fxch
+	f2xm1				/* 2^(fract(x)) - 1 */
+	fld1
+	faddp				/* 2^(fract(x)) */
+	fscale				/* e^x */
+	fstp	%st(1)
+	ret
+
+1:	testl	$0x200, %eax		/* Test sign.  */
+	jz	2f			/* If positive, jump.  */
+	fstp	%st
+	fldz				/* Set result to 0.  */
+2:	ret
@@ -0,0 +1,39 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>.
+ * Public domain.
+ */
+
+	.file	"exp2l.S"
+	.text
+	.align 4
+.globl _exp2l
+	.def	_exp2l;	.scl	2;	.type	32;	.endef
+_exp2l:
+	fldt	4(%esp)
+/* I added the following ugly construct because exp(+-Inf) resulted
+   in NaN.  The ugliness results from the bright minds at Intel.
+   For the i686 the code can be written better.
+   -- drepper@cygnus.com.  */
+	fxam				/* Is NaN or +-Inf?  */
+	fstsw	%ax
+	movb	$0x45, %dh
+	andb	%ah, %dh
+	cmpb	$0x05, %dh
+	je	1f			/* Is +-Inf, jump.  */
+	fld	%st
+	frndint				/* int(x) */
+	fsubr	%st,%st(1)		/* fract(x) */
+	fxch
+	f2xm1				/* 2^(fract(x)) - 1 */
+	fld1
+	faddp				/* 2^(fract(x)) */
+	fscale				/* e^x */
+	fstp	%st(1)
+	ret
+
+1:	testl	$0x200, %eax		/* Test sign.  */
+	jz	2f			/* If positive, jump.  */
+	fstp	%st
+	fldz				/* Set result to 0.  */
+2:	ret
@@ -0,0 +1,3 @@
+#include <math.h>
+float expf (float x)
+  {return (float) exp (x);}
@@ -0,0 +1,71 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+/*
+ * The 8087 method for the exponential function is to calculate
+ *   exp(x) = 2^(x log2(e))
+ * after separating integer and fractional parts
+ *   x log2(e) = i + f, |f| <= .5
+ * 2^i is immediate but f needs to be precise for long double accuracy.
+ * Suppress range reduction error in computing f by the following.
+ * Separate x into integer and fractional parts
+ *   x = xi + xf, |xf| <= .5
+ * Separate log2(e) into the sum of an exact number c0 and small part c1.
+ *   c0 + c1 = log2(e) to extra precision
+ * Then
+ *   f = (c0 xi - i) + c0 xf + c1 x
+ * where c0 xi is exact and so also is (c0 xi - i).
+ * -- moshier@na-net.ornl.gov
+ */
+
+#include <math.h>
+#include "cephes_mconf.h"  /* for max and min log thresholds */
+
+static long double c0 = 1.44268798828125L;
+static long double c1 = 7.05260771340735992468e-6L;
+
+static long double
+__expl (long double x)
+{
+  long double res;
+  asm ("fldl2e\n\t"             /* 1  log2(e)         */
+       "fmul %%st(1),%%st\n\t"  /* 1  x log2(e)       */
+       "frndint\n\t"            /* 1  i               */
+       "fld %%st(1)\n\t"        /* 2  x               */
+       "frndint\n\t"            /* 2  xi              */
+       "fld %%st(1)\n\t"        /* 3  i               */
+       "fldt %2\n\t"            /* 4  c0              */
+       "fld %%st(2)\n\t"        /* 5  xi              */
+       "fmul %%st(1),%%st\n\t"  /* 5  c0 xi           */
+       "fsubp %%st,%%st(2)\n\t" /* 4  f = c0 xi  - i  */
+       "fld %%st(4)\n\t"        /* 5  x               */
+       "fsub %%st(3),%%st\n\t"  /* 5  xf = x - xi     */
+       "fmulp %%st,%%st(1)\n\t" /* 4  c0 xf           */
+       "faddp %%st,%%st(1)\n\t" /* 3  f = f + c0 xf   */
+       "fldt %3\n\t"            /* 4                  */
+       "fmul %%st(4),%%st\n\t"  /* 4  c1 * x          */
+       "faddp %%st,%%st(1)\n\t" /* 3  f = f + c1 * x  */
+       "f2xm1\n\t"		/* 3 2^(fract(x * log2(e))) - 1 */
+       "fld1\n\t"               /* 4 1.0              */
+       "faddp\n\t"		/* 3 2^(fract(x * log2(e))) */
+       "fstp	%%st(1)\n\t"    /* 2  */
+       "fscale\n\t"	        /* 2 scale factor is st(1); e^x */
+       "fstp	%%st(1)\n\t"    /* 1  */
+       "fstp	%%st(1)\n\t"    /* 0  */
+       : "=t" (res) : "0" (x), "m" (c0), "m" (c1) : "ax", "dx");
+  return res;
+}
+
+long double expl (long double x)
+{
+  if (x > MAXLOGL)
+    return INFINITY;
+  else if (x < MINLOGL)
+    return 0.0L;
+  else
+    return __expl (x);
+}
@@ -0,0 +1,28 @@
+/*
+ * Written 2005 by Gregory W. Chicares <chicares@cox.net>.
+ * Adapted to double by Danny Smith <dannysmith@users.sourceforge.net>. 
+ * Public domain.
+ *
+ * F2XM1's input is constrained to (-1, +1), so the domain of
+ * 'x * LOG2EL' is (-LOGE2L, +LOGE2L). Outside that domain,
+ * delegating to exp() handles C99 7.12.6.3/2 range errors.
+ *
+ * Constants from moshier.net, file cephes/ldouble/constl.c,
+ * are used instead of M_LN2 and M_LOG2E, which would not be
+ * visible with 'gcc std=c99'.  The use of these extended precision
+ * constants also allows gcc to replace them with x87 opcodes.
+ */
+
+#include <math.h> /* expl() */
+#include "cephes_mconf.h"
+double expm1 (double x)
+{
+  if (fabs(x) < LOGE2L)
+    {
+      x *= LOG2EL;
+      __asm__("f2xm1" : "=t" (x) : "0" (x));
+      return x;
+    }
+  else
+    return exp(x) - 1.0;
+}
@@ -0,0 +1,29 @@
+/*
+ * Written 2005 by Gregory W. Chicares <chicares@cox.net>.
+ * Adapted to float by Danny Smith <dannysmith@users.sourceforge.net>. 
+ * Public domain.
+ *
+ * F2XM1's input is constrained to (-1, +1), so the domain of
+ * 'x * LOG2EL' is (-LOGE2L, +LOGE2L). Outside that domain,
+ * delegating to exp() handles C99 7.12.6.3/2 range errors.
+ *
+ * Constants from moshier.net, file cephes/ldouble/constl.c,
+ * are used instead of M_LN2 and M_LOG2E, which would not be
+ * visible with 'gcc std=c99'. The use of these extended precision
+ * constants also allows gcc to replace them with x87 opcodes.
+ */
+
+#include <math.h> /* expl() */
+#include "cephes_mconf.h"
+
+float expm1f (float x)
+{
+  if (fabsf(x) < LOGE2L)
+    {
+      x *= LOG2EL;
+      __asm__("f2xm1" : "=t" (x) : "0" (x));
+      return x;
+    }
+  else
+    return expf(x) - 1.0F;
+}
@@ -0,0 +1,29 @@
+/*
+ * Written 2005 by Gregory W. Chicares <chicares@cox.net> with
+ * help from Danny Smith. dannysmith@users.sourceforge.net>.
+ * Public domain.
+ *
+ * F2XM1's input is constrained to (-1, +1), so the domain of
+ * 'x * LOG2EL' is (-LOGE2L, +LOGE2L). Outside that domain,
+ * delegating to expl() handles C99 7.12.6.3/2 range errors.
+ *
+ * Constants from moshier.net, file cephes/ldouble/constl.c,
+ * are used instead of M_LN2 and M_LOG2E, which would not be
+ * visible with 'gcc std=c99'.  The use of these extended precision
+ * constants also allows gcc to replace them with x87 opcodes.
+ */
+
+#include <math.h> /* expl() */
+#include "cephes_mconf.h"
+
+long double expm1l (long double x)
+{
+  if (fabsl(x) < LOGE2L)
+    {
+      x *= LOG2EL;
+      __asm__("f2xm1" : "=t" (x) : "0" (x));
+      return x;
+    }
+  else
+    return expl(x) - 1.0L;
+}
@@ -0,0 +1,10 @@
+#include <math.h>
+
+double
+fabs (double x)
+{
+  double res;
+
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
@@ -0,0 +1,9 @@
+#include <math.h>
+
+float
+fabsf (float x)
+{
+  float res;
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
@@ -0,0 +1,9 @@
+#include <math.h>
+
+long double
+fabsl (long double x)
+{
+  long double res;
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
@@ -0,0 +1,115 @@
+#ifndef _MINGWEX_FASTMATH_H_
+#define _MINGWEX_FASTMATH_H_
+
+/* Fast math inlines
+   No range or domain checks. No setting of errno.  No tweaks to
+   protect precision near range limits. */
+
+/* For now this is an internal header with just the functions that
+   are currently used in building libmingwex.a math components */
+
+/* FIXME: We really should get rid of the code duplication using euther
+   C++ templates or tgmath-type macros.  */  
+
+static __inline__ double __fast_sqrt (double x)
+{
+  double res;
+  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
+
+static __inline__ long double __fast_sqrtl (long double x)
+{
+  long double res;
+  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
+
+static __inline__ float __fast_sqrtf (float x)
+{
+  float res;
+  asm __volatile__ ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
+
+
+static __inline__ double __fast_log (double x)
+{
+   double res;
+   asm __volatile__
+     ("fldln2\n\t"
+      "fxch\n\t"
+      "fyl2x"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+static __inline__ long double __fast_logl (long double x)
+{
+  long double res;
+   asm __volatile__
+     ("fldln2\n\t"
+      "fxch\n\t"
+      "fyl2x"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+
+static __inline__ float __fast_logf (float x)
+{
+   float res;
+   asm __volatile__
+     ("fldln2\n\t"
+      "fxch\n\t"
+      "fyl2x"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+static __inline__ double __fast_log1p (double x)
+{
+  double res;
+  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
+  if (fabs (x) >= 1.0 - 0.5 * 1.41421356237309504880)
+    res = __fast_log (1.0 + x);
+  else
+    asm __volatile__
+      ("fldln2\n\t"
+       "fxch\n\t"
+       "fyl2xp1"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+static __inline__ long double __fast_log1pl (long double x)
+{
+  long double res;
+  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
+  if (fabsl (x) >= 1.0L - 0.5L * 1.41421356237309504880L)
+    res = __fast_logl (1.0L + x);
+  else
+    asm __volatile__
+      ("fldln2\n\t"
+       "fxch\n\t"
+       "fyl2xp1"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+static __inline__ float __fast_log1pf (float x)
+{
+  float res;
+  /* fyl2xp1 accurate only for |x| <= 1.0 - 0.5 * sqrt (2.0) */
+  if (fabsf (x) >= 1.0 - 0.5 * 1.41421356237309504880)
+    res = __fast_logf (1.0 + x);
+  else
+    asm __volatile__
+      ("fldln2\n\t"
+       "fxch\n\t"
+       "fyl2xp1"
+       : "=t" (res) : "0" (x) : "st(1)");
+   return res;
+}
+
+#endif
@@ -0,0 +1,7 @@
+#include <math.h>
+
+double
+fdim (double x, double y)
+{
+  return  (isgreater(x, y) ? (x - y) : 0.0);
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+float
+fdimf (float x, float y)
+{
+  return  (isgreater(x, y) ? (x - y) : 0.0F);
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+long double
+fdiml (long double x, long double y)
+{
+   return  (isgreater(x, y) ? (x - y) : 0.0L);
+}
@@ -0,0 +1,365 @@
+
+/* @(#)fdlibm.h 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+/* REDHAT LOCAL: Include files.  */
+#include <math.h>
+#include <sys/types.h>
+#include <machine/ieeefp.h>
+
+/* REDHAT LOCAL: Default to XOPEN_MODE.  */
+#define _XOPEN_MODE
+
+/* Most routines need to check whether a float is finite, infinite, or not a
+   number, and many need to know whether the result of an operation will
+   overflow.  These conditions depend on whether the largest exponent is
+   used for NaNs & infinities, or whether it's used for finite numbers.  The
+   macros below wrap up that kind of information:
+
+   FLT_UWORD_IS_FINITE(X)
+	True if a positive float with bitmask X is finite.
+
+   FLT_UWORD_IS_NAN(X)
+	True if a positive float with bitmask X is not a number.
+
+   FLT_UWORD_IS_INFINITE(X)
+	True if a positive float with bitmask X is +infinity.
+
+   FLT_UWORD_MAX
+	The bitmask of FLT_MAX.
+
+   FLT_UWORD_HALF_MAX
+	The bitmask of FLT_MAX/2.
+
+   FLT_UWORD_EXP_MAX
+	The bitmask of the largest finite exponent (129 if the largest
+	exponent is used for finite numbers, 128 otherwise).
+
+   FLT_UWORD_LOG_MAX
+	The bitmask of log(FLT_MAX), rounded down.  This value is the largest
+	input that can be passed to exp() without producing overflow.
+
+   FLT_UWORD_LOG_2MAX
+	The bitmask of log(2*FLT_MAX), rounded down.  This value is the
+	largest input than can be passed to cosh() without producing
+	overflow.
+
+   FLT_LARGEST_EXP
+	The largest biased exponent that can be used for finite numbers
+	(255 if the largest exponent is used for finite numbers, 254
+	otherwise) */
+
+#ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
+#define FLT_UWORD_IS_FINITE(x) 1
+#define FLT_UWORD_IS_NAN(x) 0
+#define FLT_UWORD_IS_INFINITE(x) 0
+#define FLT_UWORD_MAX 0x7fffffff
+#define FLT_UWORD_EXP_MAX 0x43010000
+#define FLT_UWORD_LOG_MAX 0x42b2d4fc
+#define FLT_UWORD_LOG_2MAX 0x42b437e0
+#define HUGE ((float)0X1.FFFFFEP128)
+#else
+#define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
+#define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
+#define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
+#define FLT_UWORD_MAX 0x7f7fffffL
+#define FLT_UWORD_EXP_MAX 0x43000000
+#define FLT_UWORD_LOG_MAX 0x42b17217
+#define FLT_UWORD_LOG_2MAX 0x42b2d4fc
+#define HUGE ((float)3.40282346638528860e+38)
+#endif
+#define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))
+#define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
+
+/* Many routines check for zero and subnormal numbers.  Such things depend
+   on whether the target supports denormals or not:
+
+   FLT_UWORD_IS_ZERO(X)
+	True if a positive float with bitmask X is +0.	Without denormals,
+	any float with a zero exponent is a +0 representation.	With
+	denormals, the only +0 representation is a 0 bitmask.
+
+   FLT_UWORD_IS_SUBNORMAL(X)
+	True if a non-zero positive float with bitmask X is subnormal.
+	(Routines should check for zeros first.)
+
+   FLT_UWORD_MIN
+	The bitmask of the smallest float above +0.  Call this number
+	REAL_FLT_MIN...
+
+   FLT_UWORD_EXP_MIN
+	The bitmask of the float representation of REAL_FLT_MIN's exponent.
+
+   FLT_UWORD_LOG_MIN
+	The bitmask of |log(REAL_FLT_MIN)|, rounding down.
+
+   FLT_SMALLEST_EXP
+	REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
+	-22 if they are).
+*/
+
+#ifdef _FLT_NO_DENORMALS
+#define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
+#define FLT_UWORD_IS_SUBNORMAL(x) 0
+#define FLT_UWORD_MIN 0x00800000
+#define FLT_UWORD_EXP_MIN 0x42fc0000
+#define FLT_UWORD_LOG_MIN 0x42aeac50
+#define FLT_SMALLEST_EXP 1
+#else
+#define FLT_UWORD_IS_ZERO(x) ((x)==0)
+#define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
+#define FLT_UWORD_MIN 0x00000001
+#define FLT_UWORD_EXP_MIN 0x43160000
+#define FLT_UWORD_LOG_MIN 0x42cff1b5
+#define FLT_SMALLEST_EXP -22
+#endif
+
+#ifdef __STDC__
+#undef __P
+#define	__P(p)	p
+#else
+#define	__P(p)	()
+#endif
+
+/* 
+ * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
+ * (one may replace the following line by "#include <values.h>")
+ */
+
+#define X_TLOSS		1.41484755040568800000e+16 
+
+/* Functions that are not documented, and are not in <math.h>.  */
+
+extern double logb __P((double));
+#ifdef _SCALB_INT
+extern double scalb __P((double, int));
+#else
+extern double scalb __P((double, double));
+#endif
+extern double significand __P((double));
+
+/* ieee style elementary functions */
+extern double __ieee754_sqrt __P((double));			
+extern double __ieee754_acos __P((double));			
+extern double __ieee754_acosh __P((double));			
+extern double __ieee754_log __P((double));			
+extern double __ieee754_atanh __P((double));			
+extern double __ieee754_asin __P((double));			
+extern double __ieee754_atan2 __P((double,double));			
+extern double __ieee754_exp __P((double));
+extern double __ieee754_cosh __P((double));
+extern double __ieee754_fmod __P((double,double));
+extern double __ieee754_pow __P((double,double));
+extern double __ieee754_lgamma_r __P((double,int *));
+extern double __ieee754_gamma_r __P((double,int *));
+extern double __ieee754_log10 __P((double));
+extern double __ieee754_sinh __P((double));
+extern double __ieee754_hypot __P((double,double));
+extern double __ieee754_j0 __P((double));
+extern double __ieee754_j1 __P((double));
+extern double __ieee754_y0 __P((double));
+extern double __ieee754_y1 __P((double));
+extern double __ieee754_jn __P((int,double));
+extern double __ieee754_yn __P((int,double));
+extern double __ieee754_remainder __P((double,double));
+extern __int32_t __ieee754_rem_pio2 __P((double,double*));
+#ifdef _SCALB_INT
+extern double __ieee754_scalb __P((double,int));
+#else
+extern double __ieee754_scalb __P((double,double));
+#endif
+
+/* fdlibm kernel function */
+extern double __kernel_standard __P((double,double,int));
+extern double __kernel_sin __P((double,double,int));
+extern double __kernel_cos __P((double,double));
+extern double __kernel_tan __P((double,double,int));
+extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
+
+/* Undocumented float functions.  */
+extern float logbf __P((float));
+#ifdef _SCALB_INT
+extern float scalbf __P((float, int));
+#else
+extern float scalbf __P((float, float));
+#endif
+extern float significandf __P((float));
+
+/* ieee style elementary float functions */
+extern float __ieee754_sqrtf __P((float));			
+extern float __ieee754_acosf __P((float));			
+extern float __ieee754_acoshf __P((float));			
+extern float __ieee754_logf __P((float));			
+extern float __ieee754_atanhf __P((float));			
+extern float __ieee754_asinf __P((float));			
+extern float __ieee754_atan2f __P((float,float));			
+extern float __ieee754_expf __P((float));
+extern float __ieee754_coshf __P((float));
+extern float __ieee754_fmodf __P((float,float));
+extern float __ieee754_powf __P((float,float));
+extern float __ieee754_lgammaf_r __P((float,int *));
+extern float __ieee754_gammaf_r __P((float,int *));
+extern float __ieee754_log10f __P((float));
+extern float __ieee754_sinhf __P((float));
+extern float __ieee754_hypotf __P((float,float));
+extern float __ieee754_j0f __P((float));
+extern float __ieee754_j1f __P((float));
+extern float __ieee754_y0f __P((float));
+extern float __ieee754_y1f __P((float));
+extern float __ieee754_jnf __P((int,float));
+extern float __ieee754_ynf __P((int,float));
+extern float __ieee754_remainderf __P((float,float));
+extern __int32_t __ieee754_rem_pio2f __P((float,float*));
+#ifdef _SCALB_INT
+extern float __ieee754_scalbf __P((float,int));
+#else
+extern float __ieee754_scalbf __P((float,float));
+#endif
+
+/* float versions of fdlibm kernel functions */
+extern float __kernel_sinf __P((float,float,int));
+extern float __kernel_cosf __P((float,float));
+extern float __kernel_tanf __P((float,float,int));
+extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
+
+/* The original code used statements like
+	n0 = ((*(int*)&one)>>29)^1;		* index of high word *
+	ix0 = *(n0+(int*)&x);			* high word of x *
+	ix1 = *((1-n0)+(int*)&x);		* low word of x *
+   to dig two 32 bit words out of the 64 bit IEEE floating point
+   value.  That is non-ANSI, and, moreover, the gcc instruction
+   scheduler gets it wrong.  We instead use the following macros.
+   Unlike the original code, we determine the endianness at compile
+   time, not at run time; I don't see much benefit to selecting
+   endianness at run time.  */
+
+#ifndef __IEEE_BIG_ENDIAN
+#ifndef __IEEE_LITTLE_ENDIAN
+ #error Must define endianness
+#endif
+#endif
+
+/* A union which permits us to convert between a double and two 32 bit
+   ints.  */
+
+#ifdef __IEEE_BIG_ENDIAN
+
+typedef union 
+{
+  double value;
+  struct 
+  {
+    __uint32_t msw;
+    __uint32_t lsw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+#ifdef __IEEE_LITTLE_ENDIAN
+
+typedef union 
+{
+  double value;
+  struct 
+  {
+    __uint32_t lsw;
+    __uint32_t msw;
+  } parts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)				\
+do {								\
+  ieee_double_shape_type ew_u;					\
+  ew_u.value = (d);						\
+  (ix0) = ew_u.parts.msw;					\
+  (ix1) = ew_u.parts.lsw;					\
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gh_u;					\
+  gh_u.value = (d);						\
+  (i) = gh_u.parts.msw;						\
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gl_u;					\
+  gl_u.value = (d);						\
+  (i) = gl_u.parts.lsw;						\
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d,ix0,ix1)					\
+do {								\
+  ieee_double_shape_type iw_u;					\
+  iw_u.parts.msw = (ix0);					\
+  iw_u.parts.lsw = (ix1);					\
+  (d) = iw_u.value;						\
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  sh_u.parts.msw = (v);						\
+  (d) = sh_u.value;						\
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sl_u;					\
+  sl_u.value = (d);						\
+  sl_u.parts.lsw = (v);						\
+  (d) = sl_u.value;						\
+} while (0)
+
+/* A union which permits us to convert between a float and a 32 bit
+   int.  */
+
+typedef union
+{
+  float value;
+  __uint32_t word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)					\
+do {								\
+  ieee_float_shape_type gf_u;					\
+  gf_u.value = (d);						\
+  (i) = gf_u.word;						\
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)					\
+do {								\
+  ieee_float_shape_type sf_u;					\
+  sf_u.word = (i);						\
+  (d) = sf_u.value;						\
+} while (0)
@@ -0,0 +1,33 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ * 
+ */
+	.file	"floor.s"
+	.text
+	.align 4
+.globl _floor
+	.def	_floor;	.scl	2;	.type	32;	.endef
+_floor:
+	fldl	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x400,%edx		/* round towards -oo */
+	orl	4(%esp),%edx
+	andl	$0xf7ff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,35 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ * 
+ * Removed header file dependency for use in libmingwex.a by
+ *   Danny Smith <dannysmith@users.sourceforge.net>
+ */
+	.file	"floorf.S"
+	.text
+	.align 4
+.globl _floorf
+	.def	_floorf;	.scl	2;	.type	32;	.endef
+_floorf:
+	flds	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x400,%edx		/* round towards -oo */
+	orl	4(%esp),%edx
+	andl	$0xf7ff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,33 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ * 
+ */
+	.file	"floorl.S"
+	.text
+	.align 4
+.globl _floorl
+	.def	_floorl;	.scl	2;	.type	32;	.endef
+_floorl:
+	fldt	4(%esp)
+	subl	$8,%esp
+
+	fstcw	4(%esp)			/* store fpu control word */
+
+	/* We use here %edx although only the low 1 bits are defined.
+	   But none of the operations should care and they are faster
+	   than the 16 bit operations.  */
+	movl	$0x400,%edx		/* round towards -oo */
+	orl	4(%esp),%edx
+	andl	$0xf7ff,%edx
+	movl	%edx,(%esp)
+	fldcw	(%esp)			/* load modified control word */
+
+	frndint				/* round */
+
+	fldcw	4(%esp)			/* restore original control word */
+
+	addl	$8,%esp
+	ret
@@ -0,0 +1,12 @@
+	.file	"fma.S"
+	.text
+	.align 2
+	.p2align 4,,15
+.globl _fma
+	.def	_fma;	.scl	2;	.type	32;	.endef
+_fma:
+	fldl	4(%esp)
+	fmull	12(%esp)
+	fldl	20(%esp)
+	faddp
+	ret
@@ -0,0 +1,12 @@
+	.file	"fmaf.S"
+	.text
+	.align 2
+	.p2align 4,,15
+.globl _fmaf
+	.def	_fmaf;	.scl	2;	.type	32;	.endef
+_fmaf:
+	flds	4(%esp)
+	fmuls	8(%esp)
+	flds	12(%esp)
+	faddp
+	ret
@@ -0,0 +1,5 @@
+long double
+fmal ( long double _x,  long double _y,  long double _z)
+{
+  return ((_x * _y) + _z);
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+double
+fmax (double _x, double _y)
+{
+   return ( isgreaterequal (_x, _y)|| __isnan (_y) ?  _x : _y );
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+float
+fmaxf (float _x, float _y)
+{
+  return (( isgreaterequal(_x, _y) || __isnanf (_y)) ?  _x : _y );
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+long double
+fmaxl (long double _x, long double  _y)
+{
+  return (( isgreaterequal(_x, _y) || __isnanl (_y)) ?  _x : _y );
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+double
+fmin (double _x, double _y)
+{
+  return ((islessequal(_x, _y) || __isnan (_y)) ? _x : _y );
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+float
+fminf (float _x, float _y)
+{
+  return ((islessequal(_x, _y) || isnan (_y)) ? _x : _y );
+}
@@ -0,0 +1,7 @@
+#include <math.h>
+
+long double
+fminl (long double _x, long double _y)
+{
+  return ((islessequal(_x, _y) || __isnanl (_y)) ? _x : _y );
+}
@@ -0,0 +1,37 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for float type by Danny Smith
+ *  <dannysmith@users.sourceforge.net>.
+ */
+
+#include <math.h>
+
+float
+fmodf (float x, float y)
+{
+  float res;
+
+  asm ("1:\tfprem\n\t"
+       "fstsw   %%ax\n\t"
+       "sahf\n\t"
+       "jp      1b\n\t"
+       "fstp    %%st(1)"
+       : "=t" (res) : "0" (x), "u" (y) : "ax", "st(1)");
+  return res;
+}
+
+double
+fmod (double x, double y)
+{
+  float res;
+
+  asm ("1:\tfprem\n\t"
+       "fstsw   %%ax\n\t"
+       "sahf\n\t"
+       "jp      1b\n\t"
+       "fstp    %%st(1)"
+       : "=t" (res) : "0" (x), "u" (y) : "ax", "st(1)");
+  return res;
+}
@@ -0,0 +1,22 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
+ */
+
+#include <math.h>
+
+long double
+fmodl (long double x, long double y)
+{
+  long double res;
+
+  asm ("1:\tfprem\n\t"
+       "fstsw   %%ax\n\t"
+       "sahf\n\t"
+       "jp      1b\n\t"
+       "fstp    %%st(1)"
+       : "=t" (res) : "0" (x), "u" (y) : "ax", "st(1)");
+  return res;
+}
@@ -0,0 +1,14 @@
+
+#include "fp_consts.h"
+const union _ieee_rep __QNAN = { __DOUBLE_QNAN_REP };
+const union _ieee_rep __SNAN = { __DOUBLE_SNAN_REP };
+const union _ieee_rep __INF = { __DOUBLE_INF_REP };
+const union _ieee_rep __DENORM = { __DOUBLE_DENORM_REP };
+
+/* ISO C99 */
+#undef nan
+/* FIXME */
+double nan (const char * tagp __attribute__((unused)) )
+	{ return __QNAN.double_val; } 
+
+
@@ -0,0 +1,48 @@
+#ifndef _FP_CONSTS_H
+#define _FP_CONSTS_H
+
+/*
+According to IEEE 754 a QNaN has exponent bits of all 1 values and
+initial significand bit of 1.  A SNaN has has an exponent of all 1
+values and initial significand bit of 0 (with one or more other
+significand bits of 1). An Inf has significand of 0 and
+exponent of all 1 values. A denormal value has all exponent bits of 0.
+
+The following does _not_ follow those rules, but uses values
+equal to those exported from MS C++ runtime lib, msvcprt.dll
+for float and double. MSVC however, does not have long doubles.
+*/
+
+
+#define __DOUBLE_INF_REP { 0, 0, 0, 0x7ff0 }
+#define __DOUBLE_QNAN_REP { 0, 0, 0, 0xfff8 }  /* { 0, 0, 0, 0x7ff8 }  */
+#define __DOUBLE_SNAN_REP { 0, 0, 0, 0xfff0 }  /* { 1, 0, 0, 0x7ff0 }  */
+#define __DOUBLE_DENORM_REP {1, 0, 0, 0}
+
+#define D_NAN_MASK 0x7ff0000000000000LL /* this will mask NaN's and Inf's */
+
+#define __FLOAT_INF_REP { 0, 0x7f80 }
+#define __FLOAT_QNAN_REP { 0, 0xffc0 }  /* { 0, 0x7fc0 }  */
+#define __FLOAT_SNAN_REP { 0, 0xff80 }  /* { 1, 0x7f80 }  */
+#define __FLOAT_DENORM_REP {1,0}
+
+#define F_NAN_MASK 0x7f800000
+
+/*
+   This assumes no implicit (hidden) bit in extended mode.
+   Padded to 96 bits
+ */
+#define __LONG_DOUBLE_INF_REP { 0, 0, 0, 0x8000, 0x7fff, 0 }
+#define __LONG_DOUBLE_QNAN_REP { 0, 0, 0, 0xc000, 0xffff, 0 } 
+#define __LONG_DOUBLE_SNAN_REP { 0, 0, 0, 0x8000, 0xffff, 0 }
+#define __LONG_DOUBLE_DENORM_REP {1, 0, 0, 0, 0, 0}
+
+union _ieee_rep
+{
+  unsigned short rep[6];
+  float float_val;
+  double double_val;
+  long double ldouble_val;
+} ;
+
+#endif
@@ -0,0 +1,12 @@
+#include "fp_consts.h"
+
+const union _ieee_rep __QNANF = { __FLOAT_QNAN_REP };   
+const union _ieee_rep __SNANF = { __FLOAT_SNAN_REP };   
+const union _ieee_rep __INFF = { __FLOAT_INF_REP };   
+const union _ieee_rep __DENORMF = { __FLOAT_DENORM_REP };   
+
+/* ISO C99 */
+#undef nanf
+/* FIXME */
+float nanf(const char * tagp __attribute__((unused)) )
+  { return __QNANF.float_val;}
@@ -0,0 +1,12 @@
+#include "fp_consts.h"
+
+const union _ieee_rep __QNANL = { __LONG_DOUBLE_QNAN_REP };
+const union _ieee_rep __SNANL  = { __LONG_DOUBLE_SNAN_REP };
+const union _ieee_rep __INFL = { __LONG_DOUBLE_INF_REP };
+const union _ieee_rep __DENORML = { __LONG_DOUBLE_DENORM_REP };
+
+
+#undef nanl
+/* FIXME */
+long double nanl (const char * tagp __attribute__((unused)) )
+  { return __QNANL.ldouble_val; } 
@@ -0,0 +1,20 @@
+#include <math.h>
+
+/* 'fxam' sets FPU flags C3,C2,C0   'fstsw' stores: 
+ FP_NAN			001		0x0100
+ FP_NORMAL		010		0x0400
+ FP_INFINITE		011		0x0500
+ FP_ZERO		100		0x4000
+ FP_SUBNORMAL		110		0x4400
+
+and sets C1 flag (signbit) if neg */
+
+int __fpclassify (double _x){
+  unsigned short sw;
+  __asm__ (
+	"fxam; fstsw %%ax;"
+	: "=a" (sw)
+	: "t" (_x)
+	);
+  return sw & (FP_NAN | FP_NORMAL | FP_ZERO );
+}
@@ -0,0 +1,10 @@
+#include <math.h>
+int __fpclassifyf (float _x){
+  unsigned short sw;
+  __asm__ (
+	"fxam; fstsw %%ax;"
+	: "=a" (sw)
+	: "t" (_x)
+	);
+  return sw & (FP_NAN | FP_NORMAL | FP_ZERO );
+}
@@ -0,0 +1,10 @@
+#include <math.h>
+int __fpclassifyl (long double _x){
+  unsigned short sw;
+  __asm__ (
+	"fxam; fstsw %%ax;"
+	: "=a" (sw)
+	: "t" (_x)
+	);
+  return sw & (FP_NAN | FP_NORMAL | FP_ZERO );
+}
@@ -0,0 +1,74 @@
+/* ix87 specific frexp implementation for double.
+   Copyright (C) 1997, 2000, 2005 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, write to the Free
+   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+   02111-1307 USA.  */
+
+
+  .file "frexp.s"
+	.text
+
+.align 4
+two54:	.byte 0, 0, 0, 0, 0, 0, 0x50, 0x43
+
+	.p2align 4,,15
+
+.globl _frexp
+.def   _frexp;  .scl  2;  .type 32; .endef
+
+_frexp:
+
+      movl  4(%esp), %ecx
+      movl  8(%esp), %eax
+      movl  %eax, %edx
+      andl  $0x7fffffff, %eax
+      orl %eax, %ecx
+      jz  1f
+
+      xorl  %ecx, %ecx
+      cmpl  $0x7ff00000, %eax
+      jae 1f
+
+      cmpl  $0x00100000, %eax
+      jae 2f
+
+      fld  4(%esp)
+
+      fmull two54
+      movl  $-54, %ecx
+      fstpl 4(%esp)
+      fwait
+      movl  8(%esp), %eax
+      movl  %eax, %edx
+      andl  $0x7fffffff, %eax
+
+2:
+      shrl  $20, %eax
+      andl  $0x800fffff, %edx
+      subl  $1022, %eax
+      orl $0x3fe00000, %edx
+      addl  %eax, %ecx
+      movl  %edx, 8(%esp)
+
+	/* Store %ecx in the variable pointed to by the second argument,
+	   get the factor from the stack and return.  */
+1:
+      movl  12(%esp), %eax
+      fldl  4(%esp)
+      movl  %ecx, (%eax)
+
+      ret
@@ -0,0 +1,3 @@
+#include <math.h>
+float frexpf (float x, int* expn)
+  {return (float)frexp(x, expn);}
@@ -0,0 +1,71 @@
+/*
+  Cephes Math Library Release 2.7:  May, 1998
+  Copyright 1984, 1987, 1988, 1992, 1998 by Stephen L. Moshier
+
+  Extracted from floorl.387 for use in libmingwex.a by
+  Danny Smith <dannysmith@users.sourceforge.net>
+  2002-06-20
+*/
+
+/*
+ * frexpl(long double x, int* expnt) extracts the exponent from x.
+ * It returns an integer power of two to expnt and the significand
+ * between 0.5 and 1 to y.  Thus  x = y * 2**expn.
+ */ 
+	.align 2
+.globl _frexpl
+_frexpl:
+	pushl %ebp
+	movl %esp,%ebp
+	subl $24,%esp
+	pushl %esi
+	pushl %ebx
+	fldt 8(%ebp)
+	movl 20(%ebp),%ebx
+	fld %st(0)
+	fstpt -12(%ebp)
+	leal -4(%ebp),%ecx
+	movw -4(%ebp),%dx
+	andl $32767,%edx
+	jne L25
+	fldz
+	fucompp
+	fnstsw %ax
+	andb $68,%ah
+	xorb $64,%ah
+	jne L21
+	movl $0,(%ebx)
+	fldz
+	jmp L24
+	.align 2,0x90
+	.align 2,0x90
+L21:
+	fldt -12(%ebp)
+	fadd %st(0),%st
+	fstpt -12(%ebp)
+	decl %edx
+	movw (%ecx),%si
+	andl $32767,%esi
+	jne L22
+	cmpl $-66,%edx
+	jg L21
+L22:
+	addl %esi,%edx
+	jmp L19
+	.align 2,0x90
+L25:
+	fstp %st(0)
+L19:
+	addl $-16382,%edx
+	movl %edx,(%ebx)
+	movw (%ecx),%ax
+	andl $-32768,%eax
+	orl $16382,%eax
+	movw %ax,(%ecx)
+	fldt -12(%ebp)
+L24:
+	leal -32(%ebp),%esp
+	popl %ebx
+	popl %esi
+	leave
+	ret
@@ -0,0 +1,11 @@
+int 
+__fp_unordered_compare (long double x,  long double y){
+  unsigned short retval;
+  __asm__ (
+	"fucom %%st(1);"
+	"fnstsw;"
+	: "=a" (retval)
+	: "t" (x), "u" (y)
+	);
+  return retval;
+}
@@ -0,0 +1,75 @@
+#include <math.h>
+#include "fdlibm.h"
+
+float hypotf (float x, float y)
+{
+	double a=x,b=y,t1,t2,y1,y2,w;
+	__int32_t j,k,ha,hb;
+
+	GET_HIGH_WORD(ha,x);
+	ha &= 0x7fffffff;
+	GET_HIGH_WORD(hb,y);
+	hb &= 0x7fffffff;
+	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
+	SET_HIGH_WORD(a,ha);	/* a <- |a| */
+	SET_HIGH_WORD(b,hb);	/* b <- |b| */
+	if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
+	k=0;
+	if(ha > 0x5f300000) {	/* a>2**500 */
+	   if(ha >= 0x7ff00000) {	/* Inf or NaN */
+	       __uint32_t low;
+	       w = a+b;			/* for sNaN */
+	       GET_LOW_WORD(low,a);
+	       if(((ha&0xfffff)|low)==0) w = a;
+	       GET_LOW_WORD(low,b);
+	       if(((hb^0x7ff00000)|low)==0) w = b;
+	       return w;
+	   }
+	   /* scale a and b by 2**-600 */
+	   ha -= 0x25800000; hb -= 0x25800000;	k += 600;
+	   SET_HIGH_WORD(a,ha);
+	   SET_HIGH_WORD(b,hb);
+	}
+	if(hb < 0x20b00000) {	/* b < 2**-500 */
+	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */	
+	        __uint32_t low;
+		GET_LOW_WORD(low,b);
+		if((hb|low)==0) return a;
+		t1=0;
+		SET_HIGH_WORD(t1,0x7fd00000);	/* t1=2^1022 */
+		b *= t1;
+		a *= t1;
+		k -= 1022;
+	    } else {		/* scale a and b by 2^600 */
+	        ha += 0x25800000; 	/* a *= 2^600 */
+		hb += 0x25800000;	/* b *= 2^600 */
+		k -= 600;
+		SET_HIGH_WORD(a,ha);
+		SET_HIGH_WORD(b,hb);
+	    }
+	}
+    /* medium size a and b */
+	w = a-b;
+	if (w>b) {
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha);
+	    t2 = a-t1;
+	    w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
+	} else {
+	    a  = a+a;
+	    y1 = 0;
+	    SET_HIGH_WORD(y1,hb);
+	    y2 = b - y1;
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha+0x00100000);
+	    t2 = a - t1;
+	    w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+	}
+	if(k!=0) {
+	    __uint32_t high;
+	    t1 = 1.0;
+	    GET_HIGH_WORD(high,t1);
+	    SET_HIGH_WORD(t1,high+(k<<20));
+	    return t1*w;
+	} else return w;
+}
@@ -0,0 +1,73 @@
+#include <math.h>
+#include <float.h>
+#include <errno.h>
+
+/*
+This implementation is based largely on Cephes library
+function cabsl (cmplxl.c), which bears the following notice:
+
+Cephes Math Library Release 2.1:  January, 1989
+Copyright 1984, 1987, 1989 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/*
+   Modified for use in libmingwex.a
+   02 Sept 2002  Danny Smith  <dannysmith@users.sourceforege.net>
+   Calls to ldexpl replaced by logbl and calls to frexpl replaced
+   by scalbnl to avoid duplicated range checks.
+*/
+
+extern long double __INFL;
+#define PRECL 32
+
+long double
+hypotl (long double x, long double y)
+{
+  int exx;
+  int eyy;
+  int  scale;
+  long double xx =fabsl(x);
+  long double yy =fabsl(y);
+  if (!isfinite(xx) || !isfinite(yy))
+    return  xx + yy; /* Return INF or NAN. */
+
+  if (xx == 0.0L)
+     return yy;
+  if (yy == 0.0L)
+     return xx;
+
+  /* Get exponents */
+  exx =  logbl (xx);
+  eyy =  logbl (yy);
+
+  /* Check if large differences in scale */
+  scale = exx - eyy;
+  if ( scale > PRECL)
+     return xx;
+  if ( scale < -PRECL)
+     return yy;
+
+  /* Exponent of approximate geometric mean (x 2) */
+  scale = (exx + eyy) >> 1;
+
+  /*  Rescale: Geometric mean is now about 2 */  
+  x = scalbnl(xx, -scale);
+  y = scalbnl(yy, -scale);
+
+  xx = sqrtl(x * x  + y * y);
+
+  /* Check for overflow and underflow */
+  exx = logbl(xx);   
+  exx += scale;
+    if (exx > LDBL_MAX_EXP)
+    {
+      errno = ERANGE; 
+      return __INFL;
+    }
+  if (exx < LDBL_MIN_EXP)
+    return 0.0L;
+
+  /* Undo scaling */
+  return (scalbnl (xx, scale));
+}
@@ -0,0 +1,37 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+
+	.file	"ilogb.S"
+	.text
+	.align 4
+.globl _ilogb
+	.def	_ilogb;	.scl	2;	.type	32;	.endef
+_ilogb:
+
+	fldl	4(%esp)
+/* I added the following ugly construct because ilogb(+-Inf) is
+   required to return INT_MAX in ISO C99.
+   -- jakub@redhat.com.  */
+	fxam			/* Is NaN or +-Inf?  */
+	fstsw   %ax
+	movb    $0x45, %dh
+	andb    %ah, %dh
+	cmpb    $0x05, %dh
+	je      1f		/* Is +-Inf, jump.  */
+
+	fxtract
+	pushl	%eax
+	fstp	%st
+
+	fistpl	(%esp)
+	fwait
+	popl	%eax
+
+	ret
+
+1:	fstp	%st
+	movl	$0x7fffffff, %eax
+	ret
@@ -0,0 +1,35 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ */
+
+	.file	"ilogbf.S"
+	.text
+	.align 4
+.globl _ilogbf
+	.def	_ilogbf;	.scl	2;	.type	32;	.endef
+_ilogbf:
+	flds	4(%esp)
+/* I added the following ugly construct because ilogb(+-Inf) is
+   required to return INT_MAX in ISO C99.
+   -- jakub@redhat.com.  */
+	fxam			/* Is NaN or +-Inf?  */
+	fstsw   %ax
+	movb    $0x45, %dh
+	andb    %ah, %dh
+	cmpb    $0x05, %dh
+	je      1f		/* Is +-Inf, jump.  */
+
+	fxtract
+	pushl	%eax
+	fstp	%st
+
+	fistpl	(%esp)
+	fwait
+	popl	%eax
+
+	ret
+
+1:	fstp	%st
+	movl	$0x7fffffff, %eax
+	ret
@@ -0,0 +1,36 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Changes for long double by Ulrich Drepper <drepper@cygnus.com>
+ * Public domain.
+ */
+
+	.file	"ilogbl.S"
+	.text
+	.align 4
+.globl _ilogbl
+	.def	_ilogbl;	.scl	2;	.type	32;	.endef
+_ilogbl:
+	fldt	4(%esp)
+/* I added the following ugly construct because ilogb(+-Inf) is
+   required to return INT_MAX in ISO C99.
+   -- jakub@redhat.com.  */
+	fxam			/* Is NaN or +-Inf?  */
+	fstsw   %ax
+	movb    $0x45, %dh
+	andb    %ah, %dh
+	cmpb    $0x05, %dh
+	je      1f		/* Is +-Inf, jump.  */
+
+	fxtract
+	pushl	%eax
+	fstp	%st
+
+	fistpl	(%esp)
+	fwait
+	popl	%eax
+
+	ret
+
+1:	fstp	%st
+	movl	$0x7fffffff, %eax
+	ret
@@ -0,0 +1,14 @@
+#include <math.h>
+
+int
+__isnan (double _x)
+{
+  unsigned short _sw;
+  __asm__ ("fxam;"
+	   "fstsw %%ax": "=a" (_sw) : "t" (_x));
+  return (_sw & (FP_NAN | FP_NORMAL | FP_INFINITE | FP_ZERO | FP_SUBNORMAL))
+    == FP_NAN;
+}
+
+#undef isnan
+int __attribute__ ((alias ("__isnan"))) isnan (double);
@@ -0,0 +1,12 @@
+#include <math.h>
+int
+__isnanf (float _x)
+{
+  unsigned short _sw;
+  __asm__ ("fxam;"
+	   "fstsw %%ax": "=a" (_sw) : "t" (_x) );
+  return (_sw & (FP_NAN | FP_NORMAL | FP_INFINITE | FP_ZERO | FP_SUBNORMAL))
+    == FP_NAN;
+}
+
+int __attribute__ ((alias ("__isnanf"))) isnanf (float);
@@ -0,0 +1,13 @@
+#include <math.h>
+
+int
+__isnanl (long double _x)
+{
+  unsigned short _sw;
+  __asm__ ("fxam;"
+	   "fstsw %%ax": "=a" (_sw) : "t" (_x));
+  return (_sw & (FP_NAN | FP_NORMAL | FP_INFINITE | FP_ZERO | FP_SUBNORMAL))
+    == FP_NAN;
+}
+
+int __attribute__ ((alias ("__isnanl"))) isnanl (long double);
@@ -0,0 +1,19 @@
+#include <math.h>
+#include <errno.h>
+
+double ldexp(double x, int expn)
+{
+  double res;
+  if (!isfinite (x) || x == 0.0L)
+    return x;
+
+  __asm__ ("fscale"
+  	    : "=t" (res)
+        : "0" (x), "u" ((double) expn));
+
+ // if (!isfinite (res) || res == 0.0L)
+ //   errno = ERANGE;
+
+  return res;
+}
+
@@ -0,0 +1,3 @@
+#include <math.h>
+float ldexpf (float x, int expn)
+  {return (float) ldexp (x, expn);}
@@ -0,0 +1,19 @@
+#include <math.h>
+#include <errno.h>
+  
+long double ldexpl(long double x, int expn)
+{
+  long double res;
+  if (!isfinite (x) || x == 0.0L)
+    return x;
+
+  __asm__ ("fscale"
+  	    : "=t" (res)
+	    : "0" (x), "u" ((long double) expn));
+
+  if (!isfinite (res) || res == 0.0L)
+    errno = ERANGE;
+
+  return res;
+}
+
@@ -0,0 +1,369 @@
+/*							lgam()
+ *
+ *	Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, __lgamma_r();
+ * int* sgngam;
+ * y = __lgamma_r( x, sgngam );
+ * 
+ * double x, y, lgamma();
+ * y = lgamma( 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 the variable referenced by sgngam.
+ *
+ * For arguments greater than 13, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula using a polynomial approximation of
+ * degree 4. Arguments between -33 and +33 are reduced by
+ * recurrence to the interval [2,3] of a rational approximation.
+ * The cosecant reflection formula is employed for arguments
+ * less than -33.
+ *
+ * Arguments greater than MAXLGM return MAXNUM and an error
+ * message.  MAXLGM = 2.035093e36 for DEC
+ * arithmetic or 2.556348e305 for IEEE arithmetic.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ * arithmetic      domain        # trials     peak         rms
+ *    DEC     0, 3                  7000     5.2e-17     1.3e-17
+ *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
+ *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
+ *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ *
+ * The following test used the relative error criterion, though
+ * at certain points the relative error could be much higher than
+ * indicated.
+ *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
+ *
+ */
+
+/*
+ * Cephes Math Library Release 2.8:  June, 2000
+ * Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+ */
+
+/*
+ * 26-11-2002 Modified for mingw.
+ * Danny Smith <dannysmith@users.sourceforge.net>
+ */
+
+
+#ifndef __MINGW32__
+#include "mconf.h"
+#ifdef ANSIPROT
+extern double pow ( double, double );
+extern double log ( double );
+extern double exp ( double );
+extern double sin ( double );
+extern double polevl ( double, void *, int );
+extern double p1evl ( double, void *, int );
+extern double floor ( double );
+extern double fabs ( double );
+extern int isnan ( double );
+extern int isfinite ( double );
+#else
+double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs();
+int isnan(), isfinite();
+#endif
+#ifdef INFINITIES
+extern double INFINITY;
+#endif
+#ifdef NANS
+extern double NAN;
+#endif
+#else /* __MINGW32__ */
+#include "cephes_mconf.h"
+#endif /* __MINGW32__ */
+
+
+/* A[]: Stirling's formula expansion of log gamma
+ * B[], C[]: log gamma function between 2 and 3
+ */
+#ifdef UNK
+static double A[] = {
+ 8.11614167470508450300E-4,
+-5.95061904284301438324E-4,
+ 7.93650340457716943945E-4,
+-2.77777777730099687205E-3,
+ 8.33333333333331927722E-2
+};
+static double B[] = {
+-1.37825152569120859100E3,
+-3.88016315134637840924E4,
+-3.31612992738871184744E5,
+-1.16237097492762307383E6,
+-1.72173700820839662146E6,
+-8.53555664245765465627E5
+};
+static double C[] = {
+/* 1.00000000000000000000E0, */
+-3.51815701436523470549E2,
+-1.70642106651881159223E4,
+-2.20528590553854454839E5,
+-1.13933444367982507207E6,
+-2.53252307177582951285E6,
+-2.01889141433532773231E6
+};
+/* log( sqrt( 2*pi ) ) */
+static double LS2PI  =  0.91893853320467274178;
+#define MAXLGM 2.556348e305
+static double LOGPI = 1.14472988584940017414;
+#endif
+
+#ifdef DEC
+static const unsigned short A[] = {
+0035524,0141201,0034633,0031405,
+0135433,0176755,0126007,0045030,
+0035520,0006371,0003342,0172730,
+0136066,0005540,0132605,0026407,
+0037252,0125252,0125252,0125132
+};
+static const unsigned short B[] = {
+0142654,0044014,0077633,0035410,
+0144027,0110641,0125335,0144760,
+0144641,0165637,0142204,0047447,
+0145215,0162027,0146246,0155211,
+0145322,0026110,0010317,0110130,
+0145120,0061472,0120300,0025363
+};
+static const unsigned short C[] = {
+/*0040200,0000000,0000000,0000000*/
+0142257,0164150,0163630,0112622,
+0143605,0050153,0156116,0135272,
+0144527,0056045,0145642,0062332,
+0145213,0012063,0106250,0001025,
+0145432,0111254,0044577,0115142,
+0145366,0071133,0050217,0005122
+};
+/* log( sqrt( 2*pi ) ) */
+static const unsigned short LS2P[] = {040153,037616,041445,0172645,};
+#define LS2PI *(double *)LS2P
+#define MAXLGM 2.035093e36
+static const unsigned short LPI[4] = {
+0040222,0103202,0043475,0006750,
+};
+#define LOGPI *(double *)LPI
+
+#endif
+
+#ifdef IBMPC
+static const unsigned short A[] = {
+0x6661,0x2733,0x9850,0x3f4a,
+0xe943,0xb580,0x7fbd,0xbf43,
+0x5ebb,0x20dc,0x019f,0x3f4a,
+0xa5a1,0x16b0,0xc16c,0xbf66,
+0x554b,0x5555,0x5555,0x3fb5
+};
+static const unsigned short B[] = {
+0x6761,0x8ff3,0x8901,0xc095,
+0xb93e,0x355b,0xf234,0xc0e2,
+0x89e5,0xf890,0x3d73,0xc114,
+0xdb51,0xf994,0xbc82,0xc131,
+0xf20b,0x0219,0x4589,0xc13a,
+0x055e,0x5418,0x0c67,0xc12a
+};
+static const unsigned short C[] = {
+/*0x0000,0x0000,0x0000,0x3ff0,*/
+0x12b2,0x1cf3,0xfd0d,0xc075,
+0xd757,0x7b89,0xaa0d,0xc0d0,
+0x4c9b,0xb974,0xeb84,0xc10a,
+0x0043,0x7195,0x6286,0xc131,
+0xf34c,0x892f,0x5255,0xc143,
+0xe14a,0x6a11,0xce4b,0xc13e
+};
+/* log( sqrt( 2*pi ) ) */
+static const union
+{
+  unsigned short  s[4];
+  double d;
+} ls2p  =  {{0xbeb5,0xc864,0x67f1,0x3fed}};
+#define LS2PI   (ls2p.d)
+#define MAXLGM 2.556348e305
+/* log (pi) */
+static const union
+{
+  unsigned short s[4];
+  double d;
+} lpi  =  {{0xa1bd,0x48e7,0x50d0,0x3ff2}};
+#define LOGPI (lpi.d)
+#endif
+
+#ifdef MIEEE
+static const unsigned short A[] = {
+0x3f4a,0x9850,0x2733,0x6661,
+0xbf43,0x7fbd,0xb580,0xe943,
+0x3f4a,0x019f,0x20dc,0x5ebb,
+0xbf66,0xc16c,0x16b0,0xa5a1,
+0x3fb5,0x5555,0x5555,0x554b
+};
+static const unsigned short B[] = {
+0xc095,0x8901,0x8ff3,0x6761,
+0xc0e2,0xf234,0x355b,0xb93e,
+0xc114,0x3d73,0xf890,0x89e5,
+0xc131,0xbc82,0xf994,0xdb51,
+0xc13a,0x4589,0x0219,0xf20b,
+0xc12a,0x0c67,0x5418,0x055e
+};
+static const unsigned short C[] = {
+0xc075,0xfd0d,0x1cf3,0x12b2,
+0xc0d0,0xaa0d,0x7b89,0xd757,
+0xc10a,0xeb84,0xb974,0x4c9b,
+0xc131,0x6286,0x7195,0x0043,
+0xc143,0x5255,0x892f,0xf34c,
+0xc13e,0xce4b,0x6a11,0xe14a
+};
+/* log( sqrt( 2*pi ) ) */
+static const union
+{
+  unsigned short  s[4];
+  double d;
+} ls2p  =  {{0x3fed,0x67f1,0xc864,0xbeb5}};
+#define LS2PI  ls2p.d
+#define MAXLGM 2.556348e305
+/* log (pi) */
+static const union
+{
+  unsigned short s[4];
+  double d;
+} lpi  =  {{0x3ff2, 0x50d0, 0x48e7, 0xa1bd}};
+#define LOGPI (lpi.d)
+#endif
+
+
+/* Logarithm of gamma function */
+/* Reentrant version */ 
+
+double __lgamma_r(double x, int* sgngam)
+{
+double p, q, u, w, z;
+int i;
+
+*sgngam = 1;
+#ifdef NANS
+if( isnan(x) )
+	return(x);
+#endif
+
+#ifdef INFINITIES
+if( !isfinite(x) )
+	return(INFINITY);
+#endif
+
+if( x < -34.0 )
+	{
+	q = -x;
+	w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
+	p = floor(q);
+	if( p == q )
+		{
+lgsing:
+		_SET_ERRNO(EDOM);
+		mtherr( "lgam", SING );
+#ifdef INFINITIES
+		return (INFINITY);
+#else
+		return (MAXNUM);
+#endif
+		}
+	i = p;
+	if( (i & 1) == 0 )
+		*sgngam = -1;
+	else
+		*sgngam = 1;
+	z = q - p;
+	if( z > 0.5 )
+		{
+		p += 1.0;
+		z = p - q;
+		}
+	z = q * sin( PI * z );
+	if( z == 0.0 )
+		goto lgsing;
+/*	z = log(PI) - log( z ) - w;*/
+	z = LOGPI - log( z ) - w;
+	return( z );
+	}
+
+if( x < 13.0 )
+	{
+	z = 1.0;
+	p = 0.0;
+	u = x;
+	while( u >= 3.0 )
+		{
+		p -= 1.0;
+		u = x + p;
+		z *= u;
+		}
+	while( u < 2.0 )
+		{
+		if( u == 0.0 )
+			goto lgsing;
+		z /= u;
+		p += 1.0;
+		u = x + p;
+		}
+	if( z < 0.0 )
+		{
+		*sgngam = -1;
+		z = -z;
+		}
+	else
+		*sgngam = 1;
+	if( u == 2.0 )
+		return( log(z) );
+	p -= 2.0;
+	x = x + p;
+	p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
+	return( log(z) + p );
+	}
+
+if( x > MAXLGM )
+	{
+	_SET_ERRNO(ERANGE);
+	mtherr( "lgamma", OVERFLOW );
+#ifdef INFINITIES
+	return( *sgngam * INFINITY );
+#else
+	return( *sgngam * MAXNUM );
+#endif
+	}
+
+q = ( x - 0.5 ) * log(x) - x + LS2PI;
+if( x > 1.0e8 )
+	return( q );
+
+p = 1.0/(x*x);
+if( x >= 1000.0 )
+	q += ((   7.9365079365079365079365e-4 * p
+		- 2.7777777777777777777778e-3) *p
+		+ 0.0833333333333333333333) / x;
+else
+	q += polevl( p, A, 4 ) / x;
+return( q );
+}
+
+/* This is the C99 version */
+
+double lgamma(double x)
+{
+  int local_sgngam=0;
+  return (__lgamma_r(x, &local_sgngam));
+}
@@ -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));
+}
@@ -0,0 +1,416 @@
+/*							lgaml()
+ *
+ *	Natural logarithm of gamma function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, __lgammal_r();
+ * int* sgngaml;
+ * y = __lgammal_r( x, sgngaml );
+ * 
+ * long double x, y, lgammal();
+ * y = lgammal( 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 the variable referenced by sgngaml.
+ *
+ * For arguments greater than 33, the logarithm of the gamma
+ * function is approximated by the logarithmic version of
+ * Stirling's formula using a polynomial approximation of
+ * degree 4. Arguments between -33 and +33 are reduced by
+ * recurrence to the interval [2,3] of a rational approximation.
+ * The cosecant reflection formula is employed for arguments
+ * less than -33.
+ *
+ * Arguments greater than MAXLGML (10^4928) return MAXNUML.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ * arithmetic      domain        # trials     peak         rms
+ *    IEEE         -40, 40        100000     2.2e-19     4.6e-20
+ *    IEEE    10^-2000,10^+2000    20000     1.6e-19     3.3e-20
+ * The error criterion was relative when the function magnitude
+ * was greater than one but absolute when it was less than one.
+ *
+ */
+
+/*
+ * Copyright 1994 by Stephen L. Moshier
+ */
+
+/*
+ * 26-11-2002 Modified for mingw.
+ * Danny Smith <dannysmith@users.sourceforge.net>
+ */
+
+#ifndef __MINGW32__
+#include "mconf.h"
+#ifdef ANSIPROT
+extern long double fabsl ( long double );
+extern long double lgaml ( long double );
+extern long double logl ( long double );
+extern long double expl ( long double );
+extern long double gammal ( long double );
+extern long double sinl ( long double );
+extern long double floorl ( long double );
+extern long double powl ( long double, long double );
+extern long double polevll ( long double, void *, int );
+extern long double p1evll ( long double, void *, int );
+extern int isnanl ( long double );
+extern int isfinitel ( long double );
+#else
+long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl();
+long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel();
+#endif
+#ifdef INFINITIES
+extern long double INFINITYL;
+#endif
+#ifdef NANS
+extern long double NANL;
+#endif
+#else /* __MINGW32__ */
+#include "cephes_mconf.h"
+#endif /* __MINGW32__ */
+
+#if UNK
+static long double S[9] = {
+-1.193945051381510095614E-3L,
+ 7.220599478036909672331E-3L,
+-9.622023360406271645744E-3L,
+-4.219773360705915470089E-2L,
+ 1.665386113720805206758E-1L,
+-4.200263503403344054473E-2L,
+-6.558780715202540684668E-1L,
+ 5.772156649015328608253E-1L,
+ 1.000000000000000000000E0L,
+};
+#endif
+#if IBMPC
+static const unsigned short S[] = {
+0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD
+0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD
+0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD
+0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD
+0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD
+0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD
+0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD
+0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
+0x0000,0x0000,0x0000,0x8000,0x3fff, XPD
+};
+#endif
+#if MIEEE
+static long S[27] = {
+0xbff50000,0x9c7e25e5,0xd6d3baeb,
+0x3ff70000,0xec9ac74e,0xceb4fe9a,
+0xbff80000,0x9da5b0e9,0xdfef9225,
+0xbffa0000,0xacd787dc,0xec1710b0,
+0x3ffc0000,0xaa891905,0x75156b8d,
+0xbffa0000,0xac0af47d,0x126bf183,
+0xbffe0000,0xa7e7a013,0x57d17bf6,
+0x3ffe0000,0x93c467e3,0x7db0c7a9,
+0x3fff0000,0x80000000,0x00000000,
+};
+#endif
+
+#if UNK
+static long double SN[9] = {
+ 1.133374167243894382010E-3L,
+ 7.220837261893170325704E-3L,
+ 9.621911155035976733706E-3L,
+-4.219773343731191721664E-2L,
+-1.665386113944413519335E-1L,
+-4.200263503402112910504E-2L,
+ 6.558780715202536547116E-1L,
+ 5.772156649015328608727E-1L,
+-1.000000000000000000000E0L,
+};
+#endif
+#if IBMPC
+static const unsigned SN[] = {
+0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD
+0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD
+0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD
+0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD
+0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD
+0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD
+0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD
+0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD
+0x0000,0x0000,0x0000,0x8000,0xbfff, XPD
+};
+#endif
+#if MIEEE
+static long SN[27] = {
+0x3ff50000,0x948db9f7,0x02de5dd1,
+0x3ff70000,0xec9cc5f1,0xdd68989b,
+0x3ff80000,0x9da5386f,0x18f02ca1,
+0xbffa0000,0xacd787d1,0x41dd783f,
+0xbffc0000,0xaa891905,0xd76d7a5b,
+0xbffa0000,0xac0af47d,0x12347f64,
+0x3ffe0000,0xa7e7a013,0x57d15e26,
+0x3ffe0000,0x93c467e3,0x7db0c7aa,
+0xbfff0000,0x80000000,0x00000000,
+};
+#endif
+
+
+/* A[]: Stirling's formula expansion of log gamma
+ * B[], C[]: log gamma function between 2 and 3
+ */
+
+
+/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2)
+ * x >= 8
+ * Peak relative error 1.51e-21
+ * Relative spread of error peaks 5.67e-21
+ */
+#if UNK
+static long double A[7] = {
+ 4.885026142432270781165E-3L,
+-1.880801938119376907179E-3L,
+ 8.412723297322498080632E-4L,
+-5.952345851765688514613E-4L,
+ 7.936507795855070755671E-4L,
+-2.777777777750349603440E-3L,
+ 8.333333333333331447505E-2L,
+};
+#endif
+#if IBMPC
+static const unsigned short A[] = {
+0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD
+0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD
+0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD
+0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD
+0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD
+0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD
+0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD
+};
+#endif
+#if MIEEE
+static long A[21] = {
+0x3ff70000,0xa01291c2,0xcc08d984,
+0xbff50000,0xf6853da1,0x03043d91,
+0x3ff40000,0xdc88d492,0xaad13bdc,
+0xbff40000,0x9c09844e,0x9fce8b20,
+0x3ff40000,0xd00d0092,0x30e5f8f2,
+0xbff60000,0xb60b60b6,0x03a84d88,
+0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc,
+};
+#endif
+
+/* log gamma(x+2) = x B(x)/C(x)
+ * 0 <= x <= 1
+ * Peak relative error 7.16e-22
+ * Relative spread of error peaks 4.78e-20
+ */
+#if UNK
+static long double B[7] = {
+-2.163690827643812857640E3L,
+-8.723871522843511459790E4L,
+-1.104326814691464261197E6L,
+-6.111225012005214299996E6L,
+-1.625568062543700591014E7L,
+-2.003937418103815175475E7L,
+-8.875666783650703802159E6L,
+};
+static long double C[7] = {
+/* 1.000000000000000000000E0L,*/
+-5.139481484435370143617E2L,
+-3.403570840534304670537E4L,
+-6.227441164066219501697E5L,
+-4.814940379411882186630E6L,
+-1.785433287045078156959E7L,
+-3.138646407656182662088E7L,
+-2.099336717757895876142E7L,
+};
+#endif
+#if IBMPC
+static const unsigned short B[] = {
+0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD
+0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD
+0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD
+0x259a,0x258c,0xf206,0xba7f,0xc015, XPD
+0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD
+0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD
+0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD
+};
+static const unsigned short C[] = {
+/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
+0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD
+0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD
+0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD
+0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD
+0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD
+0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD
+0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD
+};
+#endif
+#if MIEEE
+static long B[21] = {
+0xc00a0000,0x873b0da1,0x49959557,
+0xc00f0000,0xaa635b8c,0x9af8fe44,
+0xc0130000,0x86ce3684,0x7cf55aa8,
+0xc0150000,0xba7ff206,0x258c259a,
+0xc0160000,0xf80ac0a0,0x1ca3be18,
+0xc0170000,0x98e36717,0x2c42168f,
+0xc0160000,0x876e92c8,0x9d552051,
+};
+static long C[21] = {
+/*0x3fff0000,0x80000000,0x00000000,*/
+0xc0080000,0x807cae76,0xcf2faa77,
+0xc00e0000,0x84f3b55a,0x0d74b280,
+0xc0120000,0x980981dc,0xcd30a505,
+0xc0150000,0x92f0b8c2,0x42463369,
+0xc0170000,0x8837be6f,0x6aee63cf,
+0xc0170000,0xef75b009,0xccc726bb,
+0xc0170000,0xa02aab96,0xbae8462b,
+};
+#endif
+
+/* log( sqrt( 2*pi ) ) */
+static const long double LS2PI  =  0.91893853320467274178L;
+#define MAXLGM 1.04848146839019521116e+4928L
+
+
+/* Logarithm of gamma function */
+/* Reentrant version */ 
+
+long double __lgammal_r(long double x, int* sgngaml)
+{
+long double p, q, w, z, f, nx;
+int i;
+
+*sgngaml = 1;
+#ifdef NANS
+if( isnanl(x) )
+	return(NANL);
+#endif
+#ifdef INFINITIES
+if( !isfinitel(x) )
+	return(INFINITYL);
+#endif
+if( x < -34.0L )
+	{
+	q = -x;
+	w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */
+	p = floorl(q);
+	if( p == q )
+		{
+lgsing:
+		_SET_ERRNO(EDOM);
+		mtherr( "lgammal", SING );
+#ifdef INFINITIES
+		return (INFINITYL);
+#else
+		return (MAXNUML);
+#endif
+		}
+	i = p;
+	if( (i & 1) == 0 )
+		*sgngaml = -1;
+	else
+		*sgngaml = 1;
+	z = q - p;
+	if( z > 0.5L )
+		{
+		p += 1.0L;
+		z = p - q;
+		}
+	z = q * sinl( PIL * z );
+	if( z == 0.0L )
+		goto lgsing;
+/*	z = LOGPI - logl( z ) - w; */
+	z = logl( PIL/z ) - w;
+	return( z );
+	}
+
+if( x < 13.0L )
+	{
+	z = 1.0L;
+	nx = floorl( x +  0.5L );
+	f = x - nx;
+	while( x >= 3.0L )
+		{
+		nx -= 1.0L;
+		x = nx + f;
+		z *= x;
+		}
+	while( x < 2.0L )
+		{
+		if( fabsl(x) <= 0.03125 )
+			goto lsmall;
+		z /= nx +  f;
+		nx += 1.0L;
+		x = nx + f;
+		}
+	if( z < 0.0L )
+		{
+		*sgngaml = -1;
+		z = -z;
+		}
+	else
+		*sgngaml = 1;
+	if( x == 2.0L )
+		return( logl(z) );
+	x = (nx - 2.0L) + f;
+	p = x * polevll( x, B, 6 ) / p1evll( x, C, 7);
+	return( logl(z) + p );
+	}
+
+if( x > MAXLGM )
+	{
+	_SET_ERRNO(ERANGE);
+	mtherr( "lgammal", OVERFLOW );
+#ifdef INFINITIES
+	return( *sgngaml * INFINITYL );
+#else
+	return( *sgngaml * MAXNUML );
+#endif
+	}
+
+q = ( x - 0.5L ) * logl(x) - x + LS2PI;
+if( x > 1.0e10L )
+	return(q);
+p = 1.0L/(x*x);
+q += polevll( p, A, 6 ) / x;
+return( q );
+
+
+lsmall:
+if( x == 0.0L )
+	goto lgsing;
+if( x < 0.0L )
+	{
+	x = -x;
+	q = z / (x * polevll( x, SN, 8 ));
+	}
+else
+	q = z / (x * polevll( x, S, 8 ));
+if( q < 0.0L )
+	{
+	*sgngaml = -1;
+	q = -q;
+	}
+else
+	*sgngaml = 1;
+q = logl( q );
+return(q);
+}
+
+/* This is the C99 version */
+
+long double lgammal(long double x)
+{
+  int local_sgngaml=0;
+  return (__lgammal_r(x, &local_sgngaml));
+}
@@ -0,0 +1,10 @@
+#include <math.h>
+
+long long llrint (double x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
@@ -0,0 +1,9 @@
+#include <math.h>
+
+long long llrintf (float x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
@@ -0,0 +1,10 @@
+#include <math.h>
+
+long long llrintl (long double x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
@@ -0,0 +1,38 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ *
+ * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
+ */
+
+	.file	"log.s"
+	.text
+	.align 4
+one:	.double 1.0
+	/* It is not important that this constant is precise.  It is only
+	   a value which is known to be on the safe side for using the
+	   fyl2xp1 instruction.  */
+limit:	.double 0.29
+
+	.text
+	.align 4
+.globl _log
+	.def	_log;	.scl	2;	.type	32;	.endef
+_log:
+	fldln2			/* log(2) */
+	fldl	4(%esp)		/* x : log(2) */
+	fld	%st		/* x : x : log(2) */
+	fsubl	one		/* x-1 : x : log(2) */
+	fld	%st		/* x-1 : x-1 : x : log(2) */
+	fabs			/* |x-1| : x-1 : x : log(2) */
+	fcompl	limit		/* x-1 : x : log(2) */
+	fnstsw			/* x-1 : x : log(2) */
+	andb	$0x45, %ah
+	jz	2f
+	fstp	%st(1)		/* x-1 : log(2) */
+	fyl2xp1			/* log(x) */
+	ret
+
+2:	fstp	%st(0)		/* x : log(2) */
+	fyl2x			/* log(x) */
+	ret
@@ -0,0 +1,48 @@
+/*
+ * Written by J.T. Conklin <jtc@netbsd.org>.
+ * Public domain.
+ * Adapted for float type by Ulrich Drepper <drepper@cygnus.com>.
+ *
+ * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
+ */
+
+	.file	"log10.s"
+	.text
+	.align 4
+one:	.double 1.0
+	/* It is not important that this constant is precise.  It is only
+	   a value which is known to be on the safe side for using the
+	   fyl2xp1 instruction.  */
+limit:	.double 0.29
+
+	.text
+	.align 4
+.globl _log10
+	.def	_log10;	.scl	2;	.type	32;	.endef
+_log10:
+	fldlg2			/* log10(2) */
+	fldl	4(%esp)		/* x : log10(2) */
+	fxam
+	fnstsw
+	fld	%st		/* x : x : log10(2) */
+	sahf
+	jc	3f		/* in case x is NaN or ±Inf */
+4:	fsubl	one		/* x-1 : x : log10(2) */
+	fld	%st		/* x-1 : x-1 : x : log10(2) */
+	fabs			/* |x-1| : x-1 : x : log10(2) */
+	fcompl	limit		/* x-1 : x : log10(2) */
+	fnstsw			/* x-1 : x : log10(2) */
+	andb	$0x45, %ah
+	jz	2f
+	fstp	%st(1)		/* x-1 : log10(2) */
+	fyl2xp1			/* log10(x) */
+	ret
+
+2:	fstp	%st(0)		/* x : log10(2) */
+	fyl2x			/* log10(x) */
+	ret
+
+3:	jp	4b		/* in case x is ±Inf */
+	fstp	%st(1)
+	fstp	%st(1)
+	ret
--- a/Show More
+++ b/Show More