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