newlib/libm/common/fdlibm.h

   1
   2 /* @(#)fdlibm.h 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 /* REDHAT LOCAL: Include files.  */
  15 #include <math.h>
  16 #include <sys/types.h>
  17 #include <machine/ieeefp.h>
  18 #include "math_config.h"
  19
  20 /* Most routines need to check whether a float is finite, infinite, or not a
  21    number, and many need to know whether the result of an operation will
  22    overflow.  These conditions depend on whether the largest exponent is
  23    used for NaNs & infinities, or whether it's used for finite numbers.  The
  24    macros below wrap up that kind of information:
  25
  26    FLT_UWORD_IS_FINITE(X)
  27         True if a positive float with bitmask X is finite.
  28
  29    FLT_UWORD_IS_NAN(X)
  30         True if a positive float with bitmask X is not a number.
  31
  32    FLT_UWORD_IS_INFINITE(X)
  33         True if a positive float with bitmask X is +infinity.
  34
  35    FLT_UWORD_MAX
  36         The bitmask of FLT_MAX.
  37
  38    FLT_UWORD_HALF_MAX
  39         The bitmask of FLT_MAX/2.
  40
  41    FLT_UWORD_EXP_MAX
  42         The bitmask of the largest finite exponent (129 if the largest
  43         exponent is used for finite numbers, 128 otherwise).
  44
  45    FLT_UWORD_LOG_MAX
  46         The bitmask of log(FLT_MAX), rounded down.  This value is the largest
  47         input that can be passed to exp() without producing overflow.
  48
  49    FLT_UWORD_LOG_2MAX
  50         The bitmask of log(2*FLT_MAX), rounded down.  This value is the
  51         largest input than can be passed to cosh() without producing
  52         overflow.
  53
  54    FLT_LARGEST_EXP
  55         The largest biased exponent that can be used for finite numbers
  56         (255 if the largest exponent is used for finite numbers, 254
  57         otherwise) */
  58
  59 #ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
  60 #define FLT_UWORD_IS_FINITE(x) 1
  61 #define FLT_UWORD_IS_NAN(x) 0
  62 #define FLT_UWORD_IS_INFINITE(x) 0
  63 #define FLT_UWORD_MAX 0x7fffffff
  64 #define FLT_UWORD_EXP_MAX 0x43010000
  65 #define FLT_UWORD_LOG_MAX 0x42b2d4fc
  66 #define FLT_UWORD_LOG_2MAX 0x42b437e0
  67 #define HUGE ((float)0X1.FFFFFEP128)
  68 #else
  69 #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
  70 #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
  71 #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
  72 #define FLT_UWORD_MAX 0x7f7fffffL
  73 #define FLT_UWORD_EXP_MAX 0x43000000
  74 #define FLT_UWORD_LOG_MAX 0x42b17217
  75 #define FLT_UWORD_LOG_2MAX 0x42b2d4fc
  76 #define HUGE ((float)3.40282346638528860e+38)
  77 #endif
  78 #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23))
  79 #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
  80
  81 /* Many routines check for zero and subnormal numbers.  Such things depend
  82    on whether the target supports denormals or not:
  83
  84    FLT_UWORD_IS_ZERO(X)
  85         True if a positive float with bitmask X is +0.  Without denormals,
  86         any float with a zero exponent is a +0 representation.  With
  87         denormals, the only +0 representation is a 0 bitmask.
  88
  89    FLT_UWORD_IS_SUBNORMAL(X)
  90         True if a non-zero positive float with bitmask X is subnormal.
  91         (Routines should check for zeros first.)
  92
  93    FLT_UWORD_MIN
  94         The bitmask of the smallest float above +0.  Call this number
  95         REAL_FLT_MIN...
  96
  97    FLT_UWORD_EXP_MIN
  98         The bitmask of the float representation of REAL_FLT_MIN's exponent.
  99
 100    FLT_UWORD_LOG_MIN
 101         The bitmask of |log(REAL_FLT_MIN)|, rounding down.
 102
 103    FLT_SMALLEST_EXP
 104         REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
 105         -22 if they are).
 106 */
 107
 108 #ifdef _FLT_NO_DENORMALS
 109 #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
 110 #define FLT_UWORD_IS_SUBNORMAL(x) 0
 111 #define FLT_UWORD_MIN 0x00800000
 112 #define FLT_UWORD_EXP_MIN 0x42fc0000
 113 #define FLT_UWORD_LOG_MIN 0x42aeac50
 114 #define FLT_SMALLEST_EXP 1
 115 #else
 116 #define FLT_UWORD_IS_ZERO(x) ((x)==0)
 117 #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
 118 #define FLT_UWORD_MIN 0x00000001
 119 #define FLT_UWORD_EXP_MIN 0x43160000
 120 #define FLT_UWORD_LOG_MIN 0x42cff1b5
 121 #define FLT_SMALLEST_EXP -22
 122 #endif
 123
 124 #ifdef __STDC__
 125 #undef __P
 126 #define __P(p)  p
 127 #else
 128 #define __P(p)  ()
 129 #endif
 130
 131 /*
 132  * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
 133  * (one may replace the following line by "#include <values.h>")
 134  */
 135
 136 #define X_TLOSS         1.41484755040568800000e+16
 137
 138 /* Functions that are not documented, and are not in <math.h>.  */
 139
 140 #ifdef _SCALB_INT
 141 extern double scalb __P((double, int));
 142 #else
 143 extern double scalb __P((double, double));
 144 #endif
 145 extern double significand __P((double));
 146
 147 extern long double __ieee754_hypotl __P((long double, long double));
 148
 149 /* ieee style elementary functions */
 150 extern double __ieee754_sqrt __P((double));
 151 extern double __ieee754_acos __P((double));
 152 extern double __ieee754_acosh __P((double));
 153 extern double __ieee754_log __P((double));
 154 extern double __ieee754_atanh __P((double));
 155 extern double __ieee754_asin __P((double));
 156 extern double __ieee754_atan2 __P((double,double));
 157 extern double __ieee754_exp __P((double));
 158 extern double __ieee754_cosh __P((double));
 159 extern double __ieee754_fmod __P((double,double));
 160 extern double __ieee754_pow __P((double,double));
 161 extern double __ieee754_lgamma_r __P((double,int *));
 162 extern double __ieee754_gamma_r __P((double,int *));
 163 extern double __ieee754_tgamma __P((double));
 164 extern double __ieee754_log10 __P((double));
 165 extern double __ieee754_sinh __P((double));
 166 extern double __ieee754_hypot __P((double,double));
 167 extern double __ieee754_j0 __P((double));
 168 extern double __ieee754_j1 __P((double));
 169 extern double __ieee754_y0 __P((double));
 170 extern double __ieee754_y1 __P((double));
 171 extern double __ieee754_jn __P((int,double));
 172 extern double __ieee754_yn __P((int,double));
 173 extern double __ieee754_remainder __P((double,double));
 174 extern __int32_t __ieee754_rem_pio2 __P((double,double*));
 175 #ifdef _SCALB_INT
 176 extern double __ieee754_scalb __P((double,int));
 177 #else
 178 extern double __ieee754_scalb __P((double,double));
 179 #endif
 180
 181 /* fdlibm kernel function */
 182 extern double __kernel_standard __P((double,double,int));
 183 extern double __kernel_sin __P((double,double,int));
 184 extern double __kernel_cos __P((double,double));
 185 extern double __kernel_tan __P((double,double,int));
 186 extern int    __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*));
 187
 188 /* Undocumented float functions.  */
 189 #ifdef _SCALB_INT
 190 extern float scalbf __P((float, int));
 191 #else
 192 extern float scalbf __P((float, float));
 193 #endif
 194 extern float significandf __P((float));
 195
 196 /* ieee style elementary float functions */
 197 extern float __ieee754_sqrtf __P((float));
 198 extern float __ieee754_acosf __P((float));
 199 extern float __ieee754_acoshf __P((float));
 200 extern float __ieee754_logf __P((float));
 201 extern float __ieee754_atanhf __P((float));
 202 extern float __ieee754_asinf __P((float));
 203 extern float __ieee754_atan2f __P((float,float));
 204 extern float __ieee754_expf __P((float));
 205 extern float __ieee754_coshf __P((float));
 206 extern float __ieee754_fmodf __P((float,float));
 207 extern float __ieee754_powf __P((float,float));
 208 extern float __ieee754_lgammaf_r __P((float,int *));
 209 extern float __ieee754_gammaf_r __P((float,int *));
 210 extern float __ieee754_tgammaf __P((float));
 211 extern float __ieee754_log10f __P((float));
 212 extern float __ieee754_sinhf __P((float));
 213 extern float __ieee754_hypotf __P((float,float));
 214 extern float __ieee754_j0f __P((float));
 215 extern float __ieee754_j1f __P((float));
 216 extern float __ieee754_y0f __P((float));
 217 extern float __ieee754_y1f __P((float));
 218 extern float __ieee754_jnf __P((int,float));
 219 extern float __ieee754_ynf __P((int,float));
 220 extern float __ieee754_remainderf __P((float,float));
 221 extern __int32_t __ieee754_rem_pio2f __P((float,float*));
 222 #ifdef _SCALB_INT
 223 extern float __ieee754_scalbf __P((float,int));
 224 #else
 225 extern float __ieee754_scalbf __P((float,float));
 226 #endif
 227
 228 #if !__OBSOLETE_MATH
 229 /* The new math code does not provide separate wrapper function
 230    for error handling, so the extern symbol is called directly.
 231    This is valid as long as there are no namespace issues (the
 232    extern symbol is reserved whenever the caller is reserved)
 233    and there are no observable error handling side effects.  */
 234 # define __ieee754_exp(x) exp(x)
 235 # define __ieee754_log(x) log(x)
 236 # define __ieee754_pow(x,y) pow(x,y)
 237 # define __ieee754_expf(x) expf(x)
 238 # define __ieee754_logf(x) logf(x)
 239 # define __ieee754_powf(x,y) powf(x,y)
 240 #endif
 241
 242 /* float versions of fdlibm kernel functions */
 243 extern float __kernel_sinf __P((float,float,int));
 244 extern float __kernel_cosf __P((float,float));
 245 extern float __kernel_tanf __P((float,float,int));
 246 extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
 247
 248 /* The original code used statements like
 249         n0 = ((*(int*)&one)>>29)^1;             * index of high word *
 250         ix0 = *(n0+(int*)&x);                   * high word of x *
 251         ix1 = *((1-n0)+(int*)&x);               * low word of x *
 252    to dig two 32 bit words out of the 64 bit IEEE floating point
 253    value.  That is non-ANSI, and, moreover, the gcc instruction
 254    scheduler gets it wrong.  We instead use the following macros.
 255    Unlike the original code, we determine the endianness at compile
 256    time, not at run time; I don't see much benefit to selecting
 257    endianness at run time.  */
 258
 259 #ifndef __IEEE_BIG_ENDIAN
 260 #ifndef __IEEE_LITTLE_ENDIAN
 261  #error Must define endianness
 262 #endif
 263 #endif
 264
 265 /* A union which permits us to convert between a double and two 32 bit
 266    ints.  */
 267
 268 #ifdef __IEEE_BIG_ENDIAN
 269
 270 typedef union
 271 {
 272   double value;
 273   struct
 274   {
 275     __uint32_t msw;
 276     __uint32_t lsw;
 277   } parts;
 278 } ieee_double_shape_type;
 279
 280 #endif
 281
 282 #ifdef __IEEE_LITTLE_ENDIAN
 283
 284 typedef union
 285 {
 286   double value;
 287   struct
 288   {
 289     __uint32_t lsw;
 290     __uint32_t msw;
 291   } parts;
 292 } ieee_double_shape_type;
 293
 294 #endif
 295
 296 /* Get two 32 bit ints from a double.  */
 297
 298 #define EXTRACT_WORDS(ix0,ix1,d)                                \
 299 do {                                                            \
 300   ieee_double_shape_type ew_u;                                  \
 301   ew_u.value = (d);                                             \
 302   (ix0) = ew_u.parts.msw;                                       \
 303   (ix1) = ew_u.parts.lsw;                                       \
 304 } while (0)
 305
 306 /* Get the more significant 32 bit int from a double.  */
 307
 308 #define GET_HIGH_WORD(i,d)                                      \
 309 do {                                                            \
 310   ieee_double_shape_type gh_u;                                  \
 311   gh_u.value = (d);                                             \
 312   (i) = gh_u.parts.msw;                                         \
 313 } while (0)
 314
 315 /* Get the less significant 32 bit int from a double.  */
 316
 317 #define GET_LOW_WORD(i,d)                                       \
 318 do {                                                            \
 319   ieee_double_shape_type gl_u;                                  \
 320   gl_u.value = (d);                                             \
 321   (i) = gl_u.parts.lsw;                                         \
 322 } while (0)
 323
 324 /* Set a double from two 32 bit ints.  */
 325
 326 #define INSERT_WORDS(d,ix0,ix1)                                 \
 327 do {                                                            \
 328   ieee_double_shape_type iw_u;                                  \
 329   iw_u.parts.msw = (ix0);                                       \
 330   iw_u.parts.lsw = (ix1);                                       \
 331   (d) = iw_u.value;                                             \
 332 } while (0)
 333
 334 /* Set the more significant 32 bits of a double from an int.  */
 335
 336 #define SET_HIGH_WORD(d,v)                                      \
 337 do {                                                            \
 338   ieee_double_shape_type sh_u;                                  \
 339   sh_u.value = (d);                                             \
 340   sh_u.parts.msw = (v);                                         \
 341   (d) = sh_u.value;                                             \
 342 } while (0)
 343
 344 /* Set the less significant 32 bits of a double from an int.  */
 345
 346 #define SET_LOW_WORD(d,v)                                       \
 347 do {                                                            \
 348   ieee_double_shape_type sl_u;                                  \
 349   sl_u.value = (d);                                             \
 350   sl_u.parts.lsw = (v);                                         \
 351   (d) = sl_u.value;                                             \
 352 } while (0)
 353
 354 /* A union which permits us to convert between a float and a 32 bit
 355    int.  */
 356
 357 typedef union
 358 {
 359   float value;
 360   __uint32_t word;
 361 } ieee_float_shape_type;
 362
 363 /* Get a 32 bit int from a float.  */
 364
 365 #define GET_FLOAT_WORD(i,d)                                     \
 366 do {                                                            \
 367   ieee_float_shape_type gf_u;                                   \
 368   gf_u.value = (d);                                             \
 369   (i) = gf_u.word;                                              \
 370 } while (0)
 371
 372 /* Set a float from a 32 bit int.  */
 373
 374 #define SET_FLOAT_WORD(d,i)                                     \
 375 do {                                                            \
 376   ieee_float_shape_type sf_u;                                   \
 377   sf_u.word = (i);                                              \
 378   (d) = sf_u.value;                                             \
 379 } while (0)
 380
 381 /* Macros to avoid undefined behaviour that can arise if the amount
 382    of a shift is exactly equal to the size of the shifted operand.  */
 383
 384 #define SAFE_LEFT_SHIFT(op,amt)                                 \
 385   (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0)
 386
 387 #define SAFE_RIGHT_SHIFT(op,amt)                                \
 388   (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0)
 389
 390 #ifdef  _COMPLEX_H
 391
 392 /*
 393  * Quoting from ISO/IEC 9899:TC2:
 394  *
 395  * 6.2.5.13 Types
 396  * Each complex type has the same representation and alignment requirements as
 397  * an array type containing exactly two elements of the corresponding real type;
 398  * the first element is equal to the real part, and the second element to the
 399  * imaginary part, of the complex number.
 400  */
 401 typedef union {
 402         float complex z;
 403         float parts[2];
 404 } float_complex;
 405
 406 typedef union {
 407         double complex z;
 408         double parts[2];
 409 } double_complex;
 410
 411 typedef union {
 412         long double complex z;
 413         long double parts[2];
 414 } long_double_complex;
 415
 416 #define REAL_PART(z)    ((z).parts[0])
 417 #define IMAG_PART(z)    ((z).parts[1])
 418
 419 #endif  /* _COMPLEX_H */
 420