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