* include/math.h (ashinh, asinhf, asinhl, acosh, acoshf, acoshl,

	atanh, atanhf, atanhl): Add prototypes.
	* mingwex/Makefile.in (MATH_OBJS): Add objects for above to list.
	(MATH_DISTFILES): Add sources for above and fastmath.h to list.
	Specify dependency on fastmath.h for new objects.
	* mingwex/math/fastmath.h: New file.
	* mingwex/math/ashinh.c: New file.
	* mingwex/math/asinhf.c: New file.
	* mingwex/math/asinhl.c: New file.
	* mingwex/math/acosh.c: New file.
	* mingwex/math/acoshf.c: New file.
	* mingwex/math/acoshl.c: New file.
	* mingwex/math/atanh.c: New file.
	* mingwex/math/atanhf.c: New file.
	* mingwex/math/atanhl.c: New file.

10 files changed, 367 insertions, 0 deletions

diff --git a/winsup/mingw/mingwex/math/acosh.c b/winsup/mingw/mingwex/math/acosh.c
new file mode 100755
index 000000000..1497883cf
--- /dev/null
+++ b/winsup/mingw/mingwex/math/acosh.c
@@ -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)));
+}
diff --git a/winsup/mingw/mingwex/math/acoshf.c b/winsup/mingw/mingwex/math/acoshf.c
new file mode 100755
index 000000000..08f190fcb
--- /dev/null
+++ b/winsup/mingw/mingwex/math/acoshf.c
@@ -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)));
+}
diff --git a/winsup/mingw/mingwex/math/acoshl.c b/winsup/mingw/mingwex/math/acoshl.c
new file mode 100755
index 000000000..c461176bb
--- /dev/null
+++ b/winsup/mingw/mingwex/math/acoshl.c
@@ -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)));
+}
diff --git a/winsup/mingw/mingwex/math/asinh.c b/winsup/mingw/mingwex/math/asinh.c
new file mode 100755
index 000000000..30404497d
--- /dev/null
+++ b/winsup/mingw/mingwex/math/asinh.c
@@ -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);
+}
+
diff --git a/winsup/mingw/mingwex/math/asinhf.c b/winsup/mingw/mingwex/math/asinhf.c
new file mode 100755
index 000000000..080a9278d
--- /dev/null
+++ b/winsup/mingw/mingwex/math/asinhf.c
@@ -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);
+}
diff --git a/winsup/mingw/mingwex/math/asinhl.c b/winsup/mingw/mingwex/math/asinhl.c
new file mode 100755
index 000000000..8f027e83d
--- /dev/null
+++ b/winsup/mingw/mingwex/math/asinhl.c
@@ -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);
+}
diff --git a/winsup/mingw/mingwex/math/atanh.c b/winsup/mingw/mingwex/math/atanh.c
new file mode 100755
index 000000000..b5d9ce100
--- /dev/null
+++ b/winsup/mingw/mingwex/math/atanh.c
@@ -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;
+}
diff --git a/winsup/mingw/mingwex/math/atanhf.c b/winsup/mingw/mingwex/math/atanhf.c
new file mode 100755
index 000000000..b7c30823e
--- /dev/null
+++ b/winsup/mingw/mingwex/math/atanhf.c
@@ -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;
+}
diff --git a/winsup/mingw/mingwex/math/atanhl.c b/winsup/mingw/mingwex/math/atanhl.c
new file mode 100755
index 000000000..2d5fec02a
--- /dev/null
+++ b/winsup/mingw/mingwex/math/atanhl.c
@@ -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;
+}
diff --git a/winsup/mingw/mingwex/math/fastmath.h b/winsup/mingw/mingwex/math/fastmath.h
new file mode 100755
index 000000000..01b06b3eb
--- /dev/null
+++ b/winsup/mingw/mingwex/math/fastmath.h
@@ -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