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 0x7f7fffffL 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-(1L<<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 #ifdef _SCALB_INT 143 extern double scalb __P((double, int)); 144 #else 145 extern double scalb __P((double, double)); 146 #endif 147 extern double significand __P((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_log10 __P((double)); 164 extern double __ieee754_sinh __P((double)); 165 extern double __ieee754_hypot __P((double,double)); 166 extern double __ieee754_j0 __P((double)); 167 extern double __ieee754_j1 __P((double)); 168 extern double __ieee754_y0 __P((double)); 169 extern double __ieee754_y1 __P((double)); 170 extern double __ieee754_jn __P((int,double)); 171 extern double __ieee754_yn __P((int,double)); 172 extern double __ieee754_remainder __P((double,double)); 173 extern __int32_t __ieee754_rem_pio2 __P((double,double*)); 174 #ifdef _SCALB_INT 175 extern double __ieee754_scalb __P((double,int)); 176 #else 177 extern double __ieee754_scalb __P((double,double)); 178 #endif 179 180 /* fdlibm kernel function */ 181 extern double __kernel_standard __P((double,double,int)); 182 extern double __kernel_sin __P((double,double,int)); 183 extern double __kernel_cos __P((double,double)); 184 extern double __kernel_tan __P((double,double,int)); 185 extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*)); 186 187 /* Undocumented float functions. */ 188 #ifdef _SCALB_INT 189 extern float scalbf __P((float, int)); 190 #else 191 extern float scalbf __P((float, float)); 192 #endif 193 extern float significandf __P((float)); 194 195 /* ieee style elementary float functions */ 196 extern float __ieee754_sqrtf __P((float)); 197 extern float __ieee754_acosf __P((float)); 198 extern float __ieee754_acoshf __P((float)); 199 extern float __ieee754_logf __P((float)); 200 extern float __ieee754_atanhf __P((float)); 201 extern float __ieee754_asinf __P((float)); 202 extern float __ieee754_atan2f __P((float,float)); 203 extern float __ieee754_expf __P((float)); 204 extern float __ieee754_coshf __P((float)); 205 extern float __ieee754_fmodf __P((float,float)); 206 extern float __ieee754_powf __P((float,float)); 207 extern float __ieee754_lgammaf_r __P((float,int *)); 208 extern float __ieee754_gammaf_r __P((float,int *)); 209 extern float __ieee754_log10f __P((float)); 210 extern float __ieee754_sinhf __P((float)); 211 extern float __ieee754_hypotf __P((float,float)); 212 extern float __ieee754_j0f __P((float)); 213 extern float __ieee754_j1f __P((float)); 214 extern float __ieee754_y0f __P((float)); 215 extern float __ieee754_y1f __P((float)); 216 extern float __ieee754_jnf __P((int,float)); 217 extern float __ieee754_ynf __P((int,float)); 218 extern float __ieee754_remainderf __P((float,float)); 219 extern __int32_t __ieee754_rem_pio2f __P((float,float*)); 220 #ifdef _SCALB_INT 221 extern float __ieee754_scalbf __P((float,int)); 222 #else 223 extern float __ieee754_scalbf __P((float,float)); 224 #endif 225 226 /* float versions of fdlibm kernel functions */ 227 extern float __kernel_sinf __P((float,float,int)); 228 extern float __kernel_cosf __P((float,float)); 229 extern float __kernel_tanf __P((float,float,int)); 230 extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const __int32_t*)); 231 232 /* The original code used statements like 233 n0 = ((*(int*)&one)>>29)^1; * index of high word * 234 ix0 = *(n0+(int*)&x); * high word of x * 235 ix1 = *((1-n0)+(int*)&x); * low word of x * 236 to dig two 32 bit words out of the 64 bit IEEE floating point 237 value. That is non-ANSI, and, moreover, the gcc instruction 238 scheduler gets it wrong. We instead use the following macros. 239 Unlike the original code, we determine the endianness at compile 240 time, not at run time; I don't see much benefit to selecting 241 endianness at run time. */ 242 243 #ifndef __IEEE_BIG_ENDIAN 244 #ifndef __IEEE_LITTLE_ENDIAN 245 #error Must define endianness 246 #endif 247 #endif 248 249 /* A union which permits us to convert between a double and two 32 bit 250 ints. */ 251 252 #ifdef __IEEE_BIG_ENDIAN 253 254 typedef union 255 { 256 double value; 257 struct 258 { 259 __uint32_t msw; 260 __uint32_t lsw; 261 } parts; 262 } ieee_double_shape_type; 263 264 #endif 265 266 #ifdef __IEEE_LITTLE_ENDIAN 267 268 typedef union 269 { 270 double value; 271 struct 272 { 273 __uint32_t lsw; 274 __uint32_t msw; 275 } parts; 276 } ieee_double_shape_type; 277 278 #endif 279 280 /* Get two 32 bit ints from a double. */ 281 282 #define EXTRACT_WORDS(ix0,ix1,d) \ 283 do { \ 284 ieee_double_shape_type ew_u; \ 285 ew_u.value = (d); \ 286 (ix0) = ew_u.parts.msw; \ 287 (ix1) = ew_u.parts.lsw; \ 288 } while (0) 289 290 /* Get the more significant 32 bit int from a double. */ 291 292 #define GET_HIGH_WORD(i,d) \ 293 do { \ 294 ieee_double_shape_type gh_u; \ 295 gh_u.value = (d); \ 296 (i) = gh_u.parts.msw; \ 297 } while (0) 298 299 /* Get the less significant 32 bit int from a double. */ 300 301 #define GET_LOW_WORD(i,d) \ 302 do { \ 303 ieee_double_shape_type gl_u; \ 304 gl_u.value = (d); \ 305 (i) = gl_u.parts.lsw; \ 306 } while (0) 307 308 /* Set a double from two 32 bit ints. */ 309 310 #define INSERT_WORDS(d,ix0,ix1) \ 311 do { \ 312 ieee_double_shape_type iw_u; \ 313 iw_u.parts.msw = (ix0); \ 314 iw_u.parts.lsw = (ix1); \ 315 (d) = iw_u.value; \ 316 } while (0) 317 318 /* Set the more significant 32 bits of a double from an int. */ 319 320 #define SET_HIGH_WORD(d,v) \ 321 do { \ 322 ieee_double_shape_type sh_u; \ 323 sh_u.value = (d); \ 324 sh_u.parts.msw = (v); \ 325 (d) = sh_u.value; \ 326 } while (0) 327 328 /* Set the less significant 32 bits of a double from an int. */ 329 330 #define SET_LOW_WORD(d,v) \ 331 do { \ 332 ieee_double_shape_type sl_u; \ 333 sl_u.value = (d); \ 334 sl_u.parts.lsw = (v); \ 335 (d) = sl_u.value; \ 336 } while (0) 337 338 /* A union which permits us to convert between a float and a 32 bit 339 int. */ 340 341 typedef union 342 { 343 float value; 344 __uint32_t word; 345 } ieee_float_shape_type; 346 347 /* Get a 32 bit int from a float. */ 348 349 #define GET_FLOAT_WORD(i,d) \ 350 do { \ 351 ieee_float_shape_type gf_u; \ 352 gf_u.value = (d); \ 353 (i) = gf_u.word; \ 354 } while (0) 355 356 /* Set a float from a 32 bit int. */ 357 358 #define SET_FLOAT_WORD(d,i) \ 359 do { \ 360 ieee_float_shape_type sf_u; \ 361 sf_u.word = (i); \ 362 (d) = sf_u.value; \ 363 } while (0) 364 365 /* Macros to avoid undefined behaviour that can arise if the amount 366 of a shift is exactly equal to the size of the shifted operand. */ 367 368 #define SAFE_LEFT_SHIFT(op,amt) \ 369 (((amt) < 8 * sizeof(op)) ? ((op) << (amt)) : 0) 370 371 #define SAFE_RIGHT_SHIFT(op,amt) \ 372 (((amt) < 8 * sizeof(op)) ? ((op) >> (amt)) : 0) 373 374 #ifdef _COMPLEX_H 375 376 /* 377 * Quoting from ISO/IEC 9899:TC2: 378 * 379 * 6.2.5.13 Types 380 * Each complex type has the same representation and alignment requirements as 381 * an array type containing exactly two elements of the corresponding real type; 382 * the first element is equal to the real part, and the second element to the 383 * imaginary part, of the complex number. 384 */ 385 typedef union { 386 float complex z; 387 float parts[2]; 388 } float_complex; 389 390 typedef union { 391 double complex z; 392 double parts[2]; 393 } double_complex; 394 395 typedef union { 396 long double complex z; 397 long double parts[2]; 398 } long_double_complex; 399 400 #define REAL_PART(z) ((z).parts[0]) 401 #define IMAG_PART(z) ((z).parts[1]) 402 403 #endif /* _COMPLEX_H */ 404 405