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 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 /* Get expsign as a 16 bit int from a long double. */
192
193 #define GET_LDBL_EXPSIGN(i,d) \
194 do { \
195 union IEEEl2bits ge_u; \
196 ge_u.e = (d); \
197 (i) = ge_u.xbits.expsign; \
198 } while (0)
199
200 /* Set expsign of a long double from a 16 bit int. */
201
202 #define SET_LDBL_EXPSIGN(d,v) \
203 do { \
204 union IEEEl2bits se_u; \
205 se_u.e = (d); \
206 se_u.xbits.expsign = (v); \
207 (d) = se_u.e; \
208 } while (0)
209
210 #ifdef __i386__
211 /* Long double constants are broken on i386. */
212 #define LD80C(m, ex, v) { \
213 .xbits.man = __CONCAT(m, ULL), \
214 .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \
215 }
216 #else
217 /* The above works on non-i386 too, but we use this to check v. */
218 #define LD80C(m, ex, v) { .e = (v), }
219 #endif
220
221 #ifdef FLT_EVAL_METHOD
222 /*
223 * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
224 */
225 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
226 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
227 #else
228 #define STRICT_ASSIGN(type, lval, rval) do { \
229 volatile type __lval; \
230 \
231 if (sizeof(type) >= sizeof(long double)) \
232 (lval) = (rval); \
233 else { \
234 __lval = (rval); \
235 (lval) = __lval; \
236 } \
237 } while (0)
238 #endif
239 #endif /* FLT_EVAL_METHOD */
240
241 /* Support switching the mode to FP_PE if necessary. */
242 #if defined(__i386__) && !defined(NO_FPSETPREC)
243 #define ENTERI() \
244 long double __retval; \
245 fp_prec_t __oprec; \
246 \
247 if ((__oprec = fpgetprec()) != FP_PE) \
248 fpsetprec(FP_PE)
249 #define RETURNI(x) do { \
250 __retval = (x); \
251 if (__oprec != FP_PE) \
252 fpsetprec(__oprec); \
253 RETURNF(__retval); \
254 } while (0)
255 #else
256 #define ENTERI(x)
257 #define RETURNI(x) RETURNF(x)
258 #endif
259
260 /* Default return statement if hack*_t() is not used. */
261 #define RETURNF(v) return (v)
262
263 /*
264 * Common routine to process the arguments to nan(), nanf(), and nanl().
265 */
266 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
267
268 #ifdef _COMPLEX_H
269
270 /*
271 * C99 specifies that complex numbers have the same representation as
272 * an array of two elements, where the first element is the real part
273 * and the second element is the imaginary part.
274 */
275 typedef union {
276 float complex f;
277 float a[2];
278 } float_complex;
279 typedef union {
280 double complex f;
281 double a[2];
282 } double_complex;
283 typedef union {
284 long double complex f;
285 long double a[2];
286 } long_double_complex;
287 #define REALPART(z) ((z).a[0])
288 #define IMAGPART(z) ((z).a[1])
289
290 /*
291 * Inline functions that can be used to construct complex values.
292 *
293 * The C99 standard intends x+I*y to be used for this, but x+I*y is
294 * currently unusable in general since gcc introduces many overflow,
295 * underflow, sign and efficiency bugs by rewriting I*y as
296 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
297 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
298 * to -0.0+I*0.0.
299 */
300 static __inline float complex
cpackf(float x,float y)301 cpackf(float x, float y)
302 {
303 float_complex z;
304
305 REALPART(z) = x;
306 IMAGPART(z) = y;
307 return (z.f);
308 }
309
310 static __inline double complex
cpack(double x,double y)311 cpack(double x, double y)
312 {
313 double_complex z;
314
315 REALPART(z) = x;
316 IMAGPART(z) = y;
317 return (z.f);
318 }
319
320 static __inline long double complex
cpackl(long double x,long double y)321 cpackl(long double x, long double y)
322 {
323 long_double_complex z;
324
325 REALPART(z) = x;
326 IMAGPART(z) = y;
327 return (z.f);
328 }
329 #endif /* _COMPLEX_H */
330
331 #ifdef __GNUCLIKE_ASM
332
333 /* Asm versions of some functions. */
334
335 #ifdef __amd64__
336 static __inline int
irint(double x)337 irint(double x)
338 {
339 int n;
340
341 asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x));
342 return (n);
343 }
344 #define HAVE_EFFICIENT_IRINT
345 #endif
346
347 #ifdef __i386__
348 static __inline int
irint(double x)349 irint(double x)
350 {
351 int n;
352
353 asm("fistl %0" : "=m" (n) : "t" (x));
354 return (n);
355 }
356 #define HAVE_EFFICIENT_IRINT
357 #endif
358
359 #if defined(__amd64__) || defined(__i386__)
360 static __inline int
irintl(long double x)361 irintl(long double x)
362 {
363 int n;
364
365 asm("fistl %0" : "=m" (n) : "t" (x));
366 return (n);
367 }
368 #define HAVE_EFFICIENT_IRINTL
369 #endif
370
371 #endif /* __GNUCLIKE_ASM */
372
373 /*
374 * ieee style elementary functions
375 *
376 * We rename functions here to improve other sources' diffability
377 * against fdlibm.
378 */
379 #define __ieee754_sqrt sqrt
380 #define __ieee754_acos acos
381 #define __ieee754_acosh acosh
382 #define __ieee754_log log
383 #define __ieee754_log2 log2
384 #define __ieee754_atanh atanh
385 #define __ieee754_asin asin
386 #define __ieee754_atan2 atan2
387 #define __ieee754_exp exp
388 #define __ieee754_cosh cosh
389 #define __ieee754_fmod fmod
390 #define __ieee754_pow pow
391 #define __ieee754_lgamma lgamma
392 #define __ieee754_gamma gamma
393 #define __ieee754_lgamma_r lgamma_r
394 #define __ieee754_gamma_r gamma_r
395 #define __ieee754_log10 log10
396 #define __ieee754_sinh sinh
397 #define __ieee754_hypot hypot
398 #define __ieee754_j0 j0
399 #define __ieee754_j1 j1
400 #define __ieee754_y0 y0
401 #define __ieee754_y1 y1
402 #define __ieee754_jn jn
403 #define __ieee754_yn yn
404 #define __ieee754_remainder remainder
405 #define __ieee754_scalb scalb
406 #define __ieee754_sqrtf sqrtf
407 #define __ieee754_acosf acosf
408 #define __ieee754_acoshf acoshf
409 #define __ieee754_logf logf
410 #define __ieee754_atanhf atanhf
411 #define __ieee754_asinf asinf
412 #define __ieee754_atan2f atan2f
413 #define __ieee754_expf expf
414 #define __ieee754_coshf coshf
415 #define __ieee754_fmodf fmodf
416 #define __ieee754_powf powf
417 #define __ieee754_lgammaf lgammaf
418 #define __ieee754_gammaf gammaf
419 #define __ieee754_lgammaf_r lgammaf_r
420 #define __ieee754_gammaf_r gammaf_r
421 #define __ieee754_log10f log10f
422 #define __ieee754_log2f log2f
423 #define __ieee754_sinhf sinhf
424 #define __ieee754_hypotf hypotf
425 #define __ieee754_j0f j0f
426 #define __ieee754_j1f j1f
427 #define __ieee754_y0f y0f
428 #define __ieee754_y1f y1f
429 #define __ieee754_jnf jnf
430 #define __ieee754_ynf ynf
431 #define __ieee754_remainderf remainderf
432 #define __ieee754_scalbf scalbf
433
434 /* fdlibm kernel function */
435 int __kernel_rem_pio2(double*,double*,int,int,int);
436
437 /* double precision kernel functions */
438 #ifndef INLINE_REM_PIO2
439 int __ieee754_rem_pio2(double,double*);
440 #endif
441 double __kernel_sin(double,double,int);
442 double __kernel_cos(double,double);
443 double __kernel_tan(double,double,int);
444 double __ldexp_exp(double,int);
445 #ifdef _COMPLEX_H
446 double complex __ldexp_cexp(double complex,int);
447 #endif
448
449 /* float precision kernel functions */
450 #ifndef INLINE_REM_PIO2F
451 int __ieee754_rem_pio2f(float,double*);
452 #endif
453 #ifndef INLINE_KERNEL_SINDF
454 float __kernel_sindf(double);
455 #endif
456 #ifndef INLINE_KERNEL_COSDF
457 float __kernel_cosdf(double);
458 #endif
459 #ifndef INLINE_KERNEL_TANDF
460 float __kernel_tandf(double,int);
461 #endif
462 float __ldexp_expf(float,int);
463 #ifdef _COMPLEX_H
464 float complex __ldexp_cexpf(float complex,int);
465 #endif
466
467 /* long double precision kernel functions */
468 long double __kernel_sinl(long double, long double, int);
469 long double __kernel_cosl(long double, long double);
470 long double __kernel_tanl(long double, long double, int);
471
472 #endif /* !_MATH_PRIVATE_H_ */
473