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 #ifndef _FDLIBM_H_ 15 #define _FDLIBM_H_ 16 17 /* REDHAT LOCAL: Include files. */ 18 #ifndef _DEFAULT_SOURCE 19 #define _DEFAULT_SOURCE 20 #endif 21 #include <math.h> 22 #include <sys/types.h> 23 #include <machine/ieeefp.h> 24 #include <fenv.h> 25 #include "math_config.h" 26 27 /* Most routines need to check whether a float is finite, infinite, or not a 28 number, and many need to know whether the result of an operation will 29 overflow. These conditions depend on whether the largest exponent is 30 used for NaNs & infinities, or whether it's used for finite numbers. The 31 macros below wrap up that kind of information: 32 33 FLT_UWORD_IS_FINITE(X) 34 True if a positive float with bitmask X is finite. 35 36 FLT_UWORD_IS_NAN(X) 37 True if a positive float with bitmask X is not a number. 38 39 FLT_UWORD_IS_INFINITE(X) 40 True if a positive float with bitmask X is +infinity. 41 42 FLT_UWORD_MAX 43 The bitmask of FLT_MAX. 44 45 FLT_UWORD_HALF_MAX 46 The bitmask of FLT_MAX/2. 47 48 FLT_UWORD_EXP_MAX 49 The bitmask of the largest finite exponent (129 if the largest 50 exponent is used for finite numbers, 128 otherwise). 51 52 FLT_UWORD_LOG_MAX 53 The bitmask of log(FLT_MAX), rounded down. This value is the largest 54 input that can be passed to exp() without producing overflow. 55 56 FLT_UWORD_LOG_2MAX 57 The bitmask of log(2*FLT_MAX), rounded down. This value is the 58 largest input than can be passed to cosh() without producing 59 overflow. 60 61 FLT_LARGEST_EXP 62 The largest biased exponent that can be used for finite numbers 63 (255 if the largest exponent is used for finite numbers, 254 64 otherwise) */ 65 66 #ifdef _FLT_LARGEST_EXPONENT_IS_NORMAL 67 #define FLT_UWORD_IS_FINITE(x) 1 68 #define FLT_UWORD_IS_NAN(x) 0 69 #define FLT_UWORD_IS_INFINITE(x) 0 70 #define FLT_UWORD_MAX 0x7fffffff 71 #define FLT_UWORD_EXP_MAX 0x43010000 72 #define FLT_UWORD_LOG_MAX 0x42b2d4fc 73 #define FLT_UWORD_LOG_2MAX 0x42b437e0 74 #define HUGE ((float)0X1.FFFFFEP128) 75 #else 76 #define FLT_UWORD_IS_FINITE(x) ((x)<0x7f800000L) 77 #define FLT_UWORD_IS_NAN(x) ((x)>0x7f800000L) 78 #define FLT_UWORD_IS_INFINITE(x) ((x)==0x7f800000L) 79 #define FLT_UWORD_MAX 0x7f7fffffL 80 #define FLT_UWORD_EXP_MAX 0x43000000 81 #define FLT_UWORD_LOG_MAX 0x42b17217 82 #define FLT_UWORD_LOG_2MAX 0x42b2d4fc 83 #define HUGE ((float)3.40282346638528860e+38) 84 #endif 85 #define FLT_UWORD_HALF_MAX (FLT_UWORD_MAX-(1L<<23)) 86 #define FLT_LARGEST_EXP (FLT_UWORD_MAX>>23) 87 88 /* rounding mode tests; nearest if not set. Assumes hardware 89 * without rounding mode support uses nearest 90 */ 91 92 /* If there are rounding modes other than FE_TONEAREST defined, then 93 * add code to check which is active 94 */ 95 #if (defined(FE_UPWARD) + defined(FE_DOWNWARD) + defined(FE_TOWARDZERO)) >= 1 96 #define FE_DECL_ROUND(v) int v = fegetround() 97 #define __is_nearest(r) ((r) == FE_TONEAREST) 98 #else 99 #define FE_DECL_ROUND(v) 100 #define __is_nearest(r) 1 101 #endif 102 103 #ifdef FE_UPWARD 104 #define __is_upward(r) ((r) == FE_UPWARD) 105 #else 106 #define __is_upward(r) 0 107 #endif 108 109 #ifdef FE_DOWNWARD 110 #define __is_downward(r) ((r) == FE_DOWNWARD) 111 #else 112 #define __is_downward(r) 0 113 #endif 114 115 #ifdef FE_TOWARDZERO 116 #define __is_towardzero(r) ((r) == FE_TOWARDZERO) 117 #else 118 #define __is_towardzero(r) 0 119 #endif 120 121 /* Many routines check for zero and subnormal numbers. Such things depend 122 on whether the target supports denormals or not: 123 124 FLT_UWORD_IS_ZERO(X) 125 True if a positive float with bitmask X is +0. Without denormals, 126 any float with a zero exponent is a +0 representation. With 127 denormals, the only +0 representation is a 0 bitmask. 128 129 FLT_UWORD_IS_SUBNORMAL(X) 130 True if a non-zero positive float with bitmask X is subnormal. 131 (Routines should check for zeros first.) 132 133 FLT_UWORD_MIN 134 The bitmask of the smallest float above +0. Call this number 135 REAL_FLT_MIN... 136 137 FLT_UWORD_EXP_MIN 138 The bitmask of the float representation of REAL_FLT_MIN's exponent. 139 140 FLT_UWORD_LOG_MIN 141 The bitmask of |log(REAL_FLT_MIN)|, rounding down. 142 143 FLT_SMALLEST_EXP 144 REAL_FLT_MIN's exponent - EXP_BIAS (1 if denormals are not supported, 145 -22 if they are). 146 */ 147 148 #ifdef _FLT_NO_DENORMALS 149 #define FLT_UWORD_IS_ZERO(x) ((x)<0x00800000L) 150 #define FLT_UWORD_IS_SUBNORMAL(x) 0 151 #define FLT_UWORD_MIN 0x00800000 152 #define FLT_UWORD_EXP_MIN 0x42fc0000 153 #define FLT_UWORD_LOG_MIN 0x42aeac50 154 #define FLT_SMALLEST_EXP 1 155 #else 156 #define FLT_UWORD_IS_ZERO(x) ((x)==0) 157 #define FLT_UWORD_IS_SUBNORMAL(x) ((x)<0x00800000L) 158 #define FLT_UWORD_MIN 0x00000001 159 #define FLT_UWORD_EXP_MIN 0x43160000 160 #define FLT_UWORD_LOG_MIN 0x42cff1b5 161 #define FLT_SMALLEST_EXP -22 162 #endif 163 164 /* 165 * set X_TLOSS = pi*2**52, which is possibly defined in <values.h> 166 * (one may replace the following line by "#include <values.h>") 167 */ 168 169 #define X_TLOSS 1.41484755040568800000e+16 170 171 #ifdef _NEED_FLOAT64 172 extern __int32_t __rem_pio2 (__float64,__float64*); 173 174 /* fdlibm kernel function */ 175 extern __float64 __kernel_sin (__float64,__float64,int); 176 extern __float64 __kernel_cos (__float64,__float64); 177 extern __float64 __kernel_tan (__float64,__float64,int); 178 extern int __kernel_rem_pio2 (__float64*,__float64*,int,int,int,const __int32_t*); 179 #endif 180 181 extern __int32_t __rem_pio2f (float,float*); 182 183 /* float versions of fdlibm kernel functions */ 184 extern float __kernel_sinf (float,float,int); 185 extern float __kernel_cosf (float,float); 186 extern float __kernel_tanf (float,float,int); 187 extern int __kernel_rem_pio2f (float*,float*,int,int,int,const __int32_t*); 188 189 /* The original code used statements like 190 n0 = ((*(int*)&one)>>29)^1; * index of high word * 191 ix0 = *(n0+(int*)&x); * high word of x * 192 ix1 = *((1-n0)+(int*)&x); * low word of x * 193 to dig two 32 bit words out of the 64 bit IEEE floating point 194 value. That is non-ANSI, and, moreover, the gcc instruction 195 scheduler gets it wrong. We instead use the following macros. 196 Unlike the original code, we determine the endianness at compile 197 time, not at run time; I don't see much benefit to selecting 198 endianness at run time. */ 199 200 #ifndef __IEEE_BIG_ENDIAN 201 #ifndef __IEEE_LITTLE_ENDIAN 202 #error Must define endianness 203 #endif 204 #endif 205 206 /* A union which permits us to convert between a double and two 32 bit 207 ints. */ 208 209 #ifdef __IEEE_BIG_ENDIAN 210 211 typedef union 212 { 213 uint64_t bits; 214 struct 215 { 216 __uint32_t msw; 217 __uint32_t lsw; 218 } parts; 219 } ieee_double_shape_type; 220 221 #endif 222 223 #ifdef __IEEE_LITTLE_ENDIAN 224 225 typedef union 226 { 227 uint64_t bits; 228 struct 229 { 230 __uint32_t lsw; 231 __uint32_t msw; 232 } parts; 233 } ieee_double_shape_type; 234 235 #endif 236 237 /* Get two 32 bit ints from a double. */ 238 239 #define EXTRACT_WORDS(ix0,ix1,d) \ 240 do { \ 241 ieee_double_shape_type ew_u; \ 242 ew_u.bits = asuint64(d); \ 243 (ix0) = ew_u.parts.msw; \ 244 (ix1) = ew_u.parts.lsw; \ 245 } while (0) 246 247 /* Get the more significant 32 bit int from a double. */ 248 249 #define GET_HIGH_WORD(i,d) \ 250 do { \ 251 ieee_double_shape_type gh_u; \ 252 gh_u.bits = asuint64(d); \ 253 (i) = gh_u.parts.msw; \ 254 } while (0) 255 256 /* Get the less significant 32 bit int from a double. */ 257 258 #define GET_LOW_WORD(i,d) \ 259 do { \ 260 ieee_double_shape_type gl_u; \ 261 gl_u.bits = asuint64(d); \ 262 (i) = gl_u.parts.lsw; \ 263 } while (0) 264 265 /* Set a double from two 32 bit ints. */ 266 267 #define INSERT_WORDS(d,ix0,ix1) \ 268 do { \ 269 ieee_double_shape_type iw_u; \ 270 iw_u.parts.msw = (ix0); \ 271 iw_u.parts.lsw = (ix1); \ 272 (d) = asfloat64(iw_u.bits); \ 273 } while (0) 274 275 /* Set the more significant 32 bits of a double from an int. */ 276 277 #define SET_HIGH_WORD(d,v) \ 278 do { \ 279 ieee_double_shape_type sh_u; \ 280 sh_u.bits = asuint64(d); \ 281 sh_u.parts.msw = (v); \ 282 (d) = asfloat64(sh_u.bits); \ 283 } while (0) 284 285 /* Set the less significant 32 bits of a double from an int. */ 286 287 #define SET_LOW_WORD(d,v) \ 288 do { \ 289 ieee_double_shape_type sl_u; \ 290 sl_u.bits = asuint64(d); \ 291 sl_u.parts.lsw = (v); \ 292 (d) = asfloat64(sl_u.bits); \ 293 } while (0) 294 295 /* A union which permits us to convert between a float and a 32 bit 296 int. */ 297 298 /* Get a 32 bit int from a float. */ 299 300 #define GET_FLOAT_WORD(i,d) ((i) = asuint(d)) 301 302 /* Set a float from a 32 bit int. */ 303 304 #define SET_FLOAT_WORD(d,i) ((d) = asfloat(i)) 305 306 /* Macros to avoid undefined behaviour that can arise if the amount 307 of a shift is exactly equal to the size of the shifted operand. */ 308 309 #define SAFE_LEFT_SHIFT(op,amt) \ 310 (((amt) < (int) (8 * sizeof(op))) ? ((op) << (amt)) : 0) 311 312 #define SAFE_RIGHT_SHIFT(op,amt) \ 313 (((amt) < (int) (8 * sizeof(op))) ? ((op) >> (amt)) : 0) 314 315 #ifdef _COMPLEX_H 316 317 /* 318 * Quoting from ISO/IEC 9899:TC2: 319 * 320 * 6.2.5.13 Types 321 * Each complex type has the same representation and alignment requirements as 322 * an array type containing exactly two elements of the corresponding real type; 323 * the first element is equal to the real part, and the second element to the 324 * imaginary part, of the complex number. 325 */ 326 typedef union { 327 float complex z; 328 float parts[2]; 329 } float_complex; 330 331 typedef union { 332 double complex z; 333 double parts[2]; 334 } double_complex; 335 336 typedef union { 337 long double complex z; 338 long double parts[2]; 339 } long_double_complex; 340 341 #define REAL_PART(z) ((z).parts[0]) 342 #define IMAG_PART(z) ((z).parts[1]) 343 344 #endif /* _COMPLEX_H */ 345 346 #endif /* _FDLIBM_H_ */ 347
