lib/msun/src/math_private.h

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /*
  13  * from: @(#)fdlibm.h 5.1 93/09/24
  14  * $FreeBSD$
  15  */
  16
  17 #ifndef _MATH_PRIVATE_H_
  18 #define _MATH_PRIVATE_H_
  19
  20 #include <sys/types.h>
  21 #include <machine/endian.h>
  22
  23 /*
  24  * The original fdlibm code used statements like:
  25  *      n0 = ((*(int*)&one)>>29)^1;             * index of high word *
  26  *      ix0 = *(n0+(int*)&x);                   * high word of x *
  27  *      ix1 = *((1-n0)+(int*)&x);               * low word of x *
  28  * to dig two 32 bit words out of the 64 bit IEEE floating point
  29  * value.  That is non-ANSI, and, moreover, the gcc instruction
  30  * scheduler gets it wrong.  We instead use the following macros.
  31  * Unlike the original code, we determine the endianness at compile
  32  * time, not at run time; I don't see much benefit to selecting
  33  * endianness at run time.
  34  */
  35
  36 /*
  37  * A union which permits us to convert between a double and two 32 bit
  38  * ints.
  39  */
  40
  41 #ifdef __arm__
  42 #if defined(__VFP_FP__)
  43 #define IEEE_WORD_ORDER BYTE_ORDER
  44 #else
  45 #define IEEE_WORD_ORDER BIG_ENDIAN
  46 #endif
  47 #else /* __arm__ */
  48 #define IEEE_WORD_ORDER BYTE_ORDER
  49 #endif
  50
  51 #if IEEE_WORD_ORDER == BIG_ENDIAN
  52
  53 typedef union
  54 {
  55   double value;
  56   struct
  57   {
  58     u_int32_t msw;
  59     u_int32_t lsw;
  60   } parts;
  61   struct
  62   {
  63     u_int64_t w;
  64   } xparts;
  65 } ieee_double_shape_type;
  66
  67 #endif
  68
  69 #if IEEE_WORD_ORDER == LITTLE_ENDIAN
  70
  71 typedef union
  72 {
  73   double value;
  74   struct
  75   {
  76     u_int32_t lsw;
  77     u_int32_t msw;
  78   } parts;
  79   struct
  80   {
  81     u_int64_t w;
  82   } xparts;
  83 } ieee_double_shape_type;
  84
  85 #endif
  86
  87 /* Get two 32 bit ints from a double.  */
  88
  89 #define EXTRACT_WORDS(ix0,ix1,d)                                \
  90 do {                                                            \
  91   ieee_double_shape_type ew_u;                                  \
  92   ew_u.value = (d);                                             \
  93   (ix0) = ew_u.parts.msw;                                       \
  94   (ix1) = ew_u.parts.lsw;                                       \
  95 } while (0)
  96
  97 /* Get a 64-bit int from a double. */
  98 #define EXTRACT_WORD64(ix,d)                                    \
  99 do {                                                            \
 100   ieee_double_shape_type ew_u;                                  \
 101   ew_u.value = (d);                                             \
 102   (ix) = ew_u.xparts.w;                                         \
 103 } while (0)
 104
 105 /* Get the more significant 32 bit int from a double.  */
 106
 107 #define GET_HIGH_WORD(i,d)                                      \
 108 do {                                                            \
 109   ieee_double_shape_type gh_u;                                  \
 110   gh_u.value = (d);                                             \
 111   (i) = gh_u.parts.msw;                                         \
 112 } while (0)
 113
 114 /* Get the less significant 32 bit int from a double.  */
 115
 116 #define GET_LOW_WORD(i,d)                                       \
 117 do {                                                            \
 118   ieee_double_shape_type gl_u;                                  \
 119   gl_u.value = (d);                                             \
 120   (i) = gl_u.parts.lsw;                                         \
 121 } while (0)
 122
 123 /* Set a double from two 32 bit ints.  */
 124
 125 #define INSERT_WORDS(d,ix0,ix1)                                 \
 126 do {                                                            \
 127   ieee_double_shape_type iw_u;                                  \
 128   iw_u.parts.msw = (ix0);                                       \
 129   iw_u.parts.lsw = (ix1);                                       \
 130   (d) = iw_u.value;                                             \
 131 } while (0)
 132
 133 /* Set a double from a 64-bit int. */
 134 #define INSERT_WORD64(d,ix)                                     \
 135 do {                                                            \
 136   ieee_double_shape_type iw_u;                                  \
 137   iw_u.xparts.w = (ix);                                         \
 138   (d) = iw_u.value;                                             \
 139 } while (0)
 140
 141 /* Set the more significant 32 bits of a double from an int.  */
 142
 143 #define SET_HIGH_WORD(d,v)                                      \
 144 do {                                                            \
 145   ieee_double_shape_type sh_u;                                  \
 146   sh_u.value = (d);                                             \
 147   sh_u.parts.msw = (v);                                         \
 148   (d) = sh_u.value;                                             \
 149 } while (0)
 150
 151 /* Set the less significant 32 bits of a double from an int.  */
 152
 153 #define SET_LOW_WORD(d,v)                                       \
 154 do {                                                            \
 155   ieee_double_shape_type sl_u;                                  \
 156   sl_u.value = (d);                                             \
 157   sl_u.parts.lsw = (v);                                         \
 158   (d) = sl_u.value;                                             \
 159 } while (0)
 160
 161 /*
 162  * A union which permits us to convert between a float and a 32 bit
 163  * int.
 164  */
 165
 166 typedef union
 167 {
 168   float value;
 169   /* FIXME: Assumes 32 bit int.  */
 170   unsigned int word;
 171 } ieee_float_shape_type;
 172
 173 /* Get a 32 bit int from a float.  */
 174
 175 #define GET_FLOAT_WORD(i,d)                                     \
 176 do {                                                            \
 177   ieee_float_shape_type gf_u;                                   \
 178   gf_u.value = (d);                                             \
 179   (i) = gf_u.word;                                              \
 180 } while (0)
 181
 182 /* Set a float from a 32 bit int.  */
 183
 184 #define SET_FLOAT_WORD(d,i)                                     \
 185 do {                                                            \
 186   ieee_float_shape_type sf_u;                                   \
 187   sf_u.word = (i);                                              \
 188   (d) = sf_u.value;                                             \
 189 } while (0)
 190
 191 /* Get expsign as a 16 bit int from a long double.  */
 192
 193 #define GET_LDBL_EXPSIGN(i,d)                                   \
 194 do {                                                            \
 195   union IEEEl2bits ge_u;                                        \
 196   ge_u.e = (d);                                                 \
 197   (i) = ge_u.xbits.expsign;                                     \
 198 } while (0)
 199
 200 /* Set expsign of a long double from a 16 bit int.  */
 201
 202 #define SET_LDBL_EXPSIGN(d,v)                                   \
 203 do {                                                            \
 204   union IEEEl2bits se_u;                                        \
 205   se_u.e = (d);                                                 \
 206   se_u.xbits.expsign = (v);                                     \
 207   (d) = se_u.e;                                                 \
 208 } while (0)
 209
 210 #ifdef __i386__
 211 /* Long double constants are broken on i386. */
 212 #define LD80C(m, ex, v) {                                               \
 213         .xbits.man = __CONCAT(m, ULL),                                  \
 214         .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0),      \
 215 }
 216 #else
 217 /* The above works on non-i386 too, but we use this to check v. */
 218 #define LD80C(m, ex, v) { .e = (v), }
 219 #endif
 220
 221 #ifdef FLT_EVAL_METHOD
 222 /*
 223  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
 224  */
 225 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
 226 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
 227 #else
 228 #define STRICT_ASSIGN(type, lval, rval) do {    \
 229         volatile type __lval;                   \
 230                                                 \
 231         if (sizeof(type) >= sizeof(long double))        \
 232                 (lval) = (rval);                \
 233         else {                                  \
 234                 __lval = (rval);                \
 235                 (lval) = __lval;                \
 236         }                                       \
 237 } while (0)
 238 #endif
 239 #endif /* FLT_EVAL_METHOD */
 240
 241 /* Support switching the mode to FP_PE if necessary. */
 242 #if defined(__i386__) && !defined(NO_FPSETPREC)
 243 #define ENTERI()                                \
 244         long double __retval;                   \
 245         fp_prec_t __oprec;                      \
 246                                                 \
 247         if ((__oprec = fpgetprec()) != FP_PE)   \
 248                 fpsetprec(FP_PE)
 249 #define RETURNI(x) do {                         \
 250         __retval = (x);                         \
 251         if (__oprec != FP_PE)                   \
 252                 fpsetprec(__oprec);             \
 253         RETURNF(__retval);                      \
 254 } while (0)
 255 #else
 256 #define ENTERI(x)
 257 #define RETURNI(x)      RETURNF(x)
 258 #endif
 259
 260 /* Default return statement if hack*_t() is not used. */
 261 #define      RETURNF(v)      return (v)
 262
 263 /*
 264  * Common routine to process the arguments to nan(), nanf(), and nanl().
 265  */
 266 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
 267
 268 #ifdef _COMPLEX_H
 269
 270 /*
 271  * C99 specifies that complex numbers have the same representation as
 272  * an array of two elements, where the first element is the real part
 273  * and the second element is the imaginary part.
 274  */
 275 typedef union {
 276         float complex f;
 277         float a[2];
 278 } float_complex;
 279 typedef union {
 280         double complex f;
 281         double a[2];
 282 } double_complex;
 283 typedef union {
 284         long double complex f;
 285         long double a[2];
 286 } long_double_complex;
 287 #define REALPART(z)     ((z).a[0])
 288 #define IMAGPART(z)     ((z).a[1])
 289
 290 /*
 291  * Inline functions that can be used to construct complex values.
 292  *
 293  * The C99 standard intends x+I*y to be used for this, but x+I*y is
 294  * currently unusable in general since gcc introduces many overflow,
 295  * underflow, sign and efficiency bugs by rewriting I*y as
 296  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
 297  * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
 298  * to -0.0+I*0.0.
 299  */
 300 static __inline float complex
 301 cpackf(float x, float y)
 302 {
 303         float_complex z;
 304
 305         REALPART(z) = x;
 306         IMAGPART(z) = y;
 307         return (z.f);
 308 }
 309
 310 static __inline double complex
 311 cpack(double x, double y)
 312 {
 313         double_complex z;
 314
 315         REALPART(z) = x;
 316         IMAGPART(z) = y;
 317         return (z.f);
 318 }
 319
 320 static __inline long double complex
 321 cpackl(long double x, long double y)
 322 {
 323         long_double_complex z;
 324
 325         REALPART(z) = x;
 326         IMAGPART(z) = y;
 327         return (z.f);
 328 }
 329 #endif /* _COMPLEX_H */
 330
 331 #ifdef __GNUCLIKE_ASM
 332
 333 /* Asm versions of some functions. */
 334
 335 #ifdef __amd64__
 336 static __inline int
 337 irint(double x)
 338 {
 339         int n;
 340
 341         asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
 342         return (n);
 343 }
 344 #define HAVE_EFFICIENT_IRINT
 345 #endif
 346
 347 #ifdef __i386__
 348 static __inline int
 349 irint(double x)
 350 {
 351         int n;
 352
 353         asm("fistl %0" : "=m" (n) : "t" (x));
 354         return (n);
 355 }
 356 #define HAVE_EFFICIENT_IRINT
 357 #endif
 358
 359 #if defined(__amd64__) || defined(__i386__)
 360 static __inline int
 361 irintl(long double x)
 362 {
 363         int n;
 364
 365         asm("fistl %0" : "=m" (n) : "t" (x));
 366         return (n);
 367 }
 368 #define HAVE_EFFICIENT_IRINTL
 369 #endif
 370
 371 #endif /* __GNUCLIKE_ASM */
 372
 373 /*
 374  * ieee style elementary functions
 375  *
 376  * We rename functions here to improve other sources' diffability
 377  * against fdlibm.
 378  */
 379 #define __ieee754_sqrt  sqrt
 380 #define __ieee754_acos  acos
 381 #define __ieee754_acosh acosh
 382 #define __ieee754_log   log
 383 #define __ieee754_log2  log2
 384 #define __ieee754_atanh atanh
 385 #define __ieee754_asin  asin
 386 #define __ieee754_atan2 atan2
 387 #define __ieee754_exp   exp
 388 #define __ieee754_cosh  cosh
 389 #define __ieee754_fmod  fmod
 390 #define __ieee754_pow   pow
 391 #define __ieee754_lgamma lgamma
 392 #define __ieee754_gamma gamma
 393 #define __ieee754_lgamma_r lgamma_r
 394 #define __ieee754_gamma_r gamma_r
 395 #define __ieee754_log10 log10
 396 #define __ieee754_sinh  sinh
 397 #define __ieee754_hypot hypot
 398 #define __ieee754_j0    j0
 399 #define __ieee754_j1    j1
 400 #define __ieee754_y0    y0
 401 #define __ieee754_y1    y1
 402 #define __ieee754_jn    jn
 403 #define __ieee754_yn    yn
 404 #define __ieee754_remainder remainder
 405 #define __ieee754_scalb scalb
 406 #define __ieee754_sqrtf sqrtf
 407 #define __ieee754_acosf acosf
 408 #define __ieee754_acoshf acoshf
 409 #define __ieee754_logf  logf
 410 #define __ieee754_atanhf atanhf
 411 #define __ieee754_asinf asinf
 412 #define __ieee754_atan2f atan2f
 413 #define __ieee754_expf  expf
 414 #define __ieee754_coshf coshf
 415 #define __ieee754_fmodf fmodf
 416 #define __ieee754_powf  powf
 417 #define __ieee754_lgammaf lgammaf
 418 #define __ieee754_gammaf gammaf
 419 #define __ieee754_lgammaf_r lgammaf_r
 420 #define __ieee754_gammaf_r gammaf_r
 421 #define __ieee754_log10f log10f
 422 #define __ieee754_log2f log2f
 423 #define __ieee754_sinhf sinhf
 424 #define __ieee754_hypotf hypotf
 425 #define __ieee754_j0f   j0f
 426 #define __ieee754_j1f   j1f
 427 #define __ieee754_y0f   y0f
 428 #define __ieee754_y1f   y1f
 429 #define __ieee754_jnf   jnf
 430 #define __ieee754_ynf   ynf
 431 #define __ieee754_remainderf remainderf
 432 #define __ieee754_scalbf scalbf
 433
 434 /* fdlibm kernel function */
 435 int     __kernel_rem_pio2(double*,double*,int,int,int);
 436
 437 /* double precision kernel functions */
 438 #ifndef INLINE_REM_PIO2
 439 int     __ieee754_rem_pio2(double,double*);
 440 #endif
 441 double  __kernel_sin(double,double,int);
 442 double  __kernel_cos(double,double);
 443 double  __kernel_tan(double,double,int);
 444 double  __ldexp_exp(double,int);
 445 #ifdef _COMPLEX_H
 446 double complex __ldexp_cexp(double complex,int);
 447 #endif
 448
 449 /* float precision kernel functions */
 450 #ifndef INLINE_REM_PIO2F
 451 int     __ieee754_rem_pio2f(float,double*);
 452 #endif
 453 #ifndef INLINE_KERNEL_SINDF
 454 float   __kernel_sindf(double);
 455 #endif
 456 #ifndef INLINE_KERNEL_COSDF
 457 float   __kernel_cosdf(double);
 458 #endif
 459 #ifndef INLINE_KERNEL_TANDF
 460 float   __kernel_tandf(double,int);
 461 #endif
 462 float   __ldexp_expf(float,int);
 463 #ifdef _COMPLEX_H
 464 float complex __ldexp_cexpf(float complex,int);
 465 #endif
 466
 467 /* long double precision kernel functions */
 468 long double __kernel_sinl(long double, long double, int);
 469 long double __kernel_cosl(long double, long double);
 470 long double __kernel_tanl(long double, long double, int);
 471
 472 #endif /* !_MATH_PRIVATE_H_ */