toolchain/gcc/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
  19 /* REDHAT LOCAL: Default to XOPEN_MODE.  */
  20 #define _XOPEN_MODE
  21
  22 /* Most routines need to check whether a float is finite, infinite, or not a
  23    number, and many need to know whether the result of an operation will
  24    overflow.  These conditions depend on whether the largest exponent is
  25    used for NaNs & infinities, or whether it's used for finite numbers.  The
  26    macros below wrap up that kind of information:
  27
  28    FLT_UWORD_IS_FINITE(X)
  29         True if a positive float with bitmask X is finite.
  30
  31    FLT_UWORD_IS_NAN(X)
  32         True if a positive float with bitmask X is not a number.
  33
  34    FLT_UWORD_IS_INFINITE(X)
  35         True if a positive float with bitmask X is +infinity.
  36
  37    FLT_UWORD_MAX
  38         The bitmask of FLT_MAX.
  39
  40    FLT_UWORD_HALF_MAX
  41         The bitmask of FLT_MAX/2.
  42
  43    FLT_UWORD_EXP_MAX
  44         The bitmask of the largest finite exponent (129 if the largest
  45         exponent is used for finite numbers, 128 otherwise).
  46
  47    FLT_UWORD_LOG_MAX
  48         The bitmask of log(FLT_MAX), rounded down.  This value is the largest
  49         input that can be passed to exp() without producing overflow.
  50
  51    FLT_UWORD_LOG_2MAX
  52         The bitmask of log(2*FLT_MAX), rounded down.  This value is the
  53         largest input than can be passed to cosh() without producing
  54         overflow.
  55
  56    FLT_LARGEST_EXP
  57         The largest biased exponent that can be used for finite numbers
  58         (255 if the largest exponent is used for finite numbers, 254
  59         otherwise) */
  60
  61 #ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL
  62 #define FLT_UWORD_IS_FINITE(x) 1
  63 #define FLT_UWORD_IS_NAN(x) 0
  64 #define FLT_UWORD_IS_INFINITE(x) 0
  65 #define FLT_UWORD_MAX 0x7fffffff
  66 #define FLT_UWORD_EXP_MAX 0x43010000
  67 #define FLT_UWORD_LOG_MAX 0x42b2d4fc
  68 #define FLT_UWORD_LOG_2MAX 0x42b437e0
  69 #define HUGE ((float)0X1.FFFFFEP128)
  70 #else
  71 #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L)
  72 #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L)
  73 #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L)
  74 #define FLT_UWORD_MAX 0x7f7fffff
  75 #define FLT_UWORD_EXP_MAX 0x43000000
  76 #define FLT_UWORD_LOG_MAX 0x42b17217
  77 #define FLT_UWORD_LOG_2MAX 0x42b2d4fc
  78 #define HUGE ((float)3.40282346638528860e+38)
  79 #endif
  80 #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1<<23))
  81 #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23)
  82
  83 /* Many routines check for zero and subnormal numbers.  Such things depend
  84    on whether the target supports denormals or not:
  85
  86    FLT_UWORD_IS_ZERO(X)
  87         True if a positive float with bitmask X is +0.  Without denormals,
  88         any float with a zero exponent is a +0 representation.  With
  89         denormals, the only +0 representation is a 0 bitmask.
  90
  91    FLT_UWORD_IS_SUBNORMAL(X)
  92         True if a non-zero positive float with bitmask X is subnormal.
  93         (Routines should check for zeros first.)
  94
  95    FLT_UWORD_MIN
  96         The bitmask of the smallest float above +0.  Call this number
  97         REAL_FLT_MIN...
  98
  99    FLT_UWORD_EXP_MIN
 100         The bitmask of the float representation of REAL_FLT_MIN's exponent.
 101
 102    FLT_UWORD_LOG_MIN
 103         The bitmask of |log(REAL_FLT_MIN)|, rounding down.
 104
 105    FLT_SMALLEST_EXP
 106         REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported,
 107         -22 if they are).
 108 */
 109
 110 #ifdef _FLT_NO_DENORMALS
 111 #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L)
 112 #define FLT_UWORD_IS_SUBNORMAL(x) 0
 113 #define FLT_UWORD_MIN 0x00800000
 114 #define FLT_UWORD_EXP_MIN 0x42fc0000
 115 #define FLT_UWORD_LOG_MIN 0x42aeac50
 116 #define FLT_SMALLEST_EXP 1
 117 #else
 118 #define FLT_UWORD_IS_ZERO(x) ((x)==0)
 119 #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L)
 120 #define FLT_UWORD_MIN 0x00000001
 121 #define FLT_UWORD_EXP_MIN 0x43160000
 122 #define FLT_UWORD_LOG_MIN 0x42cff1b5
 123 #define FLT_SMALLEST_EXP -22
 124 #endif
 125
 126 #ifdef __STDC__
 127 #undef __P
 128 #define __P(p)  p
 129 #else
 130 #define __P(p)  ()
 131 #endif
 132
 133 /*
 134  * set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
 135  * (one may replace the following line by "#include <values.h>")
 136  */
 137
 138 #define X_TLOSS         1.41484755040568800000e+16
 139
 140 /* Functions that are not documented, and are not in <math.h>.  */
 141
 142 extern double logb __P((double));
 143 #ifdef _SCALB_INT
 144 extern double scalb __P((double, int));
 145 #else
 146 extern double scalb __P((double, double));
 147 #endif
 148 extern double significand __P((double));
 149
 150 /* ieee style elementary functions */
 151 extern double __ieee754_sqrt __P((double));
 152 extern double __ieee754_acos __P((double));
 153 extern double __ieee754_acosh __P((double));
 154 extern double __ieee754_log __P((double));
 155 extern double __ieee754_atanh __P((double));
 156 extern double __ieee754_asin __P((double));
 157 extern double __ieee754_atan2 __P((double,double));
 158 extern double __ieee754_exp __P((double));
 159 extern double __ieee754_cosh __P((double));
 160 extern double __ieee754_fmod __P((double,double));
 161 extern double __ieee754_pow __P((double,double));
 162 extern double __ieee754_lgamma_r __P((double,int *));
 163 extern double __ieee754_gamma_r __P((double,int *));
 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 extern float logbf __P((float));
 190 #ifdef _SCALB_INT
 191 extern float scalbf __P((float, int));
 192 #else
 193 extern float scalbf __P((float, float));
 194 #endif
 195 extern float significandf __P((float));
 196
 197 /* ieee style elementary float functions */
 198 extern float __ieee754_sqrtf __P((float));
 199 extern float __ieee754_acosf __P((float));
 200 extern float __ieee754_acoshf __P((float));
 201 extern float __ieee754_logf __P((float));
 202 extern float __ieee754_atanhf __P((float));
 203 extern float __ieee754_asinf __P((float));
 204 extern float __ieee754_atan2f __P((float,float));
 205 extern float __ieee754_expf __P((float));
 206 extern float __ieee754_coshf __P((float));
 207 extern float __ieee754_fmodf __P((float,float));
 208 extern float __ieee754_powf __P((float,float));
 209 extern float __ieee754_lgammaf_r __P((float,int *));
 210 extern float __ieee754_gammaf_r __P((float,int *));
 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 /* float versions of fdlibm kernel functions */
 229 extern float __kernel_sinf __P((float,float,int));
 230 extern float __kernel_cosf __P((float,float));
 231 extern float __kernel_tanf __P((float,float,int));
 232 extern int   __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*));
 233
 234 /* The original code used statements like
 235         n0 = ((*(int*)&one)>>29)^1;             * index of high word *
 236         ix0 = *(n0+(int*)&x);                   * high word of x *
 237         ix1 = *((1-n0)+(int*)&x);               * low word of x *
 238    to dig two 32 bit words out of the 64 bit IEEE floating point
 239    value.  That is non-ANSI, and, moreover, the gcc instruction
 240    scheduler gets it wrong.  We instead use the following macros.
 241    Unlike the original code, we determine the endianness at compile
 242    time, not at run time; I don't see much benefit to selecting
 243    endianness at run time.  */
 244
 245 #ifndef __IEEE_BIG_ENDIAN
 246 #ifndef __IEEE_LITTLE_ENDIAN
 247  #error Must define endianness
 248 #endif
 249 #endif
 250
 251 /* A union which permits us to convert between a double and two 32 bit
 252    ints.  */
 253
 254 #ifdef __IEEE_BIG_ENDIAN
 255
 256 typedef union
 257 {
 258   double value;
 259   struct
 260   {
 261     __uint32_t msw;
 262     __uint32_t lsw;
 263   } parts;
 264 } ieee_double_shape_type;
 265
 266 #endif
 267
 268 #ifdef __IEEE_LITTLE_ENDIAN
 269
 270 typedef union
 271 {
 272   double value;
 273   struct
 274   {
 275     __uint32_t lsw;
 276     __uint32_t msw;
 277   } parts;
 278 } ieee_double_shape_type;
 279
 280 #endif
 281
 282 /* Get two 32 bit ints from a double.  */
 283
 284 #define EXTRACT_WORDS(ix0,ix1,d)                                \
 285 do {                                                            \
 286   ieee_double_shape_type ew_u;                                  \
 287   ew_u.value = (d);                                             \
 288   (ix0) = ew_u.parts.msw;                                       \
 289   (ix1) = ew_u.parts.lsw;                                       \
 290 } while (0)
 291
 292 /* Get the more significant 32 bit int from a double.  */
 293
 294 #define GET_HIGH_WORD(i,d)                                      \
 295 do {                                                            \
 296   ieee_double_shape_type gh_u;                                  \
 297   gh_u.value = (d);                                             \
 298   (i) = gh_u.parts.msw;                                         \
 299 } while (0)
 300
 301 /* Get the less significant 32 bit int from a double.  */
 302
 303 #define GET_LOW_WORD(i,d)                                       \
 304 do {                                                            \
 305   ieee_double_shape_type gl_u;                                  \
 306   gl_u.value = (d);                                             \
 307   (i) = gl_u.parts.lsw;                                         \
 308 } while (0)
 309
 310 /* Set a double from two 32 bit ints.  */
 311
 312 #define INSERT_WORDS(d,ix0,ix1)                                 \
 313 do {                                                            \
 314   ieee_double_shape_type iw_u;                                  \
 315   iw_u.parts.msw = (ix0);                                       \
 316   iw_u.parts.lsw = (ix1);                                       \
 317   (d) = iw_u.value;                                             \
 318 } while (0)
 319
 320 /* Set the more significant 32 bits of a double from an int.  */
 321
 322 #define SET_HIGH_WORD(d,v)                                      \
 323 do {                                                            \
 324   ieee_double_shape_type sh_u;                                  \
 325   sh_u.value = (d);                                             \
 326   sh_u.parts.msw = (v);                                         \
 327   (d) = sh_u.value;                                             \
 328 } while (0)
 329
 330 /* Set the less significant 32 bits of a double from an int.  */
 331
 332 #define SET_LOW_WORD(d,v)                                       \
 333 do {                                                            \
 334   ieee_double_shape_type sl_u;                                  \
 335   sl_u.value = (d);                                             \
 336   sl_u.parts.lsw = (v);                                         \
 337   (d) = sl_u.value;                                             \
 338 } while (0)
 339
 340 /* A union which permits us to convert between a float and a 32 bit
 341    int.  */
 342
 343 typedef union
 344 {
 345   float value;
 346   __uint32_t word;
 347 } ieee_float_shape_type;
 348
 349 /* Get a 32 bit int from a float.  */
 350
 351 #define GET_FLOAT_WORD(i,d)                                     \
 352 do {                                                            \
 353   ieee_float_shape_type gf_u;                                   \
 354   gf_u.value = (d);                                             \
 355   (i) = gf_u.word;                                              \
 356 } while (0)
 357
 358 /* Set a float from a 32 bit int.  */
 359
 360 #define SET_FLOAT_WORD(d,i)                                     \
 361 do {                                                            \
 362   ieee_float_shape_type sf_u;                                   \
 363   sf_u.word = (i);                                              \
 364   (d) = sf_u.value;                                             \
 365 } while (0)