nsz Git - libm/blob - src/internal/libm.h

   1 /* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 #ifndef _LIBM_H
  14 #define _LIBM_H
  15
  16 #include <stdint.h>
  17 #include <float.h>
  18 #include <math.h>
  19 #if 0
  20 #include <complex.h>
  21 #endif
  22
  23 // FIXME
  24 #include "ldhack.h"
  25
  26 /*
  27  * The original fdlibm code used statements like:
  28  *      n0 = ((*(int*)&one)>>29)^1;             * index of high word *
  29  *      ix0 = *(n0+(int*)&x);                   * high word of x *
  30  *      ix1 = *((1-n0)+(int*)&x);               * low word of x *
  31  * to dig two 32 bit words out of the 64 bit IEEE floating point
  32  * value.  That is non-ANSI, and, moreover, the gcc instruction
  33  * scheduler gets it wrong.  We instead use the following macros.
  34  * Unlike the original code, we determine the endianness at compile
  35  * time, not at run time; I don't see much benefit to selecting
  36  * endianness at run time.
  37  */
  38
  39 union fshape {
  40         float value;
  41         uint32_t bits;
  42 };
  43
  44 union dshape {
  45         double value;
  46         uint64_t bits;
  47 };
  48
  49 /* Get two 32 bit ints from a double.  */
  50 #define EXTRACT_WORDS(hi,lo,d)                                  \
  51 do {                                                            \
  52   union dshape __u;                                             \
  53   __u.value = (d);                                              \
  54   (hi) = __u.bits >> 32;                                        \
  55   (lo) = (uint32_t)__u.bits;                                    \
  56 } while (0)
  57
  58 /* Get a 64 bit int from a double.  */
  59 #define EXTRACT_WORD64(i,d)                                     \
  60 do {                                                            \
  61   union dshape __u;                                             \
  62   __u.value = (d);                                              \
  63   (i) = __u.bits;                                               \
  64 } while (0)
  65
  66 /* Get the more significant 32 bit int from a double.  */
  67 #define GET_HIGH_WORD(i,d)                                      \
  68 do {                                                            \
  69   union dshape __u;                                             \
  70   __u.value = (d);                                              \
  71   (i) = __u.bits >> 32;                                         \
  72 } while (0)
  73
  74 /* Get the less significant 32 bit int from a double.  */
  75 #define GET_LOW_WORD(i,d)                                       \
  76 do {                                                            \
  77   union dshape __u;                                             \
  78   __u.value = (d);                                              \
  79   (i) = (uint32_t)__u.bits;                                     \
  80 } while (0)
  81
  82 /* Set a double from two 32 bit ints.  */
  83 #define INSERT_WORDS(d,hi,lo)                                   \
  84 do {                                                            \
  85   union dshape __u;                                             \
  86   __u.bits = ((uint64_t)(hi) << 32) | (uint32_t)(lo);           \
  87   (d) = __u.value;                                              \
  88 } while (0)
  89
  90 /* Set a double from a 64 bit int.  */
  91 #define INSERT_WORD64(d,i)                                      \
  92 do {                                                            \
  93   union dshape __u;                                             \
  94   __u.bits = (i);                                               \
  95   (d) = __u.value;                                              \
  96 } while (0)
  97
  98 /* Set the more significant 32 bits of a double from an int.  */
  99 #define SET_HIGH_WORD(d,hi)                                     \
 100 do {                                                            \
 101   union dshape __u;                                             \
 102   __u.value = (d);                                              \
 103   __u.bits &= 0xffffffff;                                       \
 104   __u.bits |= (uint64_t)(hi) << 32;                             \
 105   (d) = __u.value;                                              \
 106 } while (0)
 107
 108 /* Set the less significant 32 bits of a double from an int.  */
 109 #define SET_LOW_WORD(d,lo)                                      \
 110 do {                                                            \
 111   union dshape __u;                                             \
 112   __u.value = (d);                                              \
 113   __u.bits &= 0xffffffff00000000ull;                            \
 114   __u.bits |= (uint32_t)(lo);                                   \
 115   (d) = __u.value;                                              \
 116 } while (0)
 117
 118 /* Get a 32 bit int from a float.  */
 119 #define GET_FLOAT_WORD(i,d)                                     \
 120 do {                                                            \
 121   union fshape __u;                                             \
 122   __u.value = (d);                                              \
 123   (i) = __u.bits;                                               \
 124 } while (0)
 125
 126 /* Set a float from a 32 bit int.  */
 127 #define SET_FLOAT_WORD(d,i)                                     \
 128 do {                                                            \
 129   union fshape __u;                                             \
 130   __u.bits = (i);                                               \
 131   (d) = __u.value;                                              \
 132 } while (0)
 133
 134 /* fdlibm kernel functions */
 135 int    __rem_pio2_slow(double*,double*,int,int,int);
 136
 137 int    __rem_pio2(double,double*);
 138 double __sin(double,double,int);
 139 double __cos(double,double);
 140 double __tan(double,double,int);
 141 double __ldexp_exp(double,int);
 142 #if 0
 143 double complex __ldexp_cexp(double complex,int);
 144 #endif
 145
 146 int    __rem_pio2f(float,double*);
 147 float  __sindf(double);
 148 float  __cosdf(double);
 149 float  __tandf(double,int);
 150 float  __ldexp_expf(float,int);
 151 #if 0
 152 float complex __ldexp_cexpf(float complex,int);
 153 #endif
 154
 155 /* long double precision kernel functions */
 156 long double __sinl(long double, long double, int);
 157 long double __cosl(long double, long double);
 158 long double __tanl(long double, long double, int);
 159
 160 /* polynomial evaluation */
 161 long double __polevll(long double, long double *, int);
 162 long double __p1evll(long double, long double *, int);
 163
 164 // FIXME: nan
 165 /*
 166  * Common routine to process the arguments to nan(), nanf(), and nanl().
 167  */
 168 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
 169
 170 // TODO: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
 171 /*
 172  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
 173  */
 174 #if 1
 175 #define STRICT_ASSIGN(type, lval, rval) do {    \
 176         volatile type __v = (rval);             \
 177         (lval) = __v;                           \
 178 } while (0)
 179 #else
 180 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
 181 #endif
 182
 183
 184 #if 0
 185 /*
 186  * C99 specifies that complex numbers have the same representation as
 187  * an array of two elements, where the first element is the real part
 188  * and the second element is the imaginary part.
 189  */
 190 typedef union {
 191         float complex f;
 192         float a[2];
 193 } float_complex;
 194 typedef union {
 195         double complex f;
 196         double a[2];
 197 } double_complex;
 198 typedef union {
 199         long double complex f;
 200         long double a[2];
 201 } long_double_complex;
 202 #define REALPART(z)     ((z).a[0])
 203 #define IMAGPART(z)     ((z).a[1])
 204
 205 /*
 206  * Inline functions that can be used to construct complex values.
 207  *
 208  * The C99 standard intends x+I*y to be used for this, but x+I*y is
 209  * currently unusable in general since gcc introduces many overflow,
 210  * underflow, sign and efficiency bugs by rewriting I*y as
 211  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
 212  * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
 213  * to -0.0+I*0.0.
 214  */
 215 static inline float complex
 216 cpackf(float x, float y)
 217 {
 218         float_complex z;
 219
 220         REALPART(z) = x;
 221         IMAGPART(z) = y;
 222         return (z.f);
 223 }
 224
 225 static inline double complex
 226 cpack(double x, double y)
 227 {
 228         double_complex z;
 229
 230         REALPART(z) = x;
 231         IMAGPART(z) = y;
 232         return (z.f);
 233 }
 234
 235 static inline long double complex
 236 cpackl(long double x, long double y)
 237 {
 238         long_double_complex z;
 239
 240         REALPART(z) = x;
 241         IMAGPART(z) = y;
 242         return (z.f);
 243 }
 244 #endif
 245
 246 #endif