13a8617a8SJordan K. Hubbard /* 23a8617a8SJordan K. Hubbard * ==================================================== 33a8617a8SJordan K. Hubbard * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 43a8617a8SJordan K. Hubbard * 53a8617a8SJordan K. Hubbard * Developed at SunPro, a Sun Microsystems, Inc. business. 63a8617a8SJordan K. Hubbard * Permission to use, copy, modify, and distribute this 73a8617a8SJordan K. Hubbard * software is freely granted, provided that this notice 83a8617a8SJordan K. Hubbard * is preserved. 93a8617a8SJordan K. Hubbard * ==================================================== 103a8617a8SJordan K. Hubbard */ 113a8617a8SJordan K. Hubbard 123a8617a8SJordan K. Hubbard /* 133a8617a8SJordan K. Hubbard * from: @(#)fdlibm.h 5.1 93/09/24 147f3dea24SPeter Wemm * $FreeBSD$ 153a8617a8SJordan K. Hubbard */ 163a8617a8SJordan K. Hubbard 173a8617a8SJordan K. Hubbard #ifndef _MATH_PRIVATE_H_ 183a8617a8SJordan K. Hubbard #define _MATH_PRIVATE_H_ 193a8617a8SJordan K. Hubbard 203a8617a8SJordan K. Hubbard #include <sys/types.h> 21bcfa1759SBrian Somers #include <machine/endian.h> 223a8617a8SJordan K. Hubbard 23ef1ee63eSAlexey Zelkin /* 24ef1ee63eSAlexey Zelkin * The original fdlibm code used statements like: 25ef1ee63eSAlexey Zelkin * n0 = ((*(int*)&one)>>29)^1; * index of high word * 26ef1ee63eSAlexey Zelkin * ix0 = *(n0+(int*)&x); * high word of x * 27ef1ee63eSAlexey Zelkin * ix1 = *((1-n0)+(int*)&x); * low word of x * 28ef1ee63eSAlexey Zelkin * to dig two 32 bit words out of the 64 bit IEEE floating point 29ef1ee63eSAlexey Zelkin * value. That is non-ANSI, and, moreover, the gcc instruction 30ef1ee63eSAlexey Zelkin * scheduler gets it wrong. We instead use the following macros. 31ef1ee63eSAlexey Zelkin * Unlike the original code, we determine the endianness at compile 32ef1ee63eSAlexey Zelkin * time, not at run time; I don't see much benefit to selecting 33ef1ee63eSAlexey Zelkin * endianness at run time. 34ef1ee63eSAlexey Zelkin */ 353a8617a8SJordan K. Hubbard 36ef1ee63eSAlexey Zelkin /* 37ef1ee63eSAlexey Zelkin * A union which permits us to convert between a double and two 32 bit 38ef1ee63eSAlexey Zelkin * ints. 39ef1ee63eSAlexey Zelkin */ 403a8617a8SJordan K. Hubbard 413a8617a8SJordan K. Hubbard #if BYTE_ORDER == BIG_ENDIAN 423a8617a8SJordan K. Hubbard 433a8617a8SJordan K. Hubbard typedef union 443a8617a8SJordan K. Hubbard { 453a8617a8SJordan K. Hubbard double value; 463a8617a8SJordan K. Hubbard struct 473a8617a8SJordan K. Hubbard { 483a8617a8SJordan K. Hubbard u_int32_t msw; 493a8617a8SJordan K. Hubbard u_int32_t lsw; 503a8617a8SJordan K. Hubbard } parts; 513a8617a8SJordan K. Hubbard } ieee_double_shape_type; 523a8617a8SJordan K. Hubbard 533a8617a8SJordan K. Hubbard #endif 543a8617a8SJordan K. Hubbard 553a8617a8SJordan K. Hubbard #if BYTE_ORDER == LITTLE_ENDIAN 563a8617a8SJordan K. Hubbard 573a8617a8SJordan K. Hubbard typedef union 583a8617a8SJordan K. Hubbard { 593a8617a8SJordan K. Hubbard double value; 603a8617a8SJordan K. Hubbard struct 613a8617a8SJordan K. Hubbard { 623a8617a8SJordan K. Hubbard u_int32_t lsw; 633a8617a8SJordan K. Hubbard u_int32_t msw; 643a8617a8SJordan K. Hubbard } parts; 653a8617a8SJordan K. Hubbard } ieee_double_shape_type; 663a8617a8SJordan K. Hubbard 673a8617a8SJordan K. Hubbard #endif 683a8617a8SJordan K. Hubbard 693a8617a8SJordan K. Hubbard /* Get two 32 bit ints from a double. */ 703a8617a8SJordan K. Hubbard 713a8617a8SJordan K. Hubbard #define EXTRACT_WORDS(ix0,ix1,d) \ 723a8617a8SJordan K. Hubbard do { \ 733a8617a8SJordan K. Hubbard ieee_double_shape_type ew_u; \ 743a8617a8SJordan K. Hubbard ew_u.value = (d); \ 753a8617a8SJordan K. Hubbard (ix0) = ew_u.parts.msw; \ 763a8617a8SJordan K. Hubbard (ix1) = ew_u.parts.lsw; \ 773a8617a8SJordan K. Hubbard } while (0) 783a8617a8SJordan K. Hubbard 793a8617a8SJordan K. Hubbard /* Get the more significant 32 bit int from a double. */ 803a8617a8SJordan K. Hubbard 813a8617a8SJordan K. Hubbard #define GET_HIGH_WORD(i,d) \ 823a8617a8SJordan K. Hubbard do { \ 833a8617a8SJordan K. Hubbard ieee_double_shape_type gh_u; \ 843a8617a8SJordan K. Hubbard gh_u.value = (d); \ 853a8617a8SJordan K. Hubbard (i) = gh_u.parts.msw; \ 863a8617a8SJordan K. Hubbard } while (0) 873a8617a8SJordan K. Hubbard 883a8617a8SJordan K. Hubbard /* Get the less significant 32 bit int from a double. */ 893a8617a8SJordan K. Hubbard 903a8617a8SJordan K. Hubbard #define GET_LOW_WORD(i,d) \ 913a8617a8SJordan K. Hubbard do { \ 923a8617a8SJordan K. Hubbard ieee_double_shape_type gl_u; \ 933a8617a8SJordan K. Hubbard gl_u.value = (d); \ 943a8617a8SJordan K. Hubbard (i) = gl_u.parts.lsw; \ 953a8617a8SJordan K. Hubbard } while (0) 963a8617a8SJordan K. Hubbard 973a8617a8SJordan K. Hubbard /* Set a double from two 32 bit ints. */ 983a8617a8SJordan K. Hubbard 993a8617a8SJordan K. Hubbard #define INSERT_WORDS(d,ix0,ix1) \ 1003a8617a8SJordan K. Hubbard do { \ 1013a8617a8SJordan K. Hubbard ieee_double_shape_type iw_u; \ 1023a8617a8SJordan K. Hubbard iw_u.parts.msw = (ix0); \ 1033a8617a8SJordan K. Hubbard iw_u.parts.lsw = (ix1); \ 1043a8617a8SJordan K. Hubbard (d) = iw_u.value; \ 1053a8617a8SJordan K. Hubbard } while (0) 1063a8617a8SJordan K. Hubbard 1073a8617a8SJordan K. Hubbard /* Set the more significant 32 bits of a double from an int. */ 1083a8617a8SJordan K. Hubbard 1093a8617a8SJordan K. Hubbard #define SET_HIGH_WORD(d,v) \ 1103a8617a8SJordan K. Hubbard do { \ 1113a8617a8SJordan K. Hubbard ieee_double_shape_type sh_u; \ 1123a8617a8SJordan K. Hubbard sh_u.value = (d); \ 1133a8617a8SJordan K. Hubbard sh_u.parts.msw = (v); \ 1143a8617a8SJordan K. Hubbard (d) = sh_u.value; \ 1153a8617a8SJordan K. Hubbard } while (0) 1163a8617a8SJordan K. Hubbard 1173a8617a8SJordan K. Hubbard /* Set the less significant 32 bits of a double from an int. */ 1183a8617a8SJordan K. Hubbard 1193a8617a8SJordan K. Hubbard #define SET_LOW_WORD(d,v) \ 1203a8617a8SJordan K. Hubbard do { \ 1213a8617a8SJordan K. Hubbard ieee_double_shape_type sl_u; \ 1223a8617a8SJordan K. Hubbard sl_u.value = (d); \ 1233a8617a8SJordan K. Hubbard sl_u.parts.lsw = (v); \ 1243a8617a8SJordan K. Hubbard (d) = sl_u.value; \ 1253a8617a8SJordan K. Hubbard } while (0) 1263a8617a8SJordan K. Hubbard 127ef1ee63eSAlexey Zelkin /* 128ef1ee63eSAlexey Zelkin * A union which permits us to convert between a float and a 32 bit 129ef1ee63eSAlexey Zelkin * int. 130ef1ee63eSAlexey Zelkin */ 1313a8617a8SJordan K. Hubbard 1323a8617a8SJordan K. Hubbard typedef union 1333a8617a8SJordan K. Hubbard { 1343a8617a8SJordan K. Hubbard float value; 1353a8617a8SJordan K. Hubbard /* FIXME: Assumes 32 bit int. */ 1363a8617a8SJordan K. Hubbard unsigned int word; 1373a8617a8SJordan K. Hubbard } ieee_float_shape_type; 1383a8617a8SJordan K. Hubbard 1393a8617a8SJordan K. Hubbard /* Get a 32 bit int from a float. */ 1403a8617a8SJordan K. Hubbard 1413a8617a8SJordan K. Hubbard #define GET_FLOAT_WORD(i,d) \ 1423a8617a8SJordan K. Hubbard do { \ 1433a8617a8SJordan K. Hubbard ieee_float_shape_type gf_u; \ 1443a8617a8SJordan K. Hubbard gf_u.value = (d); \ 1453a8617a8SJordan K. Hubbard (i) = gf_u.word; \ 1463a8617a8SJordan K. Hubbard } while (0) 1473a8617a8SJordan K. Hubbard 1483a8617a8SJordan K. Hubbard /* Set a float from a 32 bit int. */ 1493a8617a8SJordan K. Hubbard 1503a8617a8SJordan K. Hubbard #define SET_FLOAT_WORD(d,i) \ 1513a8617a8SJordan K. Hubbard do { \ 1523a8617a8SJordan K. Hubbard ieee_float_shape_type sf_u; \ 1533a8617a8SJordan K. Hubbard sf_u.word = (i); \ 1543a8617a8SJordan K. Hubbard (d) = sf_u.value; \ 1553a8617a8SJordan K. Hubbard } while (0) 1563a8617a8SJordan K. Hubbard 15719b114daSBruce Evans #ifdef _COMPLEX_H 15819b114daSBruce Evans /* 15919b114daSBruce Evans * Inline functions that can be used to construct complex values. 16019b114daSBruce Evans * 16119b114daSBruce Evans * The C99 standard intends x+I*y to be used for this, but x+I*y is 16219b114daSBruce Evans * currently unusable in general since gcc introduces many overflow, 16319b114daSBruce Evans * underflow, sign and efficiency bugs by rewriting I*y as 16419b114daSBruce Evans * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. 16519b114daSBruce Evans * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted 16619b114daSBruce Evans * to -0.0+I*0.0. 16719b114daSBruce Evans */ 16819b114daSBruce Evans static __inline float complex 16919b114daSBruce Evans cpackf(float x, float y) 17019b114daSBruce Evans { 17119b114daSBruce Evans float complex z; 17219b114daSBruce Evans 17319b114daSBruce Evans __real__ z = x; 17419b114daSBruce Evans __imag__ z = y; 17519b114daSBruce Evans return (z); 17619b114daSBruce Evans } 17719b114daSBruce Evans 17819b114daSBruce Evans static __inline double complex 17919b114daSBruce Evans cpack(double x, double y) 18019b114daSBruce Evans { 18119b114daSBruce Evans double complex z; 18219b114daSBruce Evans 18319b114daSBruce Evans __real__ z = x; 18419b114daSBruce Evans __imag__ z = y; 18519b114daSBruce Evans return (z); 18619b114daSBruce Evans } 18719b114daSBruce Evans 18819b114daSBruce Evans static __inline long double complex 18919b114daSBruce Evans cpackl(long double x, long double y) 19019b114daSBruce Evans { 19119b114daSBruce Evans long double complex z; 19219b114daSBruce Evans 19319b114daSBruce Evans __real__ z = x; 19419b114daSBruce Evans __imag__ z = y; 19519b114daSBruce Evans return (z); 19619b114daSBruce Evans } 19719b114daSBruce Evans #endif /* _COMPLEX_H */ 19819b114daSBruce Evans 199e1b61b5bSDavid Schultz /* 200e1b61b5bSDavid Schultz * ieee style elementary functions 201e1b61b5bSDavid Schultz * 202e1b61b5bSDavid Schultz * We rename functions here to improve other sources' diffability 203e1b61b5bSDavid Schultz * against fdlibm. 204e1b61b5bSDavid Schultz */ 205e1b61b5bSDavid Schultz #define __ieee754_sqrt sqrt 206e1b61b5bSDavid Schultz #define __ieee754_acos acos 207e1b61b5bSDavid Schultz #define __ieee754_acosh acosh 208e1b61b5bSDavid Schultz #define __ieee754_log log 209e1b61b5bSDavid Schultz #define __ieee754_atanh atanh 210e1b61b5bSDavid Schultz #define __ieee754_asin asin 211e1b61b5bSDavid Schultz #define __ieee754_atan2 atan2 212e1b61b5bSDavid Schultz #define __ieee754_exp exp 213e1b61b5bSDavid Schultz #define __ieee754_cosh cosh 214e1b61b5bSDavid Schultz #define __ieee754_fmod fmod 215e1b61b5bSDavid Schultz #define __ieee754_pow pow 216e1b61b5bSDavid Schultz #define __ieee754_lgamma lgamma 217e1b61b5bSDavid Schultz #define __ieee754_gamma gamma 218e1b61b5bSDavid Schultz #define __ieee754_lgamma_r lgamma_r 219e1b61b5bSDavid Schultz #define __ieee754_gamma_r gamma_r 220e1b61b5bSDavid Schultz #define __ieee754_log10 log10 221e1b61b5bSDavid Schultz #define __ieee754_sinh sinh 222e1b61b5bSDavid Schultz #define __ieee754_hypot hypot 223e1b61b5bSDavid Schultz #define __ieee754_j0 j0 224e1b61b5bSDavid Schultz #define __ieee754_j1 j1 225e1b61b5bSDavid Schultz #define __ieee754_y0 y0 226e1b61b5bSDavid Schultz #define __ieee754_y1 y1 227e1b61b5bSDavid Schultz #define __ieee754_jn jn 228e1b61b5bSDavid Schultz #define __ieee754_yn yn 229e1b61b5bSDavid Schultz #define __ieee754_remainder remainder 230e1b61b5bSDavid Schultz #define __ieee754_scalb scalb 231e1b61b5bSDavid Schultz #define __ieee754_sqrtf sqrtf 232e1b61b5bSDavid Schultz #define __ieee754_acosf acosf 233e1b61b5bSDavid Schultz #define __ieee754_acoshf acoshf 234e1b61b5bSDavid Schultz #define __ieee754_logf logf 235e1b61b5bSDavid Schultz #define __ieee754_atanhf atanhf 236e1b61b5bSDavid Schultz #define __ieee754_asinf asinf 237e1b61b5bSDavid Schultz #define __ieee754_atan2f atan2f 238e1b61b5bSDavid Schultz #define __ieee754_expf expf 239e1b61b5bSDavid Schultz #define __ieee754_coshf coshf 240e1b61b5bSDavid Schultz #define __ieee754_fmodf fmodf 241e1b61b5bSDavid Schultz #define __ieee754_powf powf 242e1b61b5bSDavid Schultz #define __ieee754_lgammaf lgammaf 243e1b61b5bSDavid Schultz #define __ieee754_gammaf gammaf 244e1b61b5bSDavid Schultz #define __ieee754_lgammaf_r lgammaf_r 245e1b61b5bSDavid Schultz #define __ieee754_gammaf_r gammaf_r 246e1b61b5bSDavid Schultz #define __ieee754_log10f log10f 247e1b61b5bSDavid Schultz #define __ieee754_sinhf sinhf 248e1b61b5bSDavid Schultz #define __ieee754_hypotf hypotf 249e1b61b5bSDavid Schultz #define __ieee754_j0f j0f 250e1b61b5bSDavid Schultz #define __ieee754_j1f j1f 251e1b61b5bSDavid Schultz #define __ieee754_y0f y0f 252e1b61b5bSDavid Schultz #define __ieee754_y1f y1f 253e1b61b5bSDavid Schultz #define __ieee754_jnf jnf 254e1b61b5bSDavid Schultz #define __ieee754_ynf ynf 255e1b61b5bSDavid Schultz #define __ieee754_remainderf remainderf 256e1b61b5bSDavid Schultz #define __ieee754_scalbf scalbf 2573a8617a8SJordan K. Hubbard 2583a8617a8SJordan K. Hubbard /* fdlibm kernel function */ 259e02846ceSDavid Schultz int __ieee754_rem_pio2(double,double*); 26069160b1eSDavid E. O'Brien double __kernel_sin(double,double,int); 26169160b1eSDavid E. O'Brien double __kernel_cos(double,double); 26269160b1eSDavid E. O'Brien double __kernel_tan(double,double,int); 26369160b1eSDavid E. O'Brien int __kernel_rem_pio2(double*,double*,int,int,int,const int*); 2643a8617a8SJordan K. Hubbard 2653a8617a8SJordan K. Hubbard /* float versions of fdlibm kernel functions */ 266e02846ceSDavid Schultz int __ieee754_rem_pio2f(float,float*); 26759aad933SBruce Evans float __kernel_sindf(double); 26859aad933SBruce Evans float __kernel_cosdf(double); 26994a5f9beSBruce Evans float __kernel_tandf(double,int); 27069160b1eSDavid E. O'Brien int __kernel_rem_pio2f(float*,float*,int,int,int,const int*); 2713a8617a8SJordan K. Hubbard 272ef1ee63eSAlexey Zelkin #endif /* !_MATH_PRIVATE_H_ */ 273