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: stable/9/lib/msun/src/math_private.h 239529 2012-08-21 19:45:48Z dim $
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 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 #ifdef FLT_EVAL_METHOD
192 /*
193 * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
194 */
195 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
196 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
197 #else
198 #define STRICT_ASSIGN(type, lval, rval) do { \
199 volatile type __lval; \
200 \
201 if (sizeof(type) >= sizeof(double)) \
202 (lval) = (rval); \
203 else { \
204 __lval = (rval); \
205 (lval) = __lval; \
206 } \
207 } while (0)
208 #endif
209 #endif
210
211 /*
212 * Common routine to process the arguments to nan(), nanf(), and nanl().
213 */
214 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
215
216 #ifdef _COMPLEX_H
217
218 /*
219 * C99 specifies that complex numbers have the same representation as
220 * an array of two elements, where the first element is the real part
221 * and the second element is the imaginary part.
222 */
223 typedef union {
224 float complex f;
225 float a[2];
226 } float_complex;
227 typedef union {
228 double complex f;
229 double a[2];
230 } double_complex;
231 typedef union {
232 long double complex f;
233 long double a[2];
234 } long_double_complex;
235 #define REALPART(z) ((z).a[0])
236 #define IMAGPART(z) ((z).a[1])
237
238 /*
239 * Inline functions that can be used to construct complex values.
240 *
241 * The C99 standard intends x+I*y to be used for this, but x+I*y is
242 * currently unusable in general since gcc introduces many overflow,
243 * underflow, sign and efficiency bugs by rewriting I*y as
244 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
245 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
246 * to -0.0+I*0.0.
247 */
248 static __inline float complex
cpackf(float x,float y)249 cpackf(float x, float y)
250 {
251 float_complex z;
252
253 REALPART(z) = x;
254 IMAGPART(z) = y;
255 return (z.f);
256 }
257
258 static __inline double complex
cpack(double x,double y)259 cpack(double x, double y)
260 {
261 double_complex z;
262
263 REALPART(z) = x;
264 IMAGPART(z) = y;
265 return (z.f);
266 }
267
268 static __inline long double complex
cpackl(long double x,long double y)269 cpackl(long double x, long double y)
270 {
271 long_double_complex z;
272
273 REALPART(z) = x;
274 IMAGPART(z) = y;
275 return (z.f);
276 }
277 #endif /* _COMPLEX_H */
278
279 #ifdef __GNUCLIKE_ASM
280
281 /* Asm versions of some functions. */
282
283 #ifdef __amd64__
284 static __inline int
irint(double x)285 irint(double x)
286 {
287 int n;
288
289 asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
290 return (n);
291 }
292 #define HAVE_EFFICIENT_IRINT
293 #endif
294
295 #ifdef __i386__
296 static __inline int
irint(double x)297 irint(double x)
298 {
299 int n;
300
301 asm("fistl %0" : "=m" (n) : "t" (x));
302 return (n);
303 }
304 #define HAVE_EFFICIENT_IRINT
305 #endif
306
307 #endif /* __GNUCLIKE_ASM */
308
309 /*
310 * ieee style elementary functions
311 *
312 * We rename functions here to improve other sources' diffability
313 * against fdlibm.
314 */
315 #define __ieee754_sqrt sqrt
316 #define __ieee754_acos acos
317 #define __ieee754_acosh acosh
318 #define __ieee754_log log
319 #define __ieee754_log2 log2
320 #define __ieee754_atanh atanh
321 #define __ieee754_asin asin
322 #define __ieee754_atan2 atan2
323 #define __ieee754_exp exp
324 #define __ieee754_cosh cosh
325 #define __ieee754_fmod fmod
326 #define __ieee754_pow pow
327 #define __ieee754_lgamma lgamma
328 #define __ieee754_gamma gamma
329 #define __ieee754_lgamma_r lgamma_r
330 #define __ieee754_gamma_r gamma_r
331 #define __ieee754_log10 log10
332 #define __ieee754_sinh sinh
333 #define __ieee754_hypot hypot
334 #define __ieee754_j0 j0
335 #define __ieee754_j1 j1
336 #define __ieee754_y0 y0
337 #define __ieee754_y1 y1
338 #define __ieee754_jn jn
339 #define __ieee754_yn yn
340 #define __ieee754_remainder remainder
341 #define __ieee754_scalb scalb
342 #define __ieee754_sqrtf sqrtf
343 #define __ieee754_acosf acosf
344 #define __ieee754_acoshf acoshf
345 #define __ieee754_logf logf
346 #define __ieee754_atanhf atanhf
347 #define __ieee754_asinf asinf
348 #define __ieee754_atan2f atan2f
349 #define __ieee754_expf expf
350 #define __ieee754_coshf coshf
351 #define __ieee754_fmodf fmodf
352 #define __ieee754_powf powf
353 #define __ieee754_lgammaf lgammaf
354 #define __ieee754_gammaf gammaf
355 #define __ieee754_lgammaf_r lgammaf_r
356 #define __ieee754_gammaf_r gammaf_r
357 #define __ieee754_log10f log10f
358 #define __ieee754_log2f log2f
359 #define __ieee754_sinhf sinhf
360 #define __ieee754_hypotf hypotf
361 #define __ieee754_j0f j0f
362 #define __ieee754_j1f j1f
363 #define __ieee754_y0f y0f
364 #define __ieee754_y1f y1f
365 #define __ieee754_jnf jnf
366 #define __ieee754_ynf ynf
367 #define __ieee754_remainderf remainderf
368 #define __ieee754_scalbf scalbf
369
370 /* fdlibm kernel function */
371 int __kernel_rem_pio2(double*,double*,int,int,int);
372
373 /* double precision kernel functions */
374 #ifndef INLINE_REM_PIO2
375 int __ieee754_rem_pio2(double,double*);
376 #endif
377 double __kernel_sin(double,double,int);
378 double __kernel_cos(double,double);
379 double __kernel_tan(double,double,int);
380
381 /* float precision kernel functions */
382 #ifndef INLINE_REM_PIO2F
383 int __ieee754_rem_pio2f(float,double*);
384 #endif
385 #ifndef INLINE_KERNEL_SINDF
386 float __kernel_sindf(double);
387 #endif
388 #ifndef INLINE_KERNEL_COSDF
389 float __kernel_cosdf(double);
390 #endif
391 #ifndef INLINE_KERNEL_TANDF
392 float __kernel_tandf(double,int);
393 #endif
394
395 /* long double precision kernel functions */
396 long double __kernel_sinl(long double, long double, int);
397 long double __kernel_cosl(long double, long double);
398 long double __kernel_tanl(long double, long double, int);
399
400 #endif /* !_MATH_PRIVATE_H_ */
401