xref: /freebsd/lib/msun/src/math_private.h (revision db612abe8df3355d1eb23bb3b50fdd97bc21e979)
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 #if BYTE_ORDER == BIG_ENDIAN
42 
43 typedef union
44 {
45   double value;
46   struct
47   {
48     u_int32_t msw;
49     u_int32_t lsw;
50   } parts;
51 } ieee_double_shape_type;
52 
53 #endif
54 
55 #if BYTE_ORDER == LITTLE_ENDIAN
56 
57 typedef union
58 {
59   double value;
60   struct
61   {
62     u_int32_t lsw;
63     u_int32_t msw;
64   } parts;
65 } ieee_double_shape_type;
66 
67 #endif
68 
69 /* Get two 32 bit ints from a double.  */
70 
71 #define EXTRACT_WORDS(ix0,ix1,d)				\
72 do {								\
73   ieee_double_shape_type ew_u;					\
74   ew_u.value = (d);						\
75   (ix0) = ew_u.parts.msw;					\
76   (ix1) = ew_u.parts.lsw;					\
77 } while (0)
78 
79 /* Get the more significant 32 bit int from a double.  */
80 
81 #define GET_HIGH_WORD(i,d)					\
82 do {								\
83   ieee_double_shape_type gh_u;					\
84   gh_u.value = (d);						\
85   (i) = gh_u.parts.msw;						\
86 } while (0)
87 
88 /* Get the less significant 32 bit int from a double.  */
89 
90 #define GET_LOW_WORD(i,d)					\
91 do {								\
92   ieee_double_shape_type gl_u;					\
93   gl_u.value = (d);						\
94   (i) = gl_u.parts.lsw;						\
95 } while (0)
96 
97 /* Set a double from two 32 bit ints.  */
98 
99 #define INSERT_WORDS(d,ix0,ix1)					\
100 do {								\
101   ieee_double_shape_type iw_u;					\
102   iw_u.parts.msw = (ix0);					\
103   iw_u.parts.lsw = (ix1);					\
104   (d) = iw_u.value;						\
105 } while (0)
106 
107 /* Set the more significant 32 bits of a double from an int.  */
108 
109 #define SET_HIGH_WORD(d,v)					\
110 do {								\
111   ieee_double_shape_type sh_u;					\
112   sh_u.value = (d);						\
113   sh_u.parts.msw = (v);						\
114   (d) = sh_u.value;						\
115 } while (0)
116 
117 /* Set the less significant 32 bits of a double from an int.  */
118 
119 #define SET_LOW_WORD(d,v)					\
120 do {								\
121   ieee_double_shape_type sl_u;					\
122   sl_u.value = (d);						\
123   sl_u.parts.lsw = (v);						\
124   (d) = sl_u.value;						\
125 } while (0)
126 
127 /*
128  * A union which permits us to convert between a float and a 32 bit
129  * int.
130  */
131 
132 typedef union
133 {
134   float value;
135   /* FIXME: Assumes 32 bit int.  */
136   unsigned int word;
137 } ieee_float_shape_type;
138 
139 /* Get a 32 bit int from a float.  */
140 
141 #define GET_FLOAT_WORD(i,d)					\
142 do {								\
143   ieee_float_shape_type gf_u;					\
144   gf_u.value = (d);						\
145   (i) = gf_u.word;						\
146 } while (0)
147 
148 /* Set a float from a 32 bit int.  */
149 
150 #define SET_FLOAT_WORD(d,i)					\
151 do {								\
152   ieee_float_shape_type sf_u;					\
153   sf_u.word = (i);						\
154   (d) = sf_u.value;						\
155 } while (0)
156 
157 #ifdef FLT_EVAL_METHOD
158 /*
159  * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
160  */
161 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
162 #define	STRICT_ASSIGN(type, lval, rval)	((lval) = (rval))
163 #else
164 #define	STRICT_ASSIGN(type, lval, rval) do {	\
165 	volatile type __lval;			\
166 						\
167 	if (sizeof(type) >= sizeof(double))	\
168 		(lval) = (rval);		\
169 	else {					\
170 		__lval = (rval);		\
171 		(lval) = __lval;		\
172 	}					\
173 } while (0)
174 #endif
175 #endif
176 
177 /*
178  * Common routine to process the arguments to nan(), nanf(), and nanl().
179  */
180 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
181 
182 #ifdef _COMPLEX_H
183 /*
184  * Inline functions that can be used to construct complex values.
185  *
186  * The C99 standard intends x+I*y to be used for this, but x+I*y is
187  * currently unusable in general since gcc introduces many overflow,
188  * underflow, sign and efficiency bugs by rewriting I*y as
189  * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
190  * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
191  * to -0.0+I*0.0.
192  */
193 static __inline float complex
194 cpackf(float x, float y)
195 {
196 	float complex z;
197 
198 	__real__ z = x;
199 	__imag__ z = y;
200 	return (z);
201 }
202 
203 static __inline double complex
204 cpack(double x, double y)
205 {
206 	double complex z;
207 
208 	__real__ z = x;
209 	__imag__ z = y;
210 	return (z);
211 }
212 
213 static __inline long double complex
214 cpackl(long double x, long double y)
215 {
216 	long double complex z;
217 
218 	__real__ z = x;
219 	__imag__ z = y;
220 	return (z);
221 }
222 #endif /* _COMPLEX_H */
223 
224 #ifdef __GNUCLIKE_ASM
225 
226 /* Asm versions of some functions. */
227 
228 #ifdef __amd64__
229 static __inline int
230 irint(double x)
231 {
232 	int n;
233 
234 	asm("cvtsd2si %1,%0" : "=r" (n) : "Y" (x));
235 	return (n);
236 }
237 #define	HAVE_EFFICIENT_IRINT
238 #endif
239 
240 #ifdef __i386__
241 static __inline int
242 irint(double x)
243 {
244 	int n;
245 
246 	asm("fistl %0" : "=m" (n) : "t" (x));
247 	return (n);
248 }
249 #define	HAVE_EFFICIENT_IRINT
250 #endif
251 
252 #endif /* __GNUCLIKE_ASM */
253 
254 /*
255  * ieee style elementary functions
256  *
257  * We rename functions here to improve other sources' diffability
258  * against fdlibm.
259  */
260 #define	__ieee754_sqrt	sqrt
261 #define	__ieee754_acos	acos
262 #define	__ieee754_acosh	acosh
263 #define	__ieee754_log	log
264 #define	__ieee754_atanh	atanh
265 #define	__ieee754_asin	asin
266 #define	__ieee754_atan2	atan2
267 #define	__ieee754_exp	exp
268 #define	__ieee754_cosh	cosh
269 #define	__ieee754_fmod	fmod
270 #define	__ieee754_pow	pow
271 #define	__ieee754_lgamma lgamma
272 #define	__ieee754_gamma	gamma
273 #define	__ieee754_lgamma_r lgamma_r
274 #define	__ieee754_gamma_r gamma_r
275 #define	__ieee754_log10	log10
276 #define	__ieee754_sinh	sinh
277 #define	__ieee754_hypot	hypot
278 #define	__ieee754_j0	j0
279 #define	__ieee754_j1	j1
280 #define	__ieee754_y0	y0
281 #define	__ieee754_y1	y1
282 #define	__ieee754_jn	jn
283 #define	__ieee754_yn	yn
284 #define	__ieee754_remainder remainder
285 #define	__ieee754_scalb	scalb
286 #define	__ieee754_sqrtf	sqrtf
287 #define	__ieee754_acosf	acosf
288 #define	__ieee754_acoshf acoshf
289 #define	__ieee754_logf	logf
290 #define	__ieee754_atanhf atanhf
291 #define	__ieee754_asinf	asinf
292 #define	__ieee754_atan2f atan2f
293 #define	__ieee754_expf	expf
294 #define	__ieee754_coshf	coshf
295 #define	__ieee754_fmodf	fmodf
296 #define	__ieee754_powf	powf
297 #define	__ieee754_lgammaf lgammaf
298 #define	__ieee754_gammaf gammaf
299 #define	__ieee754_lgammaf_r lgammaf_r
300 #define	__ieee754_gammaf_r gammaf_r
301 #define	__ieee754_log10f log10f
302 #define	__ieee754_sinhf	sinhf
303 #define	__ieee754_hypotf hypotf
304 #define	__ieee754_j0f	j0f
305 #define	__ieee754_j1f	j1f
306 #define	__ieee754_y0f	y0f
307 #define	__ieee754_y1f	y1f
308 #define	__ieee754_jnf	jnf
309 #define	__ieee754_ynf	ynf
310 #define	__ieee754_remainderf remainderf
311 #define	__ieee754_scalbf scalbf
312 
313 /* fdlibm kernel function */
314 int	__kernel_rem_pio2(double*,double*,int,int,int);
315 
316 /* double precision kernel functions */
317 int	__ieee754_rem_pio2(double,double*);
318 double	__kernel_sin(double,double,int);
319 double	__kernel_cos(double,double);
320 double	__kernel_tan(double,double,int);
321 
322 /* float precision kernel functions */
323 int	__ieee754_rem_pio2f(float,double*);
324 float	__kernel_sindf(double);
325 float	__kernel_cosdf(double);
326 float	__kernel_tandf(double,int);
327 
328 /* long double precision kernel functions */
329 long double __kernel_sinl(long double, long double, int);
330 long double __kernel_cosl(long double, long double);
331 long double __kernel_tanl(long double, long double, int);
332 
333 #endif /* !_MATH_PRIVATE_H_ */
334