initial cmath code and minor libm.h update
[libm] / src / internal / libm.h
index 13a4378..00c2972 100644 (file)
 #include <stdint.h>
 #include <float.h>
 #include <math.h>
 #include <stdint.h>
 #include <float.h>
 #include <math.h>
-#if 0
 #include <complex.h>
 #include <complex.h>
-#endif
 
 // FIXME
 #include "ldhack.h"
 
 
 // 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;
 union fshape {
        float value;
        uint32_t bits;
@@ -132,6 +117,7 @@ do {                                                            \
 } while (0)
 
 /* fdlibm kernel functions */
 } while (0)
 
 /* fdlibm kernel functions */
+
 int    __rem_pio2_large(double*,double*,int,int,int);
 
 int    __rem_pio2(double,double*);
 int    __rem_pio2_large(double*,double*,int,int,int);
 
 int    __rem_pio2(double,double*);
@@ -139,20 +125,15 @@ double __sin(double,double,int);
 double __cos(double,double);
 double __tan(double,double,int);
 double __ldexp_exp(double,int);
 double __cos(double,double);
 double __tan(double,double,int);
 double __ldexp_exp(double,int);
-#if 0
 double complex __ldexp_cexp(double complex,int);
 double complex __ldexp_cexp(double complex,int);
-#endif
 
 int    __rem_pio2f(float,double*);
 float  __sindf(double);
 float  __cosdf(double);
 float  __tandf(double,int);
 float  __ldexp_expf(float,int);
 
 int    __rem_pio2f(float,double*);
 float  __sindf(double);
 float  __cosdf(double);
 float  __tandf(double,int);
 float  __ldexp_expf(float,int);
-#if 0
 float complex __ldexp_cexpf(float complex,int);
 float complex __ldexp_cexpf(float complex,int);
-#endif
 
 
-/* long double precision kernel functions */
 long double __sinl(long double, long double, int);
 long double __cosl(long double, long double);
 long double __tanl(long double, long double, int);
 long double __sinl(long double, long double, int);
 long double __cosl(long double, long double);
 long double __tanl(long double, long double, int);
@@ -161,13 +142,7 @@ long double __tanl(long double, long double, int);
 long double __polevll(long double, long double *, int);
 long double __p1evll(long double, long double *, int);
 
 long double __polevll(long double, long double *, int);
 long double __p1evll(long double, long double *, int);
 
-// FIXME: nan
-/*
- * Common routine to process the arguments to nan(), nanf(), and nanl().
- */
-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.
  */
 /*
  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
  */
@@ -181,66 +156,32 @@ void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
 #endif
 
 
 #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
 
 #endif