enable complex in libm.h

[libm] / src / internal / libm.h

diff --git a/src/internal/libm.h b/src/internal/libm.h
index 13a4378..fecc50a 100644 (file)
--- a/src/internal/libm.h
+++ b/src/internal/libm.h
@@ -16,26 +16,11 @@
 #include <stdint.h>
 #include <float.h>
 #include <math.h>
-#if 0
 #include <complex.h>
-#endif
 
 // FIXME
 #include "ldhack.h"
 
-/*
- * The original fdlibm 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.
- */
-
 union fshape {
        float value;
        uint32_t bits;
@@ -167,7 +152,7 @@ long double __p1evll(long double, long double *, int);
  */
 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
 
-// TODO: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
+// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
 /*
  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
  */
@@ -181,66 +166,32 @@ void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
 #endif
 
 
-#if 0
-/*
- * C99 specifies that complex numbers have the same representation as
- * an array of two elements, where the first element is the real part
- * and the second element is the imaginary part.
- */
-typedef union {
-        float complex f;
-        float a[2];
-} float_complex;
-typedef union {
-        double complex f;
-        double a[2];
-} double_complex;
-typedef union {
-        long double complex f;
-        long double a[2];
-} long_double_complex;
-#define REALPART(z)     ((z).a[0])
-#define IMAGPART(z)     ((z).a[1])
+/* complex */
 
-/*
- * Inline functions that can be used to construct complex values.
- *
- * The C99 standard intends x+I*y to be used for this, but x+I*y is
- * currently unusable in general since gcc introduces many overflow,
- * underflow, sign and efficiency bugs by rewriting I*y as
- * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
- * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
- * to -0.0+I*0.0.
- */
-static inline float complex
-cpackf(float x, float y)
-{
-        float_complex z;
-
-        REALPART(z) = x;
-        IMAGPART(z) = y;
-        return (z.f);
-}
-
-static inline double complex
-cpack(double x, double y)
-{
-        double_complex z;
-
-        REALPART(z) = x;
-        IMAGPART(z) = y;
-        return (z.f);
-}
-
-static inline long double complex
-cpackl(long double x, long double y)
-{
-        long_double_complex z;
-
-        REALPART(z) = x;
-        IMAGPART(z) = y;
-        return (z.f);
-}
-#endif
+union dcomplex {
+       double complex z;
+       double a[2];
+};
+union fcomplex {
+       float complex z;
+       float a[2];
+};
+union lcomplex {
+       long double complex z;
+       long double a[2];
+};
+
+// FIXME: move to complex.h ?
+#define creal(z) ((double)(z))
+#define crealf(z) ((float)(z))
+#define creall(z) ((long double)(z))
+#define cimag(z) ((union dcomplex){(z)}.a[1])
+#define cimagf(z) ((union fcomplex){(z)}.a[1])
+#define cimagl(z) ((union lcomplex){(z)}.a[1])
+
+/* x + y*I is not supported properly by gcc */
+#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z)
+#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z)
+#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z)
 
 #endif