xref: /freebsd/lib/msun/src/math_private.h (revision 59aad933ab2d68f9c1b68caf2344f1e81f0a13f4)
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