root/lib/math/fdlibm.h

/* [<][>][^][v][top][bottom][index][help] */

INCLUDED FROM


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

/* [<][>][^][v][top][bottom][index][help] */