Add incomplet long double math support to libmingwex.a

1 files changed, 377 insertions, 113 deletions

diff --git a/winsup/mingw/include/math.h b/winsup/mingw/include/math.h
index c87b6163b..189e93bf8 100644
--- a/winsup/mingw/include/math.h
+++ b/winsup/mingw/include/math.h
@@ -144,17 +144,30 @@ double	log (double);
 double	log10 (double);
 double	pow (double, double);
 double	sqrt (double);
+extern __inline__  double sqrt (double x)
+{
+  double res;
+  asm ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
 double	ceil (double);
 double	floor (double);
 double	fabs (double);
+extern __inline__  double fabs (double x)
+{
+  double res;
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
 double	ldexp (double, int);
 double	frexp (double, int*);
 double	modf (double, double*);
 double	fmod (double, double);
 
+
 #ifndef	__STRICT_ANSI__
 
-/* Complex number (for cabs) */
+/* Complex number (for cabs). This really belongs in complex.h */
 struct _complex
 {
 	double	x;	/* Real part */
@@ -162,6 +175,7 @@ struct _complex
 };
 
 double	_cabs (struct _complex);
+
 double	_hypot (double, double);
 double	_j0 (double);
 double	_j1 (double);
@@ -190,17 +204,15 @@ int	_isnan		(double);
 
 /* END FLOAT.H COPY */
 
-#if !defined (_NO_OLDNAMES)  \
-   || (defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L )
 
 /*
- * Non-underscored versions of non-ANSI functions. These reside in
- * liboldnames.a. They are now also ISO C99 standand names.
- * Provided for extra portability.
+ * Non-underscored versions of non-ANSI functions.
+ * These reside in liboldnames.a.
  */
 
+#if !defined (_NO_OLDNAMES)
+
 double cabs (struct _complex);
-double hypot (double, double);
 double j0 (double);
 double j1 (double);
 double jn (int, double);
@@ -208,28 +220,27 @@ double y0 (double);
 double y1 (double);
 double yn (int, double);
 
-#endif	/* Not _NO_OLDNAMES */
+double chgsign	(double);
+double scalb (double, long);
+int finite (double);
+int fpclass (double);
 
-#endif	/* Not __STRICT_ANSI__ */
+#endif	/* Not _NO_OLDNAMES */
 
-#ifdef __cplusplus
-}
-#endif
-#endif	/* Not RC_INVOKED */
+#endif /* __STRICT_ANSI__ */
 
 
 #ifndef __NO_ISOCEXT
-#ifndef RC_INVOKED
-#ifdef __cplusplus
-extern "C" {
-#endif
-
 #if (defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L) \
 	|| !defined  __STRICT_ANSI__
 
-#define INFINITY HUGE_VAL
 #define NAN (0.0F/0.0F)
+#define HUGE_VALF (1.0F/0.0F)
+#define HUGE_VALL (1.0L/0.0L)
+#define INFINITY (1.0F/0.0F)
+
 
+/* 7.12.3.1 */
 /*
    Return values for fpclassify.
    These are based on Intel x87 fpu condition codes
@@ -244,17 +255,6 @@ extern "C" {
 /* 0x0200 is signbit mask */
 
 
-/* Return a  NaN */
-double nan(const char *tagp);
-float nanf(const char *tagp);
-long double nanl(const char *tagp);
-
-#ifndef __STRICT_ANSI__
-#define nan() nan("")
-#define nanf() nanf("")
-#define nanl() nanl("")
-#endif
-
 /*
   We can't inline float or double, because we want to ensure truncation
   to semantic type before classification. 
@@ -275,8 +275,16 @@ extern __inline__ int __fpclassifyl (long double x){
 		       : sizeof (x) == sizeof (double) ? __fpclassify (x) \
 		       : __fpclassifyl (x))
 
+/* 7.12.3.2 */
+#define isfinite(x) ((fpclassify(x) & FP_NAN) == 0)
+
+/* 7.12.3.3 */
+#define isinf(x) (fpclassify(x) == FP_INFINITE)
+
+/* 7.12.3.4 */
 /* We don't need to worry about trucation here:
    A NaN stays a NaN. */
+
 extern __inline__ int __isnan (double _x)
 {
   unsigned short sw;
@@ -304,15 +312,15 @@ extern __inline__ int __isnanl (long double _x)
     == FP_NAN;
 }
 
+
 #define isnan(x) (sizeof (x) == sizeof (float) ? __isnanf (x)	\
 		  : sizeof (x) == sizeof (double) ? __isnan (x)	\
 		  : __isnanl (x))
 
-#define isfinite(x) ((fpclassify(x) & FP_NAN) == 0)
-#define isinf(x) (fpclassify(x) == FP_INFINITE)
+/* 7.12.3.5 */
 #define isnormal(x) (fpclassify(x) == FP_NORMAL)
 
-
+/* 7.12.3.6 The signbit macro */
 extern __inline__ int __signbit (double x) {
   unsigned short stw;
   __asm__ ( "fxam; fstsw %%ax;": "=a" (stw) : "t" (x));
@@ -335,151 +343,407 @@ extern  __inline__ int __signbitl (long double x) {
 		    : sizeof (x) == sizeof (double) ? __signbit (x)	\
 		    : __signbitl (x))
 
-/* 
- *  With these functions, comparisons involving quiet NaNs set the FP
- *  condition code to "unordered".  The IEEE floating-point spec
- *  dictates that the result of floating-point comparisons should be
- *  false whenever a NaN is involved, with the exception of the !=, 
- *  which always returns true: yes, (NaN != NaN) is true).
- */
+/* 7.12.4 Trigonometric functions: Double in C89 */
+extern float sinf (float);
+extern long double sinl (long double);
 
-#if __GNUC__ >= 3
+extern float cosf (float);
+extern long double cosl (long double);
 
-#define isgreater(x, y) __builtin_isgreater(x, y)
-#define isgreaterequal(x, y) __builtin_isgreaterequal(x, y)
-#define isless(x, y) __builtin_isless(x, y)
-#define islessequal(x, y) __builtin_islessequal(x, y)
-#define islessgreater(x, y) __builtin_islessgreater(x, y)
-#define isunordered(x, y) __builtin_isunordered(x, y)
+extern float tanf (float);
+extern long double tanl (long double);
 
-#else
-/*  helper  */
-extern  __inline__ 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;
+extern float asinf (float);
+extern long double asinl (long double);
+
+extern float acosf (float);
+extern long double acosl (long double);
+
+extern float atanf (float);
+extern long double atanl (long double);
+
+extern float atan2f (float, float);
+extern long double atan2l (long double, long double);
+
+/* 7.12.5 Hyperbolic functions: Double in C89  */
+extern __inline__ float sinhf (float x)
+  {return (float) sinh (x);}
+extern long double sinhl (long double);
+
+extern __inline__ float coshf (float x)
+  {return (float) cosh (x);}
+extern long double coshl (long double);
+
+extern __inline__ float tanhf (float x)
+  {return (float) tanh (x);}
+extern long double tanhl (long double);
+
+/*
+ * TODO: asinh, acosh, atanh
+ */ 
+
+/* 7.12.6.1 Double in C89 */
+extern __inline__ float expf (float x)
+  {return (float) exp (x);}
+extern long double expl (long double);
+
+/* 7.12.6.2 */
+extern double exp2(double);
+extern float exp2f(float);
+extern long double exp2l(long double);
+
+/* 7.12.6.3 The expm1 functions: TODO */
+
+/* 7.12.6.4 Double in C89 */
+extern __inline__ float frexpf (float x, int* expn)
+  {return (float) frexp (x, expn);}
+extern long double frexpl (long double, int*);
+
+/* 7.12.6.5 */
+#define FP_ILOGB0 ((int)0x80000000)
+#define FP_ILOGBNAN ((int)0x80000000)
+extern int ilogb (double);
+extern int ilogbf (float);
+extern int ilogbl (long double);
+
+/* 7.12.6.6  Double in C89 */
+extern __inline__ float ldexpf (float x, int expn)
+  {return (float) ldexp (x, expn);}
+extern long double ldexpl (long double, int);
+
+/* 7.12.6.7 Double in C89 */
+extern float logf (float);
+extern long double logl (long double);
+
+/* 7.12.6.8 Double in C89 */
+extern float log10f (float);
+extern long double log10l (long double);
+
+/* 7.12.6.9 */
+extern double log1p(double);
+extern float log1pf(float);
+extern long double log1pl(long double);
+
+/* 7.12.6.10 */
+extern double log2 (double);
+extern float log2f (float);
+extern long double log2l (long double);
+
+/* 7.12.6.11 */
+extern double logb (double);
+extern float logbf (float);
+extern long double logbl (long double);
+
+extern __inline__ double logb (double x)
+{
+  double res;
+  asm ("fxtract\n\t"
+       "fstp	%%st" : "=t" (res) : "0" (x));
+  return res;
 }
 
-#define isgreater(x, y) ((__fp_unordered_compare(x, y) \
-			   & 0x4500) == 0)
-#define isless(x, y) ((__fp_unordered_compare (y, x) \
-                       & 0x4500) == 0)
-#define isgreaterequal(x, y) ((__fp_unordered_compare (x, y) \
-                               & FP_INFINITE) == 0)
-#define islessequal(x, y) ((__fp_unordered_compare(y, x) \
-			    & FP_INFINITE) == 0)
-#define islessgreater(x, y) ((__fp_unordered_compare(x, y) \
-			      & FP_SUBNORMAL) == 0)
-#define isunordered(x, y) ((__fp_unordered_compare(x, y) \
-			    & 0x4500) == 0x4500)
+extern __inline__ float logbf (float x)
+{
+  float res;
+  asm ("fxtract\n\t"
+       "fstp	%%st" : "=t" (res) : "0" (x));
+  return res;
+}
 
-#endif
+extern __inline__ long double logbl (long double x)
+{
+  long double res;
+  asm ("fxtract\n\t"
+       "fstp	%%st" : "=t" (res) : "0" (x));
+  return res;
+}
+
+/* 7.12.6.12  Double in C89 */
+extern float modff (float, float*);
+extern long double modfl (long double, long double*);
+
+/* 7.12.6.13 */
+extern double scalbn (double, int);
+extern float scalbnf (float, int);
+extern long double scalbnl (long double, int);
+
+extern double scalbln (double, long);
+extern float scalblnf (float, long);
+extern long double scalblnl (long double, long);
+
+/* 7.12.7.1 */
+/* Implementations adapted from Cephes versions */ 
+extern double cbrt (double);
+extern float cbrtf (float);
+extern long double cbrtl (long double);
+
+/* 7.12.7.2 The fabs functions: Double in C89 */
+extern __inline__ float fabsf (float x)
+{
+  float res;
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
+
+extern __inline__ long double fabsl (long double x)
+{
+  long double res;
+  asm ("fabs;" : "=t" (res) : "0" (x));
+  return res;
+}
+
+/* 7.12.7.3  */
+extern double hypot (double, double); /* in libmoldname.a */
+extern __inline__ float hypotf (float x, float y)
+  { return (float) _hypot (x, y);}
+extern long double hypotl (long double, long double);
+
+/* 7.12.7.4 The pow functions. Double in C89 */
+extern __inline__ float powf (float x, float y)
+  {return (float) pow (x, y);}
+extern long double powl (long double, long double);
 
+/* 7.12.7.5 The sqrt functions. Double in C89. */
+extern __inline__ float sqrtf (float x)
+{
+  float res;
+  asm ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
+
+extern __inline__ long double sqrtl (long double x)
+{
+  long double res;
+  asm ("fsqrt" : "=t" (res) : "0" (x));
+  return res;
+}
+
+/* 7.12.8 Error and gamma functions: TODO */
+
+/* 7.12.9.1 Double in C89 */
+extern float ceilf (float);
+extern long double ceill (long double);
+
+/* 7.12.9.2 Double in C89 */
+extern float floorf (float);
+extern long double floorl (long double);
+
+/* 7.12.9.3 */
+extern double nearbyint ( double);
+extern float nearbyintf (float);
+extern long double nearbyintl (long double);
+
+/* 7.12.9.4 */
 /* round, using fpu control word settings */
-extern  __inline__ double rint (double x)
+extern __inline__ double rint (double x)
 {
   double retval;
   __asm__ ("frndint;": "=t" (retval) : "0" (x));
   return retval;
 }
 
-extern  __inline__ float rintf (float x)
+extern __inline__ float rintf (float x)
 {
   float retval;
   __asm__ ("frndint;" : "=t" (retval) : "0" (x) );
   return retval;
 }
 
-extern  __inline__ long double rintl (long double x)
+extern __inline__ long double rintl (long double x)
 {
   long double retval;
   __asm__ ("frndint;" : "=t" (retval) : "0" (x) );
   return retval;
 }
 
+/* 7.12.9.5 */
+extern __inline__ long lrint (double x) 
+{
+  long retval;  
+  __asm__ __volatile__							      \
+    ("fistpl %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+extern __inline__ long lrintf (float x) 
+{
+  long retval;
+  __asm__ __volatile__							      \
+    ("fistpl %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+extern __inline__ long lrintl (long double x) 
+{
+  long retval;
+  __asm__ __volatile__							      \
+    ("fistpl %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+extern __inline__ long long llrint (double x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+extern __inline__ long long llrintf (float x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+extern __inline__ long long llrintl (long double x) 
+{
+  long long retval;
+  __asm__ __volatile__							      \
+    ("fistpll %0"  : "=m" (retval) : "t" (x) : "st");				      \
+  return retval;
+}
+
+/* 7.12.9.6 */
 /* round away from zero, regardless of fpu control word settings */
 extern double round (double);
 extern float roundf (float);
 extern long double roundl (long double);
 
+/* 7.12.9.7  */
+extern long lround (double);
+extern long lroundf (float);
+extern long lroundl (long double);
+
+extern long long llround (double);
+extern long long llroundf (float);
+extern long long llroundl (long double);
 
+/* 7.12.9.8 */
 /* round towards zero, regardless of fpu control word settings */
 extern double trunc (double);
 extern float truncf (float);
 extern long double truncl (long double);
 
+/* 7.12.10.1 Double in C89 */
+extern float fmodf (float, float);
+extern long double fmodl (long double, long double);
+
+/* 7.12.10.2 */ 
+extern double remainder (double, double);
+extern float remainderf (float, float);
+extern long double remainderl (long double, long double);
+
+/* 7.12.10.3 */
+extern double remquo(double, double, int *);
+extern float remquof(float, float, int *);
+extern long double remquol(long double, long double, int *);
+
+/* 7.12.11.1 */
+extern double copysign (double, double); /* in libmoldname.a */
+extern float copysignf (float, float);
+extern long double copysignl (long double, long double);
+
+/* 7.12.11.2 Return a NaN */
+extern double nan(const char *tagp);
+extern float nanf(const char *tagp);
+extern long double nanl(const char *tagp);
+
+#ifndef __STRICT_ANSI__
+#define _nan() nan("")
+#define _nanf() nanf("")
+#define _nanl() nanl("")
+#endif
+
+/* 7.12.11.3 */
+extern double nextafter (double, double); /* in libmoldname.a */
+extern float nextafterf (float, float);
+/* TODO: Not yet implemented */
+/* extern long double nextafterl (long double, long double); */
+
+/* 7.12.11.4 The nexttoward functions: TODO */
+
+/* 7.12.12.1 */
+/*  x > y ? (x - y) : 0.0  */
+extern double fdim (double x, double y);
+extern float fdimf (float x, float y);
+extern long double fdiml (long double x, long double y);
 
 /* fmax and fmin.
    NaN arguments are treated as missing data: if one argument is a NaN
    and the other numeric, then these functions choose the numeric
    value. */
 
+/* 7.12.12.2 */
 extern double fmax  (double, double);
 extern float fmaxf (float, float);
 extern long double fmaxl (long double, long double);
 
+/* 7.12.12.3 */
 extern double fmin (double, double);
 extern float fminf (float, float);
 extern long double fminl (long double, long double);
 
+/* 7.12.13.1 */
 /* return x * y + z as a ternary op */ 
 extern double fma (double, double, double);
 extern float fmaf (float, float, float);
 extern long double fmal (long double, long double, long double);
 
-/*  x > y ? (x - y) : 0.0  */
-extern double fdim (double, double);
-extern float fdimf (float, float);
-extern long double fdiml (long double, long double);
-
-/* one lonely transcendental */
-extern double log2 (double _x);
-extern float log2f (float _x);
-extern long double log2l (long double _x);
-#endif /* __STDC_VERSION__ >= 199901L */
 
+/* 7.12.14 */
+/* 
+ *  With these functions, comparisons involving quiet NaNs set the FP
+ *  condition code to "unordered".  The IEEE floating-point spec
+ *  dictates that the result of floating-point comparisons should be
+ *  false whenever a NaN is involved, with the exception of the != op, 
+ *  which always returns true: yes, (NaN != NaN) is true).
+ */
 
-/* The underscored versions for double are in MSVCRT.dll.
-   The stubs for float and double versions are in libmingwex.a */
+#if __GNUC__ >= 3
 
-double copysign (double, double);
-float copysignf (float, float);
-long double copysignl (long double, long double);
+#define isgreater(x, y) __builtin_isgreater(x, y)
+#define isgreaterequal(x, y) __builtin_isgreaterequal(x, y)
+#define isless(x, y) __builtin_isless(x, y)
+#define islessequal(x, y) __builtin_islessequal(x, y)
+#define islessgreater(x, y) __builtin_islessgreater(x, y)
+#define isunordered(x, y) __builtin_isunordered(x, y)
+
+#else
+/*  helper  */
+extern  __inline__ 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;
+}
+
+#define isgreater(x, y) ((__fp_unordered_compare(x, y) \
+			   & 0x4500) == 0)
+#define isless(x, y) ((__fp_unordered_compare (y, x) \
+                       & 0x4500) == 0)
+#define isgreaterequal(x, y) ((__fp_unordered_compare (x, y) \
+                               & FP_INFINITE) == 0)
+#define islessequal(x, y) ((__fp_unordered_compare(y, x) \
+			    & FP_INFINITE) == 0)
+#define islessgreater(x, y) ((__fp_unordered_compare(x, y) \
+			      & FP_SUBNORMAL) == 0)
+#define isunordered(x, y) ((__fp_unordered_compare(x, y) \
+			    & 0x4500) == 0x4500)
+
+#endif
+
+
+#endif /* __STDC_VERSION__ >= 199901L */
+#endif /* __NO_ISOCEXT */
 
-double logb (double);
-float logbf (float);
-double nextafter (double, double);
-float nextafterf (float, float);
-double scalb (double, long);
-float scalbf (float, long);
-
-#if !defined (__STRICT_ANSI__)  /* inline using non-ANSI functions */
-extern  __inline__ double copysign (double x, double y)
-	{ return _copysign(x, y); }
-extern  __inline__ float copysignf (float x, float y)
-	{ return  _copysign(x, y); } 
-extern  __inline__ double logb (double x)
-	{ return _logb(x); }
-extern  __inline__ float logbf (float x)
-	{ return  _logb(x); }
-extern  __inline__ double nextafter(double x, double y)
-	{ return _nextafter(x, y); }
-extern  __inline__ float nextafterf(float x, float y)
-	{ return _nextafter(x, y); }
-extern  __inline__ double scalb (double x, long i)
-	{ return _scalb (x, i); }
-extern  __inline__ float scalbf (float x, long i)
-	{ return _scalb(x, i); }
-#endif /* (__STRICT_ANSI__)  */
 
 #ifdef __cplusplus
 }
 #endif
 #endif	/* Not RC_INVOKED */
 
-#endif /* __NO_ISOCEXT */
 
 #endif	/* Not _MATH_H_ */
-