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