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 extern __int32_t __rem_pio2 (__float64,__float64*); 172 173 /* fdlibm kernel function */ 174 extern __float64 __kernel_sin (__float64,__float64,int); 175 extern __float64 __kernel_cos (__float64,__float64); 176 extern __float64 __kernel_tan (__float64,__float64,int); 177 extern int __kernel_rem_pio2 (__float64*,__float64*,int,int,int,const __int32_t*); 178 179 extern __int32_t __rem_pio2f (float,float*); 180 181 /* float versions of fdlibm kernel functions */ 182 extern float __kernel_sinf (float,float,int); 183 extern float __kernel_cosf (float,float); 184 extern float __kernel_tanf (float,float,int); 185 extern int __kernel_rem_pio2f (float*,float*,int,int,int,const __int32_t*); 186 187 /* The original code used statements like 188 n0 = ((*(int*)&one)>>29)^1; * index of high word * 189 ix0 = *(n0+(int*)&x); * high word of x * 190 ix1 = *((1-n0)+(int*)&x); * low word of x * 191 to dig two 32 bit words out of the 64 bit IEEE floating point 192 value. That is non-ANSI, and, moreover, the gcc instruction 193 scheduler gets it wrong. We instead use the following macros. 194 Unlike the original code, we determine the endianness at compile 195 time, not at run time; I don't see much benefit to selecting 196 endianness at run time. */ 197 198 #ifndef __IEEE_BIG_ENDIAN 199 #ifndef __IEEE_LITTLE_ENDIAN 200 #error Must define endianness 201 #endif 202 #endif 203 204 /* A union which permits us to convert between a double and two 32 bit 205 ints. */ 206 207 #ifdef __IEEE_BIG_ENDIAN 208 209 typedef union 210 { 211 uint64_t bits; 212 struct 213 { 214 __uint32_t msw; 215 __uint32_t lsw; 216 } parts; 217 } ieee_double_shape_type; 218 219 #endif 220 221 #ifdef __IEEE_LITTLE_ENDIAN 222 223 typedef union 224 { 225 uint64_t bits; 226 struct 227 { 228 __uint32_t lsw; 229 __uint32_t msw; 230 } parts; 231 } ieee_double_shape_type; 232 233 #endif 234 235 /* Get two 32 bit ints from a double. */ 236 237 #define EXTRACT_WORDS(ix0,ix1,d) \ 238 do { \ 239 ieee_double_shape_type ew_u; \ 240 ew_u.bits = asuint64(d); \ 241 (ix0) = ew_u.parts.msw; \ 242 (ix1) = ew_u.parts.lsw; \ 243 } while (0) 244 245 /* Get the more significant 32 bit int from a double. */ 246 247 #define GET_HIGH_WORD(i,d) \ 248 do { \ 249 ieee_double_shape_type gh_u; \ 250 gh_u.bits = asuint64(d); \ 251 (i) = gh_u.parts.msw; \ 252 } while (0) 253 254 /* Get the less significant 32 bit int from a double. */ 255 256 #define GET_LOW_WORD(i,d) \ 257 do { \ 258 ieee_double_shape_type gl_u; \ 259 gl_u.bits = asuint64(d); \ 260 (i) = gl_u.parts.lsw; \ 261 } while (0) 262 263 /* Set a double from two 32 bit ints. */ 264 265 #define INSERT_WORDS(d,ix0,ix1) \ 266 do { \ 267 ieee_double_shape_type iw_u; \ 268 iw_u.parts.msw = (ix0); \ 269 iw_u.parts.lsw = (ix1); \ 270 (d) = asdouble(iw_u.bits); \ 271 } while (0) 272 273 /* Set the more significant 32 bits of a double from an int. */ 274 275 #define SET_HIGH_WORD(d,v) \ 276 do { \ 277 ieee_double_shape_type sh_u; \ 278 sh_u.bits = asuint64(d); \ 279 sh_u.parts.msw = (v); \ 280 (d) = asdouble(sh_u.bits); \ 281 } while (0) 282 283 /* Set the less significant 32 bits of a double from an int. */ 284 285 #define SET_LOW_WORD(d,v) \ 286 do { \ 287 ieee_double_shape_type sl_u; \ 288 sl_u.bits = asuint64(d); \ 289 sl_u.parts.lsw = (v); \ 290 (d) = asdouble(sl_u.bits); \ 291 } while (0) 292 293 /* A union which permits us to convert between a float and a 32 bit 294 int. */ 295 296 /* Get a 32 bit int from a float. */ 297 298 #define GET_FLOAT_WORD(i,d) ((i) = asuint(d)) 299 300 /* Set a float from a 32 bit int. */ 301 302 #define SET_FLOAT_WORD(d,i) ((d) = asfloat(i)) 303 304 /* Macros to avoid undefined behaviour that can arise if the amount 305 of a shift is exactly equal to the size of the shifted operand. */ 306 307 #define SAFE_LEFT_SHIFT(op,amt) \ 308 (((amt) < (int) (8 * sizeof(op))) ? ((op) << (amt)) : 0) 309 310 #define SAFE_RIGHT_SHIFT(op,amt) \ 311 (((amt) < (int) (8 * sizeof(op))) ? ((op) >> (amt)) : 0) 312 313 #ifdef _COMPLEX_H 314 315 /* 316 * Quoting from ISO/IEC 9899:TC2: 317 * 318 * 6.2.5.13 Types 319 * Each complex type has the same representation and alignment requirements as 320 * an array type containing exactly two elements of the corresponding real type; 321 * the first element is equal to the real part, and the second element to the 322 * imaginary part, of the complex number. 323 */ 324 typedef union { 325 float complex z; 326 float parts[2]; 327 } float_complex; 328 329 typedef union { 330 double complex z; 331 double parts[2]; 332 } double_complex; 333 334 typedef union { 335 long double complex z; 336 long double parts[2]; 337 } long_double_complex; 338 339 #define REAL_PART(z) ((z).parts[0]) 340 #define IMAG_PART(z) ((z).parts[1]) 341 342 #endif /* _COMPLEX_H */ 343 344 #endif /* _FDLIBM_H_ */ 345
