1 /* 2 * ==================================================== 3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4 * 5 * Developed at SunPro, a Sun Microsystems, Inc. business. 6 * Permission to use, copy, modify, and distribute this 7 * software is freely granted, provided that this notice 8 * is preserved. 9 * ==================================================== 10 */ 11 12 /* 13 * from: @(#)fdlibm.h 5.1 93/09/24 14 * $FreeBSD$ 15 */ 16 17 #ifndef _MATH_PRIVATE_H_ 18 #define _MATH_PRIVATE_H_ 19 20 #include <sys/types.h> 21 #include <machine/endian.h> 22 23 /* 24 * The original fdlibm code used statements like: 25 * n0 = ((*(int*)&one)>>29)^1; * index of high word * 26 * ix0 = *(n0+(int*)&x); * high word of x * 27 * ix1 = *((1-n0)+(int*)&x); * low word of x * 28 * to dig two 32 bit words out of the 64 bit IEEE floating point 29 * value. That is non-ANSI, and, moreover, the gcc instruction 30 * scheduler gets it wrong. We instead use the following macros. 31 * Unlike the original code, we determine the endianness at compile 32 * time, not at run time; I don't see much benefit to selecting 33 * endianness at run time. 34 */ 35 36 /* 37 * A union which permits us to convert between a double and two 32 bit 38 * ints. 39 */ 40 41 #ifdef __arm__ 42 #if defined(__VFP_FP__) 43 #define IEEE_WORD_ORDER BYTE_ORDER 44 #else 45 #define IEEE_WORD_ORDER BIG_ENDIAN 46 #endif 47 #else /* __arm__ */ 48 #define IEEE_WORD_ORDER BYTE_ORDER 49 #endif 50 51 #if IEEE_WORD_ORDER == BIG_ENDIAN 52 53 typedef union 54 { 55 double value; 56 struct 57 { 58 u_int32_t msw; 59 u_int32_t lsw; 60 } parts; 61 } ieee_double_shape_type; 62 63 #endif 64 65 #if IEEE_WORD_ORDER == LITTLE_ENDIAN 66 67 typedef union 68 { 69 double value; 70 struct 71 { 72 u_int32_t lsw; 73 u_int32_t msw; 74 } parts; 75 } ieee_double_shape_type; 76 77 #endif 78 79 /* Get two 32 bit ints from a double. */ 80 81 #define EXTRACT_WORDS(ix0,ix1,d) \ 82 do { \ 83 ieee_double_shape_type ew_u; \ 84 ew_u.value = (d); \ 85 (ix0) = ew_u.parts.msw; \ 86 (ix1) = ew_u.parts.lsw; \ 87 } while (0) 88 89 /* Get the more significant 32 bit int from a double. */ 90 91 #define GET_HIGH_WORD(i,d) \ 92 do { \ 93 ieee_double_shape_type gh_u; \ 94 gh_u.value = (d); \ 95 (i) = gh_u.parts.msw; \ 96 } while (0) 97 98 /* Get the less significant 32 bit int from a double. */ 99 100 #define GET_LOW_WORD(i,d) \ 101 do { \ 102 ieee_double_shape_type gl_u; \ 103 gl_u.value = (d); \ 104 (i) = gl_u.parts.lsw; \ 105 } while (0) 106 107 /* Set a double from two 32 bit ints. */ 108 109 #define INSERT_WORDS(d,ix0,ix1) \ 110 do { \ 111 ieee_double_shape_type iw_u; \ 112 iw_u.parts.msw = (ix0); \ 113 iw_u.parts.lsw = (ix1); \ 114 (d) = iw_u.value; \ 115 } while (0) 116 117 /* Set the more significant 32 bits of a double from an int. */ 118 119 #define SET_HIGH_WORD(d,v) \ 120 do { \ 121 ieee_double_shape_type sh_u; \ 122 sh_u.value = (d); \ 123 sh_u.parts.msw = (v); \ 124 (d) = sh_u.value; \ 125 } while (0) 126 127 /* Set the less significant 32 bits of a double from an int. */ 128 129 #define SET_LOW_WORD(d,v) \ 130 do { \ 131 ieee_double_shape_type sl_u; \ 132 sl_u.value = (d); \ 133 sl_u.parts.lsw = (v); \ 134 (d) = sl_u.value; \ 135 } while (0) 136 137 /* 138 * A union which permits us to convert between a float and a 32 bit 139 * int. 140 */ 141 142 typedef union 143 { 144 float value; 145 /* FIXME: Assumes 32 bit int. */ 146 unsigned int word; 147 } ieee_float_shape_type; 148 149 /* Get a 32 bit int from a float. */ 150 151 #define GET_FLOAT_WORD(i,d) \ 152 do { \ 153 ieee_float_shape_type gf_u; \ 154 gf_u.value = (d); \ 155 (i) = gf_u.word; \ 156 } while (0) 157 158 /* Set a float from a 32 bit int. */ 159 160 #define SET_FLOAT_WORD(d,i) \ 161 do { \ 162 ieee_float_shape_type sf_u; \ 163 sf_u.word = (i); \ 164 (d) = sf_u.value; \ 165 } while (0) 166 167 #ifdef FLT_EVAL_METHOD 168 /* 169 * Attempt to get strict C99 semantics for assignment with non-C99 compilers. 170 */ 171 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 172 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) 173 #else 174 #define STRICT_ASSIGN(type, lval, rval) do { \ 175 volatile type __lval; \ 176 \ 177 if (sizeof(type) >= sizeof(double)) \ 178 (lval) = (rval); \ 179 else { \ 180 __lval = (rval); \ 181 (lval) = __lval; \ 182 } \ 183 } while (0) 184 #endif 185 #endif 186 187 /* 188 * Common routine to process the arguments to nan(), nanf(), and nanl(). 189 */ 190 void _scan_nan(uint32_t *__words, int __num_words, const char *__s); 191 192 #ifdef _COMPLEX_H 193 194 /* 195 * C99 specifies that complex numbers have the same representation as 196 * an array of two elements, where the first element is the real part 197 * and the second element is the imaginary part. 198 */ 199 typedef union { 200 float complex f; 201 float a[2]; 202 } float_complex; 203 typedef union { 204 double complex f; 205 double a[2]; 206 } double_complex; 207 typedef union { 208 long double complex f; 209 long double a[2]; 210 } long_double_complex; 211 #define REALPART(z) ((z).a[0]) 212 #define IMAGPART(z) ((z).a[1]) 213 214 /* 215 * Inline functions that can be used to construct complex values. 216 * 217 * The C99 standard intends x+I*y to be used for this, but x+I*y is 218 * currently unusable in general since gcc introduces many overflow, 219 * underflow, sign and efficiency bugs by rewriting I*y as 220 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. 221 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted 222 * to -0.0+I*0.0. 223 */ 224 static __inline float complex 225 cpackf(float x, float y) 226 { 227 float_complex z; 228 229 REALPART(z) = x; 230 IMAGPART(z) = y; 231 return (z.f); 232 } 233 234 static __inline double complex 235 cpack(double x, double y) 236 { 237 double_complex z; 238 239 REALPART(z) = x; 240 IMAGPART(z) = y; 241 return (z.f); 242 } 243 244 static __inline long double complex 245 cpackl(long double x, long double y) 246 { 247 long_double_complex z; 248 249 REALPART(z) = x; 250 IMAGPART(z) = y; 251 return (z.f); 252 } 253 #endif /* _COMPLEX_H */ 254 255 #ifdef __GNUCLIKE_ASM 256 257 /* Asm versions of some functions. */ 258 259 #ifdef __amd64__ 260 static __inline int 261 irint(double x) 262 { 263 int n; 264 265 asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); 266 return (n); 267 } 268 #define HAVE_EFFICIENT_IRINT 269 #endif 270 271 #ifdef __i386__ 272 static __inline int 273 irint(double x) 274 { 275 int n; 276 277 asm("fistl %0" : "=m" (n) : "t" (x)); 278 return (n); 279 } 280 #define HAVE_EFFICIENT_IRINT 281 #endif 282 283 #endif /* __GNUCLIKE_ASM */ 284 285 /* 286 * ieee style elementary functions 287 * 288 * We rename functions here to improve other sources' diffability 289 * against fdlibm. 290 */ 291 #define __ieee754_sqrt sqrt 292 #define __ieee754_acos acos 293 #define __ieee754_acosh acosh 294 #define __ieee754_log log 295 #define __ieee754_atanh atanh 296 #define __ieee754_asin asin 297 #define __ieee754_atan2 atan2 298 #define __ieee754_exp exp 299 #define __ieee754_cosh cosh 300 #define __ieee754_fmod fmod 301 #define __ieee754_pow pow 302 #define __ieee754_lgamma lgamma 303 #define __ieee754_gamma gamma 304 #define __ieee754_lgamma_r lgamma_r 305 #define __ieee754_gamma_r gamma_r 306 #define __ieee754_log10 log10 307 #define __ieee754_sinh sinh 308 #define __ieee754_hypot hypot 309 #define __ieee754_j0 j0 310 #define __ieee754_j1 j1 311 #define __ieee754_y0 y0 312 #define __ieee754_y1 y1 313 #define __ieee754_jn jn 314 #define __ieee754_yn yn 315 #define __ieee754_remainder remainder 316 #define __ieee754_scalb scalb 317 #define __ieee754_sqrtf sqrtf 318 #define __ieee754_acosf acosf 319 #define __ieee754_acoshf acoshf 320 #define __ieee754_logf logf 321 #define __ieee754_atanhf atanhf 322 #define __ieee754_asinf asinf 323 #define __ieee754_atan2f atan2f 324 #define __ieee754_expf expf 325 #define __ieee754_coshf coshf 326 #define __ieee754_fmodf fmodf 327 #define __ieee754_powf powf 328 #define __ieee754_lgammaf lgammaf 329 #define __ieee754_gammaf gammaf 330 #define __ieee754_lgammaf_r lgammaf_r 331 #define __ieee754_gammaf_r gammaf_r 332 #define __ieee754_log10f log10f 333 #define __ieee754_sinhf sinhf 334 #define __ieee754_hypotf hypotf 335 #define __ieee754_j0f j0f 336 #define __ieee754_j1f j1f 337 #define __ieee754_y0f y0f 338 #define __ieee754_y1f y1f 339 #define __ieee754_jnf jnf 340 #define __ieee754_ynf ynf 341 #define __ieee754_remainderf remainderf 342 #define __ieee754_scalbf scalbf 343 344 /* fdlibm kernel function */ 345 int __kernel_rem_pio2(double*,double*,int,int,int); 346 347 /* double precision kernel functions */ 348 #ifdef INLINE_REM_PIO2 349 __inline 350 #endif 351 int __ieee754_rem_pio2(double,double*); 352 double __kernel_sin(double,double,int); 353 double __kernel_cos(double,double); 354 double __kernel_tan(double,double,int); 355 356 /* float precision kernel functions */ 357 #ifdef INLINE_REM_PIO2F 358 __inline 359 #endif 360 int __ieee754_rem_pio2f(float,double*); 361 #ifdef INLINE_KERNEL_SINDF 362 __inline 363 #endif 364 float __kernel_sindf(double); 365 #ifdef INLINE_KERNEL_COSDF 366 __inline 367 #endif 368 float __kernel_cosdf(double); 369 #ifdef INLINE_KERNEL_TANDF 370 __inline 371 #endif 372 float __kernel_tandf(double,int); 373 374 /* long double precision kernel functions */ 375 long double __kernel_sinl(long double, long double, int); 376 long double __kernel_cosl(long double, long double); 377 long double __kernel_tanl(long double, long double, int); 378 379 #endif /* !_MATH_PRIVATE_H_ */ 380