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 143extern double scalb __P((double, int)); 144#else 145extern double scalb __P((double, double)); 146#endif 147extern double significand __P((double)); 148 149/* ieee style elementary functions */ 150extern double __ieee754_sqrt __P((double)); 151extern double __ieee754_acos __P((double)); 152extern double __ieee754_acosh __P((double)); 153extern double __ieee754_log __P((double)); 154extern double __ieee754_atanh __P((double)); 155extern double __ieee754_asin __P((double)); 156extern double __ieee754_atan2 __P((double,double)); 157extern double __ieee754_exp __P((double)); 158extern double __ieee754_cosh __P((double)); 159extern double __ieee754_fmod __P((double,double)); 160extern double __ieee754_pow __P((double,double)); 161extern double __ieee754_lgamma_r __P((double,int *)); 162extern double __ieee754_gamma_r __P((double,int *)); 163extern double __ieee754_log10 __P((double)); 164extern double __ieee754_sinh __P((double)); 165extern double __ieee754_hypot __P((double,double)); 166extern double __ieee754_j0 __P((double)); 167extern double __ieee754_j1 __P((double)); 168extern double __ieee754_y0 __P((double)); 169extern double __ieee754_y1 __P((double)); 170extern double __ieee754_jn __P((int,double)); 171extern double __ieee754_yn __P((int,double)); 172extern double __ieee754_remainder __P((double,double)); 173extern __int32_t __ieee754_rem_pio2 __P((double,double*)); 174#ifdef _SCALB_INT 175extern double __ieee754_scalb __P((double,int)); 176#else 177extern double __ieee754_scalb __P((double,double)); 178#endif 179 180/* fdlibm kernel function */ 181extern double __kernel_standard __P((double,double,int)); 182extern double __kernel_sin __P((double,double,int)); 183extern double __kernel_cos __P((double,double)); 184extern double __kernel_tan __P((double,double,int)); 185extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const __int32_t*)); 186 187/* Undocumented float functions. */ 188#ifdef _SCALB_INT 189extern float scalbf __P((float, int)); 190#else 191extern float scalbf __P((float, float)); 192#endif 193extern float significandf __P((float)); 194 195/* ieee style elementary float functions */ 196extern float __ieee754_sqrtf __P((float)); 197extern float __ieee754_acosf __P((float)); 198extern float __ieee754_acoshf __P((float)); 199extern float __ieee754_logf __P((float)); 200extern float __ieee754_atanhf __P((float)); 201extern float __ieee754_asinf __P((float)); 202extern float __ieee754_atan2f __P((float,float)); 203extern float __ieee754_expf __P((float)); 204extern float __ieee754_coshf __P((float)); 205extern float __ieee754_fmodf __P((float,float)); 206extern float __ieee754_powf __P((float,float)); 207extern float __ieee754_lgammaf_r __P((float,int *)); 208extern float __ieee754_gammaf_r __P((float,int *)); 209extern float __ieee754_log10f __P((float)); 210extern float __ieee754_sinhf __P((float)); 211extern float __ieee754_hypotf __P((float,float)); 212extern float __ieee754_j0f __P((float)); 213extern float __ieee754_j1f __P((float)); 214extern float __ieee754_y0f __P((float)); 215extern float __ieee754_y1f __P((float)); 216extern float __ieee754_jnf __P((int,float)); 217extern float __ieee754_ynf __P((int,float)); 218extern float __ieee754_remainderf __P((float,float)); 219extern __int32_t __ieee754_rem_pio2f __P((float,float*)); 220#ifdef _SCALB_INT 221extern float __ieee754_scalbf __P((float,int)); 222#else 223extern float __ieee754_scalbf __P((float,float)); 224#endif 225 226/* float versions of fdlibm kernel functions */ 227extern float __kernel_sinf __P((float,float,int)); 228extern float __kernel_cosf __P((float,float)); 229extern float __kernel_tanf __P((float,float,int)); 230extern 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 254typedef 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 268typedef 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) \ 283do { \ 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) \ 293do { \ 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) \ 302do { \ 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) \ 311do { \ 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) \ 321do { \ 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) \ 331do { \ 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 341typedef 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) \ 350do { \ 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) \ 359do { \ 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 */ 385typedef union { 386 float complex z; 387 float parts[2]; 388} float_complex; 389 390typedef union { 391 double complex z; 392 double parts[2]; 393} double_complex; 394 395typedef 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