Home | History | Annotate | Download | only in src
      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__) || defined(__ARM_EABI__)
     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   struct
     62   {
     63     u_int64_t w;
     64   } xparts;
     65 } ieee_double_shape_type;
     66 
     67 #endif
     68 
     69 #if IEEE_WORD_ORDER == LITTLE_ENDIAN
     70 
     71 typedef union
     72 {
     73   double value;
     74   struct
     75   {
     76     u_int32_t lsw;
     77     u_int32_t msw;
     78   } parts;
     79   struct
     80   {
     81     u_int64_t w;
     82   } xparts;
     83 } ieee_double_shape_type;
     84 
     85 #endif
     86 
     87 /* Get two 32 bit ints from a double.  */
     88 
     89 #define EXTRACT_WORDS(ix0,ix1,d)				\
     90 do {								\
     91   ieee_double_shape_type ew_u;					\
     92   ew_u.value = (d);						\
     93   (ix0) = ew_u.parts.msw;					\
     94   (ix1) = ew_u.parts.lsw;					\
     95 } while (0)
     96 
     97 /* Get a 64-bit int from a double. */
     98 #define EXTRACT_WORD64(ix,d)					\
     99 do {								\
    100   ieee_double_shape_type ew_u;					\
    101   ew_u.value = (d);						\
    102   (ix) = ew_u.xparts.w;						\
    103 } while (0)
    104 
    105 /* Get the more significant 32 bit int from a double.  */
    106 
    107 #define GET_HIGH_WORD(i,d)					\
    108 do {								\
    109   ieee_double_shape_type gh_u;					\
    110   gh_u.value = (d);						\
    111   (i) = gh_u.parts.msw;						\
    112 } while (0)
    113 
    114 /* Get the less significant 32 bit int from a double.  */
    115 
    116 #define GET_LOW_WORD(i,d)					\
    117 do {								\
    118   ieee_double_shape_type gl_u;					\
    119   gl_u.value = (d);						\
    120   (i) = gl_u.parts.lsw;						\
    121 } while (0)
    122 
    123 /* Set a double from two 32 bit ints.  */
    124 
    125 #define INSERT_WORDS(d,ix0,ix1)					\
    126 do {								\
    127   ieee_double_shape_type iw_u;					\
    128   iw_u.parts.msw = (ix0);					\
    129   iw_u.parts.lsw = (ix1);					\
    130   (d) = iw_u.value;						\
    131 } while (0)
    132 
    133 /* Set a double from a 64-bit int. */
    134 #define INSERT_WORD64(d,ix)					\
    135 do {								\
    136   ieee_double_shape_type iw_u;					\
    137   iw_u.xparts.w = (ix);						\
    138   (d) = iw_u.value;						\
    139 } while (0)
    140 
    141 /* Set the more significant 32 bits of a double from an int.  */
    142 
    143 #define SET_HIGH_WORD(d,v)					\
    144 do {								\
    145   ieee_double_shape_type sh_u;					\
    146   sh_u.value = (d);						\
    147   sh_u.parts.msw = (v);						\
    148   (d) = sh_u.value;						\
    149 } while (0)
    150 
    151 /* Set the less significant 32 bits of a double from an int.  */
    152 
    153 #define SET_LOW_WORD(d,v)					\
    154 do {								\
    155   ieee_double_shape_type sl_u;					\
    156   sl_u.value = (d);						\
    157   sl_u.parts.lsw = (v);						\
    158   (d) = sl_u.value;						\
    159 } while (0)
    160 
    161 /*
    162  * A union which permits us to convert between a float and a 32 bit
    163  * int.
    164  */
    165 
    166 typedef union
    167 {
    168   float value;
    169   /* FIXME: Assumes 32 bit int.  */
    170   unsigned int word;
    171 } ieee_float_shape_type;
    172 
    173 /* Get a 32 bit int from a float.  */
    174 
    175 #define GET_FLOAT_WORD(i,d)					\
    176 do {								\
    177   ieee_float_shape_type gf_u;					\
    178   gf_u.value = (d);						\
    179   (i) = gf_u.word;						\
    180 } while (0)
    181 
    182 /* Set a float from a 32 bit int.  */
    183 
    184 #define SET_FLOAT_WORD(d,i)					\
    185 do {								\
    186   ieee_float_shape_type sf_u;					\
    187   sf_u.word = (i);						\
    188   (d) = sf_u.value;						\
    189 } while (0)
    190 
    191 /*
    192  * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
    193  * double.
    194  */
    195 
    196 #define	EXTRACT_LDBL80_WORDS(ix0,ix1,d)				\
    197 do {								\
    198   union IEEEl2bits ew_u;					\
    199   ew_u.e = (d);							\
    200   (ix0) = ew_u.xbits.expsign;					\
    201   (ix1) = ew_u.xbits.man;					\
    202 } while (0)
    203 
    204 /*
    205  * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
    206  * long double.
    207  */
    208 
    209 #define	EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d)			\
    210 do {								\
    211   union IEEEl2bits ew_u;					\
    212   ew_u.e = (d);							\
    213   (ix0) = ew_u.xbits.expsign;					\
    214   (ix1) = ew_u.xbits.manh;					\
    215   (ix2) = ew_u.xbits.manl;					\
    216 } while (0)
    217 
    218 /* Get expsign as a 16 bit int from a long double.  */
    219 
    220 #define	GET_LDBL_EXPSIGN(i,d)					\
    221 do {								\
    222   union IEEEl2bits ge_u;					\
    223   ge_u.e = (d);							\
    224   (i) = ge_u.xbits.expsign;					\
    225 } while (0)
    226 
    227 /*
    228  * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
    229  * mantissa.
    230  */
    231 
    232 #define	INSERT_LDBL80_WORDS(d,ix0,ix1)				\
    233 do {								\
    234   union IEEEl2bits iw_u;					\
    235   iw_u.xbits.expsign = (ix0);					\
    236   iw_u.xbits.man = (ix1);					\
    237   (d) = iw_u.e;							\
    238 } while (0)
    239 
    240 /*
    241  * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
    242  * comprising the mantissa.
    243  */
    244 
    245 #define	INSERT_LDBL128_WORDS(d,ix0,ix1,ix2)			\
    246 do {								\
    247   union IEEEl2bits iw_u;					\
    248   iw_u.xbits.expsign = (ix0);					\
    249   iw_u.xbits.manh = (ix1);					\
    250   iw_u.xbits.manl = (ix2);					\
    251   (d) = iw_u.e;							\
    252 } while (0)
    253 
    254 /* Set expsign of a long double from a 16 bit int.  */
    255 
    256 #define	SET_LDBL_EXPSIGN(d,v)					\
    257 do {								\
    258   union IEEEl2bits se_u;					\
    259   se_u.e = (d);							\
    260   se_u.xbits.expsign = (v);					\
    261   (d) = se_u.e;							\
    262 } while (0)
    263 
    264 #ifdef __i386__
    265 /* Long double constants are broken on i386. */
    266 #define	LD80C(m, ex, v) {						\
    267 	.xbits.man = __CONCAT(m, ULL),					\
    268 	.xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0),	\
    269 }
    270 #else
    271 /* The above works on non-i386 too, but we use this to check v. */
    272 #define	LD80C(m, ex, v)	{ .e = (v), }
    273 #endif
    274 
    275 #ifdef FLT_EVAL_METHOD
    276 /*
    277  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
    278  */
    279 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
    280 #define	STRICT_ASSIGN(type, lval, rval)	((lval) = (rval))
    281 #else
    282 #define	STRICT_ASSIGN(type, lval, rval) do {	\
    283 	volatile type __lval;			\
    284 						\
    285 	if (sizeof(type) >= sizeof(long double))	\
    286 		(lval) = (rval);		\
    287 	else {					\
    288 		__lval = (rval);		\
    289 		(lval) = __lval;		\
    290 	}					\
    291 } while (0)
    292 #endif
    293 #endif /* FLT_EVAL_METHOD */
    294 
    295 /* Support switching the mode to FP_PE if necessary. */
    296 #if defined(__i386__) && !defined(NO_FPSETPREC)
    297 #define	ENTERI()				\
    298 	long double __retval;			\
    299 	fp_prec_t __oprec;			\
    300 						\
    301 	if ((__oprec = fpgetprec()) != FP_PE)	\
    302 		fpsetprec(FP_PE)
    303 #define	RETURNI(x) do {				\
    304 	__retval = (x);				\
    305 	if (__oprec != FP_PE)			\
    306 		fpsetprec(__oprec);		\
    307 	RETURNF(__retval);			\
    308 } while (0)
    309 #else
    310 #define	ENTERI(x)
    311 #define	RETURNI(x)	RETURNF(x)
    312 #endif
    313 
    314 /* Default return statement if hack*_t() is not used. */
    315 #define      RETURNF(v)      return (v)
    316 
    317 /*
    318  * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
    319  * a == 0, but is slower.
    320  */
    321 #define	_2sum(a, b) do {	\
    322 	__typeof(a) __s, __w;	\
    323 				\
    324 	__w = (a) + (b);	\
    325 	__s = __w - (a);	\
    326 	(b) = ((a) - (__w - __s)) + ((b) - __s); \
    327 	(a) = __w;		\
    328 } while (0)
    329 
    330 /*
    331  * 2sumF algorithm.
    332  *
    333  * "Normalize" the terms in the infinite-precision expression a + b for
    334  * the sum of 2 floating point values so that b is as small as possible
    335  * relative to 'a'.  (The resulting 'a' is the value of the expression in
    336  * the same precision as 'a' and the resulting b is the rounding error.)
    337  * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
    338  * exponent overflow or underflow must not occur.  This uses a Theorem of
    339  * Dekker (1971).  See Knuth (1981) 4.2.2 Theorem C.  The name "TwoSum"
    340  * is apparently due to Skewchuk (1997).
    341  *
    342  * For this to always work, assignment of a + b to 'a' must not retain any
    343  * extra precision in a + b.  This is required by C standards but broken
    344  * in many compilers.  The brokenness cannot be worked around using
    345  * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
    346  * algorithm would be destroyed by non-null strict assignments.  (The
    347  * compilers are correct to be broken -- the efficiency of all floating
    348  * point code calculations would be destroyed similarly if they forced the
    349  * conversions.)
    350  *
    351  * Fortunately, a case that works well can usually be arranged by building
    352  * any extra precision into the type of 'a' -- 'a' should have type float_t,
    353  * double_t or long double.  b's type should be no larger than 'a's type.
    354  * Callers should use these types with scopes as large as possible, to
    355  * reduce their own extra-precision and efficiciency problems.  In
    356  * particular, they shouldn't convert back and forth just to call here.
    357  */
    358 #ifdef DEBUG
    359 #define	_2sumF(a, b) do {				\
    360 	__typeof(a) __w;				\
    361 	volatile __typeof(a) __ia, __ib, __r, __vw;	\
    362 							\
    363 	__ia = (a);					\
    364 	__ib = (b);					\
    365 	assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib));	\
    366 							\
    367 	__w = (a) + (b);				\
    368 	(b) = ((a) - __w) + (b);			\
    369 	(a) = __w;					\
    370 							\
    371 	/* The next 2 assertions are weak if (a) is already long double. */ \
    372 	assert((long double)__ia + __ib == (long double)(a) + (b));	\
    373 	__vw = __ia + __ib;				\
    374 	__r = __ia - __vw;				\
    375 	__r += __ib;					\
    376 	assert(__vw == (a) && __r == (b));		\
    377 } while (0)
    378 #else /* !DEBUG */
    379 #define	_2sumF(a, b) do {	\
    380 	__typeof(a) __w;	\
    381 				\
    382 	__w = (a) + (b);	\
    383 	(b) = ((a) - __w) + (b); \
    384 	(a) = __w;		\
    385 } while (0)
    386 #endif /* DEBUG */
    387 
    388 /*
    389  * Set x += c, where x is represented in extra precision as a + b.
    390  * x must be sufficiently normalized and sufficiently larger than c,
    391  * and the result is then sufficiently normalized.
    392  *
    393  * The details of ordering are that |a| must be >= |c| (so that (a, c)
    394  * can be normalized without extra work to swap 'a' with c).  The details of
    395  * the normalization are that b must be small relative to the normalized 'a'.
    396  * Normalization of (a, c) makes the normalized c tiny relative to the
    397  * normalized a, so b remains small relative to 'a' in the result.  However,
    398  * b need not ever be tiny relative to 'a'.  For example, b might be about
    399  * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
    400  * That is usually enough, and adding c (which by normalization is about
    401  * 2**53 times smaller than a) cannot change b significantly.  However,
    402  * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
    403  * significantly relative to b.  The caller must ensure that significant
    404  * cancellation doesn't occur, either by having c of the same sign as 'a',
    405  * or by having |c| a few percent smaller than |a|.  Pre-normalization of
    406  * (a, b) may help.
    407  *
    408  * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
    409  * exercise 19).  We gain considerable efficiency by requiring the terms to
    410  * be sufficiently normalized and sufficiently increasing.
    411  */
    412 #define	_3sumF(a, b, c) do {	\
    413 	__typeof(a) __tmp;	\
    414 				\
    415 	__tmp = (c);		\
    416 	_2sumF(__tmp, (a));	\
    417 	(b) += (a);		\
    418 	(a) = __tmp;		\
    419 } while (0)
    420 
    421 /*
    422  * Common routine to process the arguments to nan(), nanf(), and nanl().
    423  */
    424 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
    425 
    426 #ifdef _COMPLEX_H
    427 
    428 /*
    429  * C99 specifies that complex numbers have the same representation as
    430  * an array of two elements, where the first element is the real part
    431  * and the second element is the imaginary part.
    432  */
    433 typedef union {
    434 	float complex f;
    435 	float a[2];
    436 } float_complex;
    437 typedef union {
    438 	double complex f;
    439 	double a[2];
    440 } double_complex;
    441 typedef union {
    442 	long double complex f;
    443 	long double a[2];
    444 } long_double_complex;
    445 #define	REALPART(z)	((z).a[0])
    446 #define	IMAGPART(z)	((z).a[1])
    447 
    448 /*
    449  * Inline functions that can be used to construct complex values.
    450  *
    451  * The C99 standard intends x+I*y to be used for this, but x+I*y is
    452  * currently unusable in general since gcc introduces many overflow,
    453  * underflow, sign and efficiency bugs by rewriting I*y as
    454  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
    455  * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
    456  * to -0.0+I*0.0.
    457  */
    458 static __inline float complex
    459 cpackf(float x, float y)
    460 {
    461 	float_complex z;
    462 
    463 	REALPART(z) = x;
    464 	IMAGPART(z) = y;
    465 	return (z.f);
    466 }
    467 
    468 static __inline double complex
    469 cpack(double x, double y)
    470 {
    471 	double_complex z;
    472 
    473 	REALPART(z) = x;
    474 	IMAGPART(z) = y;
    475 	return (z.f);
    476 }
    477 
    478 static __inline long double complex
    479 cpackl(long double x, long double y)
    480 {
    481 	long_double_complex z;
    482 
    483 	REALPART(z) = x;
    484 	IMAGPART(z) = y;
    485 	return (z.f);
    486 }
    487 #endif /* _COMPLEX_H */
    488 
    489 #ifdef __GNUCLIKE_ASM
    490 
    491 /* Asm versions of some functions. */
    492 
    493 #ifdef __amd64__
    494 static __inline int
    495 irint(double x)
    496 {
    497 	int n;
    498 
    499 	asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
    500 	return (n);
    501 }
    502 #define	HAVE_EFFICIENT_IRINT
    503 #endif
    504 
    505 #ifdef __i386__
    506 static __inline int
    507 irint(double x)
    508 {
    509 	int n;
    510 
    511 	asm("fistl %0" : "=m" (n) : "t" (x));
    512 	return (n);
    513 }
    514 #define	HAVE_EFFICIENT_IRINT
    515 #endif
    516 
    517 #if defined(__amd64__) || defined(__i386__)
    518 static __inline int
    519 irintl(long double x)
    520 {
    521 	int n;
    522 
    523 	asm("fistl %0" : "=m" (n) : "t" (x));
    524 	return (n);
    525 }
    526 #define	HAVE_EFFICIENT_IRINTL
    527 #endif
    528 
    529 #endif /* __GNUCLIKE_ASM */
    530 
    531 #ifdef DEBUG
    532 #if defined(__amd64__) || defined(__i386__)
    533 #define	breakpoint()	asm("int $3")
    534 #else
    535 #include <signal.h>
    536 
    537 #define	breakpoint()	raise(SIGTRAP)
    538 #endif
    539 #endif
    540 
    541 /* Write a pari script to test things externally. */
    542 #ifdef DOPRINT
    543 #include <stdio.h>
    544 
    545 #ifndef DOPRINT_SWIZZLE
    546 #define	DOPRINT_SWIZZLE		0
    547 #endif
    548 
    549 #ifdef DOPRINT_LD80
    550 
    551 #define	DOPRINT_START(xp) do {						\
    552 	uint64_t __lx;							\
    553 	uint16_t __hx;							\
    554 									\
    555 	/* Hack to give more-problematic args. */			\
    556 	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);				\
    557 	__lx ^= DOPRINT_SWIZZLE;					\
    558 	INSERT_LDBL80_WORDS(*xp, __hx, __lx);				\
    559 	printf("x = %.21Lg; ", (long double)*xp);			\
    560 } while (0)
    561 #define	DOPRINT_END1(v)							\
    562 	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
    563 #define	DOPRINT_END2(hi, lo)						\
    564 	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
    565 	    (long double)(hi), (long double)(lo))
    566 
    567 #elif defined(DOPRINT_D64)
    568 
    569 #define	DOPRINT_START(xp) do {						\
    570 	uint32_t __hx, __lx;						\
    571 									\
    572 	EXTRACT_WORDS(__hx, __lx, *xp);					\
    573 	__lx ^= DOPRINT_SWIZZLE;					\
    574 	INSERT_WORDS(*xp, __hx, __lx);					\
    575 	printf("x = %.21Lg; ", (long double)*xp);			\
    576 } while (0)
    577 #define	DOPRINT_END1(v)							\
    578 	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
    579 #define	DOPRINT_END2(hi, lo)						\
    580 	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
    581 	    (long double)(hi), (long double)(lo))
    582 
    583 #elif defined(DOPRINT_F32)
    584 
    585 #define	DOPRINT_START(xp) do {						\
    586 	uint32_t __hx;							\
    587 									\
    588 	GET_FLOAT_WORD(__hx, *xp);					\
    589 	__hx ^= DOPRINT_SWIZZLE;					\
    590 	SET_FLOAT_WORD(*xp, __hx);					\
    591 	printf("x = %.21Lg; ", (long double)*xp);			\
    592 } while (0)
    593 #define	DOPRINT_END1(v)							\
    594 	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
    595 #define	DOPRINT_END2(hi, lo)						\
    596 	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
    597 	    (long double)(hi), (long double)(lo))
    598 
    599 #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
    600 
    601 #ifndef DOPRINT_SWIZZLE_HIGH
    602 #define	DOPRINT_SWIZZLE_HIGH	0
    603 #endif
    604 
    605 #define	DOPRINT_START(xp) do {						\
    606 	uint64_t __lx, __llx;						\
    607 	uint16_t __hx;							\
    608 									\
    609 	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);			\
    610 	__llx ^= DOPRINT_SWIZZLE;					\
    611 	__lx ^= DOPRINT_SWIZZLE_HIGH;					\
    612 	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);			\
    613 	printf("x = %.36Lg; ", (long double)*xp);					\
    614 } while (0)
    615 #define	DOPRINT_END1(v)							\
    616 	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
    617 #define	DOPRINT_END2(hi, lo)						\
    618 	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",		\
    619 	    (long double)(hi), (long double)(lo))
    620 
    621 #endif /* DOPRINT_LD80 */
    622 
    623 #else /* !DOPRINT */
    624 #define	DOPRINT_START(xp)
    625 #define	DOPRINT_END1(v)
    626 #define	DOPRINT_END2(hi, lo)
    627 #endif /* DOPRINT */
    628 
    629 #define	RETURNP(x) do {			\
    630 	DOPRINT_END1(x);		\
    631 	RETURNF(x);			\
    632 } while (0)
    633 #define	RETURNPI(x) do {		\
    634 	DOPRINT_END1(x);		\
    635 	RETURNI(x);			\
    636 } while (0)
    637 #define	RETURN2P(x, y) do {		\
    638 	DOPRINT_END2((x), (y));		\
    639 	RETURNF((x) + (y));		\
    640 } while (0)
    641 #define	RETURN2PI(x, y) do {		\
    642 	DOPRINT_END2((x), (y));		\
    643 	RETURNI((x) + (y));		\
    644 } while (0)
    645 #ifdef STRUCT_RETURN
    646 #define	RETURNSP(rp) do {		\
    647 	if (!(rp)->lo_set)		\
    648 		RETURNP((rp)->hi);	\
    649 	RETURN2P((rp)->hi, (rp)->lo);	\
    650 } while (0)
    651 #define	RETURNSPI(rp) do {		\
    652 	if (!(rp)->lo_set)		\
    653 		RETURNPI((rp)->hi);	\
    654 	RETURN2PI((rp)->hi, (rp)->lo);	\
    655 } while (0)
    656 #endif
    657 #define	SUM2P(x, y) ({			\
    658 	const __typeof (x) __x = (x);	\
    659 	const __typeof (y) __y = (y);	\
    660 					\
    661 	DOPRINT_END2(__x, __y);		\
    662 	__x + __y;			\
    663 })
    664 
    665 /*
    666  * ieee style elementary functions
    667  *
    668  * We rename functions here to improve other sources' diffability
    669  * against fdlibm.
    670  */
    671 #define	__ieee754_sqrt	sqrt
    672 #define	__ieee754_acos	acos
    673 #define	__ieee754_acosh	acosh
    674 #define	__ieee754_log	log
    675 #define	__ieee754_log2	log2
    676 #define	__ieee754_atanh	atanh
    677 #define	__ieee754_asin	asin
    678 #define	__ieee754_atan2	atan2
    679 #define	__ieee754_exp	exp
    680 #define	__ieee754_cosh	cosh
    681 #define	__ieee754_fmod	fmod
    682 #define	__ieee754_pow	pow
    683 #define	__ieee754_lgamma lgamma
    684 #define	__ieee754_gamma	gamma
    685 #define	__ieee754_lgamma_r lgamma_r
    686 #define	__ieee754_gamma_r gamma_r
    687 #define	__ieee754_log10	log10
    688 #define	__ieee754_sinh	sinh
    689 #define	__ieee754_hypot	hypot
    690 #define	__ieee754_j0	j0
    691 #define	__ieee754_j1	j1
    692 #define	__ieee754_y0	y0
    693 #define	__ieee754_y1	y1
    694 #define	__ieee754_jn	jn
    695 #define	__ieee754_yn	yn
    696 #define	__ieee754_remainder remainder
    697 #define	__ieee754_scalb	scalb
    698 #define	__ieee754_sqrtf	sqrtf
    699 #define	__ieee754_acosf	acosf
    700 #define	__ieee754_acoshf acoshf
    701 #define	__ieee754_logf	logf
    702 #define	__ieee754_atanhf atanhf
    703 #define	__ieee754_asinf	asinf
    704 #define	__ieee754_atan2f atan2f
    705 #define	__ieee754_expf	expf
    706 #define	__ieee754_coshf	coshf
    707 #define	__ieee754_fmodf	fmodf
    708 #define	__ieee754_powf	powf
    709 #define	__ieee754_lgammaf lgammaf
    710 #define	__ieee754_gammaf gammaf
    711 #define	__ieee754_lgammaf_r lgammaf_r
    712 #define	__ieee754_gammaf_r gammaf_r
    713 #define	__ieee754_log10f log10f
    714 #define	__ieee754_log2f log2f
    715 #define	__ieee754_sinhf	sinhf
    716 #define	__ieee754_hypotf hypotf
    717 #define	__ieee754_j0f	j0f
    718 #define	__ieee754_j1f	j1f
    719 #define	__ieee754_y0f	y0f
    720 #define	__ieee754_y1f	y1f
    721 #define	__ieee754_jnf	jnf
    722 #define	__ieee754_ynf	ynf
    723 #define	__ieee754_remainderf remainderf
    724 #define	__ieee754_scalbf scalbf
    725 
    726 /* fdlibm kernel function */
    727 int	__kernel_rem_pio2(double*,double*,int,int,int);
    728 
    729 /* double precision kernel functions */
    730 #ifndef INLINE_REM_PIO2
    731 int	__ieee754_rem_pio2(double,double*);
    732 #endif
    733 double	__kernel_sin(double,double,int);
    734 double	__kernel_cos(double,double);
    735 double	__kernel_tan(double,double,int);
    736 double	__ldexp_exp(double,int);
    737 #ifdef _COMPLEX_H
    738 double complex __ldexp_cexp(double complex,int);
    739 #endif
    740 
    741 /* float precision kernel functions */
    742 #ifndef INLINE_REM_PIO2F
    743 int	__ieee754_rem_pio2f(float,double*);
    744 #endif
    745 #ifndef INLINE_KERNEL_SINDF
    746 float	__kernel_sindf(double);
    747 #endif
    748 #ifndef INLINE_KERNEL_COSDF
    749 float	__kernel_cosdf(double);
    750 #endif
    751 #ifndef INLINE_KERNEL_TANDF
    752 float	__kernel_tandf(double,int);
    753 #endif
    754 float	__ldexp_expf(float,int);
    755 #ifdef _COMPLEX_H
    756 float complex __ldexp_cexpf(float complex,int);
    757 #endif
    758 
    759 /* long double precision kernel functions */
    760 long double __kernel_sinl(long double, long double, int);
    761 long double __kernel_cosl(long double, long double);
    762 long double __kernel_tanl(long double, long double, int);
    763 
    764 #endif /* !_MATH_PRIVATE_H_ */
    765