2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
9 * ====================================================
13 * from: @(#)fdlibm.h 5.1 93/09/24
14 * $FreeBSD: head/lib/msun/src/math_private.h 241051 2012-09-29 16:40:12Z kargl $
17 #ifndef _MATH_PRIVATE_H_
18 #define _MATH_PRIVATE_H_
20 #include <sys/types.h>
21 #include <machine/endian.h>
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.
37 * A union which permits us to convert between a double and two 32 bit
41 #define IEEE_WORD_ORDER BYTE_ORDER
43 #if IEEE_WORD_ORDER == BIG_ENDIAN
57 } ieee_double_shape_type;
61 #if IEEE_WORD_ORDER == LITTLE_ENDIAN
75 } ieee_double_shape_type;
79 /* Get two 32 bit ints from a double. */
81 #define EXTRACT_WORDS(ix0,ix1,d) \
83 ieee_double_shape_type ew_u; \
85 (ix0) = ew_u.parts.msw; \
86 (ix1) = ew_u.parts.lsw; \
89 /* Get a 64-bit int from a double. */
90 #define EXTRACT_WORD64(ix,d) \
92 ieee_double_shape_type ew_u; \
94 (ix) = ew_u.xparts.w; \
97 /* Get the more significant 32 bit int from a double. */
99 #define GET_HIGH_WORD(i,d) \
101 ieee_double_shape_type gh_u; \
103 (i) = gh_u.parts.msw; \
106 /* Get the less significant 32 bit int from a double. */
108 #define GET_LOW_WORD(i,d) \
110 ieee_double_shape_type gl_u; \
112 (i) = gl_u.parts.lsw; \
115 /* Set a double from two 32 bit ints. */
117 #define INSERT_WORDS(d,ix0,ix1) \
119 ieee_double_shape_type iw_u; \
120 iw_u.parts.msw = (ix0); \
121 iw_u.parts.lsw = (ix1); \
125 /* Set a double from a 64-bit int. */
126 #define INSERT_WORD64(d,ix) \
128 ieee_double_shape_type iw_u; \
129 iw_u.xparts.w = (ix); \
133 /* Set the more significant 32 bits of a double from an int. */
135 #define SET_HIGH_WORD(d,v) \
137 ieee_double_shape_type sh_u; \
139 sh_u.parts.msw = (v); \
143 /* Set the less significant 32 bits of a double from an int. */
145 #define SET_LOW_WORD(d,v) \
147 ieee_double_shape_type sl_u; \
149 sl_u.parts.lsw = (v); \
154 * A union which permits us to convert between a float and a 32 bit
161 /* FIXME: Assumes 32 bit int. */
163 } ieee_float_shape_type;
165 /* Get a 32 bit int from a float. */
167 #define GET_FLOAT_WORD(i,d) \
169 ieee_float_shape_type gf_u; \
174 /* Set a float from a 32 bit int. */
176 #define SET_FLOAT_WORD(d,i) \
178 ieee_float_shape_type sf_u; \
183 /* Get expsign as a 16 bit int from a long double. */
185 #define GET_LDBL_EXPSIGN(i,d) \
187 union IEEEl2bits ge_u; \
189 (i) = ge_u.xbits.expsign; \
192 /* Set expsign of a long double from a 16 bit int. */
194 #define SET_LDBL_EXPSIGN(d,v) \
196 union IEEEl2bits se_u; \
198 se_u.xbits.expsign = (v); \
203 /* Long double constants are broken on i386. */
204 #define LD80C(m, ex, v) { \
205 .xbits.man = __CONCAT(m, ULL), \
206 .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \
209 /* The above works on non-i386 too, but we use this to check v. */
210 #define LD80C(m, ex, v) { .e = (v), }
213 #ifdef FLT_EVAL_METHOD
215 * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
217 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
218 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
220 #define STRICT_ASSIGN(type, lval, rval) do { \
221 volatile type __lval; \
223 if (sizeof(type) >= sizeof(long double)) \
231 #endif /* FLT_EVAL_METHOD */
233 /* Support switching the mode to FP_PE if necessary. */
234 #if defined(__i386__) && !defined(NO_FPSETPREC)
236 long double __retval; \
239 if ((__oprec = fpgetprec()) != FP_PE) \
241 #define RETURNI(x) do { \
243 if (__oprec != FP_PE) \
244 fpsetprec(__oprec); \
249 #define RETURNI(x) RETURNF(x)
252 /* Default return statement if hack*_t() is not used. */
253 #define RETURNF(v) return (v)
256 * Common routine to process the arguments to nan(), nanf(), and nanl().
258 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
263 * C99 specifies that complex numbers have the same representation as
264 * an array of two elements, where the first element is the real part
265 * and the second element is the imaginary part.
276 long double complex f;
278 } long_double_complex;
279 #define REALPART(z) ((z).a[0])
280 #define IMAGPART(z) ((z).a[1])
283 * Inline functions that can be used to construct complex values.
285 * The C99 standard intends x+I*y to be used for this, but x+I*y is
286 * currently unusable in general since gcc introduces many overflow,
287 * underflow, sign and efficiency bugs by rewriting I*y as
288 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
289 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
292 static __inline float complex
293 cpackf(float x, float y)
302 static __inline double complex
303 cpack(double x, double y)
312 static __inline long double complex
313 cpackl(long double x, long double y)
315 long_double_complex z;
321 #endif /* _COMPLEX_H */
325 /* Asm versions of some functions. */
333 asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
336 #define HAVE_EFFICIENT_IRINT
345 asm("fistl %0" : "=m" (n) : "t" (x));
348 #define HAVE_EFFICIENT_IRINT
352 irintl(long double x)
356 asm("fistl %0" : "=m" (n) : "t" (x));
359 #define HAVE_EFFICIENT_IRINTL
361 #endif /* __GNUC__ */
364 * ieee style elementary functions
366 * We rename functions here to improve other sources' diffability
369 #define __ieee754_sqrt sqrt
370 #define __ieee754_acos acos
371 #define __ieee754_acosh acosh
372 #define __ieee754_log log
373 #define __ieee754_log2 log2
374 #define __ieee754_atanh atanh
375 #define __ieee754_asin asin
376 #define __ieee754_atan2 atan2
377 #define __ieee754_exp exp
378 #define __ieee754_cosh cosh
379 #define __ieee754_fmod fmod
380 #define __ieee754_pow pow
381 #define __ieee754_lgamma lgamma
382 #define __ieee754_gamma gamma
383 #define __ieee754_lgamma_r lgamma_r
384 #define __ieee754_gamma_r gamma_r
385 #define __ieee754_log10 log10
386 #define __ieee754_sinh sinh
387 #define __ieee754_hypot hypot
388 #define __ieee754_j0 j0
389 #define __ieee754_j1 j1
390 #define __ieee754_y0 y0
391 #define __ieee754_y1 y1
392 #define __ieee754_jn jn
393 #define __ieee754_yn yn
394 #define __ieee754_remainder remainder
395 #define __ieee754_scalb scalb
396 #define __ieee754_sqrtf sqrtf
397 #define __ieee754_acosf acosf
398 #define __ieee754_acoshf acoshf
399 #define __ieee754_logf logf
400 #define __ieee754_atanhf atanhf
401 #define __ieee754_asinf asinf
402 #define __ieee754_atan2f atan2f
403 #define __ieee754_expf expf
404 #define __ieee754_coshf coshf
405 #define __ieee754_fmodf fmodf
406 #define __ieee754_powf powf
407 #define __ieee754_lgammaf lgammaf
408 #define __ieee754_gammaf gammaf
409 #define __ieee754_lgammaf_r lgammaf_r
410 #define __ieee754_gammaf_r gammaf_r
411 #define __ieee754_log10f log10f
412 #define __ieee754_log2f log2f
413 #define __ieee754_sinhf sinhf
414 #define __ieee754_hypotf hypotf
415 #define __ieee754_j0f j0f
416 #define __ieee754_j1f j1f
417 #define __ieee754_y0f y0f
418 #define __ieee754_y1f y1f
419 #define __ieee754_jnf jnf
420 #define __ieee754_ynf ynf
421 #define __ieee754_remainderf remainderf
422 #define __ieee754_scalbf scalbf
424 /* fdlibm kernel function */
425 int __kernel_rem_pio2(double*,double*,int,int,int);
427 /* double precision kernel functions */
428 #ifndef INLINE_REM_PIO2
429 int __ieee754_rem_pio2(double,double*);
431 double __kernel_sin(double,double,int);
432 double __kernel_cos(double,double);
433 double __kernel_tan(double,double,int);
434 double __ldexp_exp(double,int);
436 double complex __ldexp_cexp(double complex,int);
439 /* float precision kernel functions */
440 #ifndef INLINE_REM_PIO2F
441 int __ieee754_rem_pio2f(float,double*);
443 #ifndef INLINE_KERNEL_SINDF
444 float __kernel_sindf(double);
446 #ifndef INLINE_KERNEL_COSDF
447 float __kernel_cosdf(double);
449 #ifndef INLINE_KERNEL_TANDF
450 float __kernel_tandf(double,int);
452 float __ldexp_expf(float,int);
454 float complex __ldexp_cexpf(float complex,int);
457 /* long double precision kernel functions */
458 long double __kernel_sinl(long double, long double, int);
459 long double __kernel_cosl(long double, long double);
460 long double __kernel_tanl(long double, long double, int);
462 #endif /* !_MATH_PRIVATE_H_ */