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 <cfloat>
21 #include <stdint.h>
22 #include <sys/types.h>
23
24 #include "fdlibm.h"
25
26 #include "mozilla/EndianUtils.h"
27
28 /*
29 * Emulate FreeBSD internal double types.
30 * Adapted from https://github.com/freebsd/freebsd-src/search?q=__double_t
31 */
32
33 #ifdef __LP64__
34 typedef double __double_t;
35 #ifndef __float_t
36 typedef float __float_t;
37 #endif
38 #else
39 typedef long double __double_t;
40 #endif
41 typedef __double_t double_t;
42
43 /*
44 * The original fdlibm code used statements like:
45 * n0 = ((*(int*)&one)>>29)^1; * index of high word *
46 * ix0 = *(n0+(int*)&x); * high word of x *
47 * ix1 = *((1-n0)+(int*)&x); * low word of x *
48 * to dig two 32 bit words out of the 64 bit IEEE floating point
49 * value. That is non-ANSI, and, moreover, the gcc instruction
50 * scheduler gets it wrong. We instead use the following macros.
51 * Unlike the original code, we determine the endianness at compile
52 * time, not at run time; I don't see much benefit to selecting
53 * endianness at run time.
54 */
55
56 #ifndef u_int32_t
57 #define u_int32_t uint32_t
58 #endif
59 #ifndef u_int64_t
60 #define u_int64_t uint64_t
61 #endif
62
63 /* A union which permits us to convert between a long double and
64 four 32 bit ints. */
65
66 #if MOZ_BIG_ENDIAN()
67
68 typedef union
69 {
70 long double value;
71 struct {
72 u_int32_t mswhi;
73 u_int32_t mswlo;
74 u_int32_t lswhi;
75 u_int32_t lswlo;
76 } parts32;
77 struct {
78 u_int64_t msw;
79 u_int64_t lsw;
80 } parts64;
81 } ieee_quad_shape_type;
82
83 #endif
84
85 #if MOZ_LITTLE_ENDIAN()
86
87 typedef union
88 {
89 long double value;
90 struct {
91 u_int32_t lswlo;
92 u_int32_t lswhi;
93 u_int32_t mswlo;
94 u_int32_t mswhi;
95 } parts32;
96 struct {
97 u_int64_t lsw;
98 u_int64_t msw;
99 } parts64;
100 } ieee_quad_shape_type;
101
102 #endif
103
104 /*
105 * A union which permits us to convert between a double and two 32 bit
106 * ints.
107 */
108
109 #if MOZ_BIG_ENDIAN()
110
111 typedef union
112 {
113 double value;
114 struct
115 {
116 u_int32_t msw;
117 u_int32_t lsw;
118 } parts;
119 struct
120 {
121 u_int64_t w;
122 } xparts;
123 } ieee_double_shape_type;
124
125 #endif
126
127 #if MOZ_LITTLE_ENDIAN()
128
129 typedef union
130 {
131 double value;
132 struct
133 {
134 u_int32_t lsw;
135 u_int32_t msw;
136 } parts;
137 struct
138 {
139 u_int64_t w;
140 } xparts;
141 } ieee_double_shape_type;
142
143 #endif
144
145 /* Get two 32 bit ints from a double. */
146
147 #define EXTRACT_WORDS(ix0,ix1,d) \
148 do { \
149 ieee_double_shape_type ew_u; \
150 ew_u.value = (d); \
151 (ix0) = ew_u.parts.msw; \
152 (ix1) = ew_u.parts.lsw; \
153 } while (0)
154
155 /* Get a 64-bit int from a double. */
156 #define EXTRACT_WORD64(ix,d) \
157 do { \
158 ieee_double_shape_type ew_u; \
159 ew_u.value = (d); \
160 (ix) = ew_u.xparts.w; \
161 } while (0)
162
163 /* Get the more significant 32 bit int from a double. */
164
165 #define GET_HIGH_WORD(i,d) \
166 do { \
167 ieee_double_shape_type gh_u; \
168 gh_u.value = (d); \
169 (i) = gh_u.parts.msw; \
170 } while (0)
171
172 /* Get the less significant 32 bit int from a double. */
173
174 #define GET_LOW_WORD(i,d) \
175 do { \
176 ieee_double_shape_type gl_u; \
177 gl_u.value = (d); \
178 (i) = gl_u.parts.lsw; \
179 } while (0)
180
181 /* Set a double from two 32 bit ints. */
182
183 #define INSERT_WORDS(d,ix0,ix1) \
184 do { \
185 ieee_double_shape_type iw_u; \
186 iw_u.parts.msw = (ix0); \
187 iw_u.parts.lsw = (ix1); \
188 (d) = iw_u.value; \
189 } while (0)
190
191 /* Set a double from a 64-bit int. */
192 #define INSERT_WORD64(d,ix) \
193 do { \
194 ieee_double_shape_type iw_u; \
195 iw_u.xparts.w = (ix); \
196 (d) = iw_u.value; \
197 } while (0)
198
199 /* Set the more significant 32 bits of a double from an int. */
200
201 #define SET_HIGH_WORD(d,v) \
202 do { \
203 ieee_double_shape_type sh_u; \
204 sh_u.value = (d); \
205 sh_u.parts.msw = (v); \
206 (d) = sh_u.value; \
207 } while (0)
208
209 /* Set the less significant 32 bits of a double from an int. */
210
211 #define SET_LOW_WORD(d,v) \
212 do { \
213 ieee_double_shape_type sl_u; \
214 sl_u.value = (d); \
215 sl_u.parts.lsw = (v); \
216 (d) = sl_u.value; \
217 } while (0)
218
219 /*
220 * A union which permits us to convert between a float and a 32 bit
221 * int.
222 */
223
224 typedef union
225 {
226 float value;
227 /* FIXME: Assumes 32 bit int. */
228 unsigned int word;
229 } ieee_float_shape_type;
230
231 /* Get a 32 bit int from a float. */
232
233 #define GET_FLOAT_WORD(i,d) \
234 do { \
235 ieee_float_shape_type gf_u; \
236 gf_u.value = (d); \
237 (i) = gf_u.word; \
238 } while (0)
239
240 /* Set a float from a 32 bit int. */
241
242 #define SET_FLOAT_WORD(d,i) \
243 do { \
244 ieee_float_shape_type sf_u; \
245 sf_u.word = (i); \
246 (d) = sf_u.value; \
247 } while (0)
248
249 /*
250 * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
251 * double.
252 */
253
254 #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
255 do { \
256 union IEEEl2bits ew_u; \
257 ew_u.e = (d); \
258 (ix0) = ew_u.xbits.expsign; \
259 (ix1) = ew_u.xbits.man; \
260 } while (0)
261
262 /*
263 * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
264 * long double.
265 */
266
267 #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
268 do { \
269 union IEEEl2bits ew_u; \
270 ew_u.e = (d); \
271 (ix0) = ew_u.xbits.expsign; \
272 (ix1) = ew_u.xbits.manh; \
273 (ix2) = ew_u.xbits.manl; \
274 } while (0)
275
276 /* Get expsign as a 16 bit int from a long double. */
277
278 #define GET_LDBL_EXPSIGN(i,d) \
279 do { \
280 union IEEEl2bits ge_u; \
281 ge_u.e = (d); \
282 (i) = ge_u.xbits.expsign; \
283 } while (0)
284
285 /*
286 * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
287 * mantissa.
288 */
289
290 #define INSERT_LDBL80_WORDS(d,ix0,ix1) \
291 do { \
292 union IEEEl2bits iw_u; \
293 iw_u.xbits.expsign = (ix0); \
294 iw_u.xbits.man = (ix1); \
295 (d) = iw_u.e; \
296 } while (0)
297
298 /*
299 * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
300 * comprising the mantissa.
301 */
302
303 #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
304 do { \
305 union IEEEl2bits iw_u; \
306 iw_u.xbits.expsign = (ix0); \
307 iw_u.xbits.manh = (ix1); \
308 iw_u.xbits.manl = (ix2); \
309 (d) = iw_u.e; \
310 } while (0)
311
312 /* Set expsign of a long double from a 16 bit int. */
313
314 #define SET_LDBL_EXPSIGN(d,v) \
315 do { \
316 union IEEEl2bits se_u; \
317 se_u.e = (d); \
318 se_u.xbits.expsign = (v); \
319 (d) = se_u.e; \
320 } while (0)
321
322 #ifdef __i386__
323 /* Long double constants are broken on i386. */
324 #define LD80C(m, ex, v) { \
325 .xbits.man = __CONCAT(m, ULL), \
326 .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \
327 }
328 #else
329 /* The above works on non-i386 too, but we use this to check v. */
330 #define LD80C(m, ex, v) { .e = (v), }
331 #endif
332
333 #ifdef FLT_EVAL_METHOD
334 /*
335 * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
336 */
337 #if !defined(_MSC_VER) && (FLT_EVAL_METHOD == 0 || __GNUC__ == 0)
338 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
339 #else
340 #define STRICT_ASSIGN(type, lval, rval) do { \
341 volatile type __lval; \
342 \
343 if (sizeof(type) >= sizeof(long double)) \
344 (lval) = (rval); \
345 else { \
346 __lval = (rval); \
347 (lval) = __lval; \
348 } \
349 } while (0)
350 #endif
351 #else
352 #define STRICT_ASSIGN(type, lval, rval) do { \
353 volatile type __lval; \
354 \
355 if (sizeof(type) >= sizeof(long double)) \
356 (lval) = (rval); \
357 else { \
358 __lval = (rval); \
359 (lval) = __lval; \
360 } \
361 } while (0)
362 #endif /* FLT_EVAL_METHOD */
363
364 /* Support switching the mode to FP_PE if necessary. */
365 #if defined(__i386__) && !defined(NO_FPSETPREC)
366 #define ENTERI() ENTERIT(long double)
367 #define ENTERIT(returntype) \
368 returntype __retval; \
369 fp_prec_t __oprec; \
370 \
371 if ((__oprec = fpgetprec()) != FP_PE) \
372 fpsetprec(FP_PE)
373 #define RETURNI(x) do { \
374 __retval = (x); \
375 if (__oprec != FP_PE) \
376 fpsetprec(__oprec); \
377 RETURNF(__retval); \
378 } while (0)
379 #define ENTERV() \
380 fp_prec_t __oprec; \
381 \
382 if ((__oprec = fpgetprec()) != FP_PE) \
383 fpsetprec(FP_PE)
384 #define RETURNV() do { \
385 if (__oprec != FP_PE) \
386 fpsetprec(__oprec); \
387 return; \
388 } while (0)
389 #else
390 #define ENTERI()
391 #define ENTERIT(x)
392 #define RETURNI(x) RETURNF(x)
393 #define ENTERV()
394 #define RETURNV() return
395 #endif
396
397 /* Default return statement if hack*_t() is not used. */
398 #define RETURNF(v) return (v)
399
400 /*
401 * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
402 * a == 0, but is slower.
403 */
404 #define _2sum(a, b) do { \
405 __typeof(a) __s, __w; \
406 \
407 __w = (a) + (b); \
408 __s = __w - (a); \
409 (b) = ((a) - (__w - __s)) + ((b) - __s); \
410 (a) = __w; \
411 } while (0)
412
413 /*
414 * 2sumF algorithm.
415 *
416 * "Normalize" the terms in the infinite-precision expression a + b for
417 * the sum of 2 floating point values so that b is as small as possible
418 * relative to 'a'. (The resulting 'a' is the value of the expression in
419 * the same precision as 'a' and the resulting b is the rounding error.)
420 * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
421 * exponent overflow or underflow must not occur. This uses a Theorem of
422 * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
423 * is apparently due to Skewchuk (1997).
424 *
425 * For this to always work, assignment of a + b to 'a' must not retain any
426 * extra precision in a + b. This is required by C standards but broken
427 * in many compilers. The brokenness cannot be worked around using
428 * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
429 * algorithm would be destroyed by non-null strict assignments. (The
430 * compilers are correct to be broken -- the efficiency of all floating
431 * point code calculations would be destroyed similarly if they forced the
432 * conversions.)
433 *
434 * Fortunately, a case that works well can usually be arranged by building
435 * any extra precision into the type of 'a' -- 'a' should have type float_t,
436 * double_t or long double. b's type should be no larger than 'a's type.
437 * Callers should use these types with scopes as large as possible, to
438 * reduce their own extra-precision and efficiciency problems. In
439 * particular, they shouldn't convert back and forth just to call here.
440 */
441 #ifdef DEBUG
442 #define _2sumF(a, b) do { \
443 __typeof(a) __w; \
444 volatile __typeof(a) __ia, __ib, __r, __vw; \
445 \
446 __ia = (a); \
447 __ib = (b); \
448 assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
449 \
450 __w = (a) + (b); \
451 (b) = ((a) - __w) + (b); \
452 (a) = __w; \
453 \
454 /* The next 2 assertions are weak if (a) is already long double. */ \
455 assert((long double)__ia + __ib == (long double)(a) + (b)); \
456 __vw = __ia + __ib; \
457 __r = __ia - __vw; \
458 __r += __ib; \
459 assert(__vw == (a) && __r == (b)); \
460 } while (0)
461 #else /* !DEBUG */
462 #define _2sumF(a, b) do { \
463 __typeof(a) __w; \
464 \
465 __w = (a) + (b); \
466 (b) = ((a) - __w) + (b); \
467 (a) = __w; \
468 } while (0)
469 #endif /* DEBUG */
470
471 /*
472 * Set x += c, where x is represented in extra precision as a + b.
473 * x must be sufficiently normalized and sufficiently larger than c,
474 * and the result is then sufficiently normalized.
475 *
476 * The details of ordering are that |a| must be >= |c| (so that (a, c)
477 * can be normalized without extra work to swap 'a' with c). The details of
478 * the normalization are that b must be small relative to the normalized 'a'.
479 * Normalization of (a, c) makes the normalized c tiny relative to the
480 * normalized a, so b remains small relative to 'a' in the result. However,
481 * b need not ever be tiny relative to 'a'. For example, b might be about
482 * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
483 * That is usually enough, and adding c (which by normalization is about
484 * 2**53 times smaller than a) cannot change b significantly. However,
485 * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
486 * significantly relative to b. The caller must ensure that significant
487 * cancellation doesn't occur, either by having c of the same sign as 'a',
488 * or by having |c| a few percent smaller than |a|. Pre-normalization of
489 * (a, b) may help.
490 *
491 * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
492 * exercise 19). We gain considerable efficiency by requiring the terms to
493 * be sufficiently normalized and sufficiently increasing.
494 */
495 #define _3sumF(a, b, c) do { \
496 __typeof(a) __tmp; \
497 \
498 __tmp = (c); \
499 _2sumF(__tmp, (a)); \
500 (b) += (a); \
501 (a) = __tmp; \
502 } while (0)
503
504 /*
505 * Common routine to process the arguments to nan(), nanf(), and nanl().
506 */
507 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
508
509 /*
510 * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns
511 * signaling NaNs into quiet NaNs by setting a quiet bit. We do this
512 * because we want to never return a signaling NaN, and also because we
513 * don't want the quiet bit to affect the result. Then mix the converted
514 * args using the specified operation.
515 *
516 * When one arg is NaN, the result is typically that arg quieted. When both
517 * args are NaNs, the result is typically the quietening of the arg whose
518 * mantissa is largest after quietening. When neither arg is NaN, the
519 * result may be NaN because it is indeterminate, or finite for subsequent
520 * construction of a NaN as the indeterminate 0.0L/0.0L.
521 *
522 * Technical complications: the result in bits after rounding to the final
523 * precision might depend on the runtime precision and/or on compiler
524 * optimizations, especially when different register sets are used for
525 * different precisions. Try to make the result not depend on at least the
526 * runtime precision by always doing the main mixing step in long double
527 * precision. Try to reduce dependencies on optimizations by adding the
528 * the 0's in different precisions (unless everything is in long double
529 * precision).
530 */
531 #define nan_mix(x, y) (nan_mix_op((x), (y), +))
532 #define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0))
533
534 #ifdef _COMPLEX_H
535
536 /*
537 * C99 specifies that complex numbers have the same representation as
538 * an array of two elements, where the first element is the real part
539 * and the second element is the imaginary part.
540 */
541 typedef union {
542 float complex f;
543 float a[2];
544 } float_complex;
545 typedef union {
546 double complex f;
547 double a[2];
548 } double_complex;
549 typedef union {
550 long double complex f;
551 long double a[2];
552 } long_double_complex;
553 #define REALPART(z) ((z).a[0])
554 #define IMAGPART(z) ((z).a[1])
555
556 /*
557 * Inline functions that can be used to construct complex values.
558 *
559 * The C99 standard intends x+I*y to be used for this, but x+I*y is
560 * currently unusable in general since gcc introduces many overflow,
561 * underflow, sign and efficiency bugs by rewriting I*y as
562 * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
563 * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
564 * to -0.0+I*0.0.
565 *
566 * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
567 * to construct complex values. Compilers that conform to the C99
568 * standard require the following functions to avoid the above issues.
569 */
570
571 #ifndef CMPLXF
572 static __inline float complex
CMPLXF(float x,float y)573 CMPLXF(float x, float y)
574 {
575 float_complex z;
576
577 REALPART(z) = x;
578 IMAGPART(z) = y;
579 return (z.f);
580 }
581 #endif
582
583 #ifndef CMPLX
584 static __inline double complex
CMPLX(double x,double y)585 CMPLX(double x, double y)
586 {
587 double_complex z;
588
589 REALPART(z) = x;
590 IMAGPART(z) = y;
591 return (z.f);
592 }
593 #endif
594
595 #ifndef CMPLXL
596 static __inline long double complex
CMPLXL(long double x,long double y)597 CMPLXL(long double x, long double y)
598 {
599 long_double_complex z;
600
601 REALPART(z) = x;
602 IMAGPART(z) = y;
603 return (z.f);
604 }
605 #endif
606
607 #endif /* _COMPLEX_H */
608
609 /*
610 * The rnint() family rounds to the nearest integer for a restricted range
611 * range of args (up to about 2**MANT_DIG). We assume that the current
612 * rounding mode is FE_TONEAREST so that this can be done efficiently.
613 * Extra precision causes more problems in practice, and we only centralize
614 * this here to reduce those problems, and have not solved the efficiency
615 * problems. The exp2() family uses a more delicate version of this that
616 * requires extracting bits from the intermediate value, so it is not
617 * centralized here and should copy any solution of the efficiency problems.
618 */
619
620 static inline double
rnint(__double_t x)621 rnint(__double_t x)
622 {
623 /*
624 * This casts to double to kill any extra precision. This depends
625 * on the cast being applied to a double_t to avoid compiler bugs
626 * (this is a cleaner version of STRICT_ASSIGN()). This is
627 * inefficient if there actually is extra precision, but is hard
628 * to improve on. We use double_t in the API to minimise conversions
629 * for just calling here. Note that we cannot easily change the
630 * magic number to the one that works directly with double_t, since
631 * the rounding precision is variable at runtime on x86 so the
632 * magic number would need to be variable. Assuming that the
633 * rounding precision is always the default is too fragile. This
634 * and many other complications will move when the default is
635 * changed to FP_PE.
636 */
637 return ((double)(x + 0x1.8p52) - 0x1.8p52);
638 }
639
640 static inline float
rnintf(__float_t x)641 rnintf(__float_t x)
642 {
643 /*
644 * As for rnint(), except we could just call that to handle the
645 * extra precision case, usually without losing efficiency.
646 */
647 return ((float)(x + 0x1.8p23F) - 0x1.8p23F);
648 }
649
650 #ifdef LDBL_MANT_DIG
651 /*
652 * The complications for extra precision are smaller for rnintl() since it
653 * can safely assume that the rounding precision has been increased from
654 * its default to FP_PE on x86. We don't exploit that here to get small
655 * optimizations from limiting the rangle to double. We just need it for
656 * the magic number to work with long doubles. ld128 callers should use
657 * rnint() instead of this if possible. ld80 callers should prefer
658 * rnintl() since for amd64 this avoids swapping the register set, while
659 * for i386 it makes no difference (assuming FP_PE), and for other arches
660 * it makes little difference.
661 */
662 static inline long double
rnintl(long double x)663 rnintl(long double x)
664 {
665 return (x + __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2 -
666 __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2);
667 }
668 #endif /* LDBL_MANT_DIG */
669
670 /*
671 * irint() and i64rint() give the same result as casting to their integer
672 * return type provided their arg is a floating point integer. They can
673 * sometimes be more efficient because no rounding is required.
674 */
675 #if (defined(amd64) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
676 #define irint(x) \
677 (sizeof(x) == sizeof(float) && \
678 sizeof(__float_t) == sizeof(long double) ? irintf(x) : \
679 sizeof(x) == sizeof(double) && \
680 sizeof(__double_t) == sizeof(long double) ? irintd(x) : \
681 sizeof(x) == sizeof(long double) ? irintl(x) : (int)(x))
682 #else
683 #define irint(x) ((int)(x))
684 #endif
685
686 #define i64rint(x) ((int64_t)(x)) /* only needed for ld128 so not opt. */
687
688 #if defined(__i386__) && defined(__GNUCLIKE_ASM)
689 static __inline int
irintf(float x)690 irintf(float x)
691 {
692 int n;
693
694 __asm("fistl %0" : "=m" (n) : "t" (x));
695 return (n);
696 }
697
698 static __inline int
irintd(double x)699 irintd(double x)
700 {
701 int n;
702
703 __asm("fistl %0" : "=m" (n) : "t" (x));
704 return (n);
705 }
706 #endif
707
708 #if (defined(__amd64__) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
709 static __inline int
irintl(long double x)710 irintl(long double x)
711 {
712 int n;
713
714 __asm("fistl %0" : "=m" (n) : "t" (x));
715 return (n);
716 }
717 #endif
718
719 #ifdef DEBUG
720 #if defined(__amd64__) || defined(__i386__)
721 #define breakpoint() asm("int $3")
722 #else
723 #include <signal.h>
724
725 #define breakpoint() raise(SIGTRAP)
726 #endif
727 #endif
728
729 /* Write a pari script to test things externally. */
730 #ifdef DOPRINT
731 #include <stdio.h>
732
733 #ifndef DOPRINT_SWIZZLE
734 #define DOPRINT_SWIZZLE 0
735 #endif
736
737 #ifdef DOPRINT_LD80
738
739 #define DOPRINT_START(xp) do { \
740 uint64_t __lx; \
741 uint16_t __hx; \
742 \
743 /* Hack to give more-problematic args. */ \
744 EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
745 __lx ^= DOPRINT_SWIZZLE; \
746 INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
747 printf("x = %.21Lg; ", (long double)*xp); \
748 } while (0)
749 #define DOPRINT_END1(v) \
750 printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
751 #define DOPRINT_END2(hi, lo) \
752 printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
753 (long double)(hi), (long double)(lo))
754
755 #elif defined(DOPRINT_D64)
756
757 #define DOPRINT_START(xp) do { \
758 uint32_t __hx, __lx; \
759 \
760 EXTRACT_WORDS(__hx, __lx, *xp); \
761 __lx ^= DOPRINT_SWIZZLE; \
762 INSERT_WORDS(*xp, __hx, __lx); \
763 printf("x = %.21Lg; ", (long double)*xp); \
764 } while (0)
765 #define DOPRINT_END1(v) \
766 printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
767 #define DOPRINT_END2(hi, lo) \
768 printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
769 (long double)(hi), (long double)(lo))
770
771 #elif defined(DOPRINT_F32)
772
773 #define DOPRINT_START(xp) do { \
774 uint32_t __hx; \
775 \
776 GET_FLOAT_WORD(__hx, *xp); \
777 __hx ^= DOPRINT_SWIZZLE; \
778 SET_FLOAT_WORD(*xp, __hx); \
779 printf("x = %.21Lg; ", (long double)*xp); \
780 } while (0)
781 #define DOPRINT_END1(v) \
782 printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
783 #define DOPRINT_END2(hi, lo) \
784 printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
785 (long double)(hi), (long double)(lo))
786
787 #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
788
789 #ifndef DOPRINT_SWIZZLE_HIGH
790 #define DOPRINT_SWIZZLE_HIGH 0
791 #endif
792
793 #define DOPRINT_START(xp) do { \
794 uint64_t __lx, __llx; \
795 uint16_t __hx; \
796 \
797 EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
798 __llx ^= DOPRINT_SWIZZLE; \
799 __lx ^= DOPRINT_SWIZZLE_HIGH; \
800 INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
801 printf("x = %.36Lg; ", (long double)*xp); \
802 } while (0)
803 #define DOPRINT_END1(v) \
804 printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
805 #define DOPRINT_END2(hi, lo) \
806 printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
807 (long double)(hi), (long double)(lo))
808
809 #endif /* DOPRINT_LD80 */
810
811 #else /* !DOPRINT */
812 #define DOPRINT_START(xp)
813 #define DOPRINT_END1(v)
814 #define DOPRINT_END2(hi, lo)
815 #endif /* DOPRINT */
816
817 #define RETURNP(x) do { \
818 DOPRINT_END1(x); \
819 RETURNF(x); \
820 } while (0)
821 #define RETURNPI(x) do { \
822 DOPRINT_END1(x); \
823 RETURNI(x); \
824 } while (0)
825 #define RETURN2P(x, y) do { \
826 DOPRINT_END2((x), (y)); \
827 RETURNF((x) + (y)); \
828 } while (0)
829 #define RETURN2PI(x, y) do { \
830 DOPRINT_END2((x), (y)); \
831 RETURNI((x) + (y)); \
832 } while (0)
833 #ifdef STRUCT_RETURN
834 #define RETURNSP(rp) do { \
835 if (!(rp)->lo_set) \
836 RETURNP((rp)->hi); \
837 RETURN2P((rp)->hi, (rp)->lo); \
838 } while (0)
839 #define RETURNSPI(rp) do { \
840 if (!(rp)->lo_set) \
841 RETURNPI((rp)->hi); \
842 RETURN2PI((rp)->hi, (rp)->lo); \
843 } while (0)
844 #endif
845 #define SUM2P(x, y) ({ \
846 const __typeof (x) __x = (x); \
847 const __typeof (y) __y = (y); \
848 \
849 DOPRINT_END2(__x, __y); \
850 __x + __y; \
851 })
852
853 /*
854 * ieee style elementary functions
855 *
856 * We rename functions here to improve other sources' diffability
857 * against fdlibm.
858 */
859 #define __ieee754_sqrt sqrt
860 #define __ieee754_acos acos
861 #define __ieee754_acosh acosh
862 #define __ieee754_log log
863 #define __ieee754_log2 log2
864 #define __ieee754_atanh atanh
865 #define __ieee754_asin asin
866 #define __ieee754_atan2 atan2
867 #define __ieee754_exp exp
868 #define __ieee754_cosh cosh
869 #define __ieee754_fmod fmod
870 #define __ieee754_pow pow
871 #define __ieee754_lgamma lgamma
872 #define __ieee754_gamma gamma
873 #define __ieee754_lgamma_r lgamma_r
874 #define __ieee754_gamma_r gamma_r
875 #define __ieee754_log10 log10
876 #define __ieee754_sinh sinh
877 #define __ieee754_hypot hypot
878 #define __ieee754_j0 j0
879 #define __ieee754_j1 j1
880 #define __ieee754_y0 y0
881 #define __ieee754_y1 y1
882 #define __ieee754_jn jn
883 #define __ieee754_yn yn
884 #define __ieee754_remainder remainder
885 #define __ieee754_scalb scalb
886 #define __ieee754_sqrtf sqrtf
887 #define __ieee754_acosf acosf
888 #define __ieee754_acoshf acoshf
889 #define __ieee754_logf logf
890 #define __ieee754_atanhf atanhf
891 #define __ieee754_asinf asinf
892 #define __ieee754_atan2f atan2f
893 #define __ieee754_expf expf
894 #define __ieee754_coshf coshf
895 #define __ieee754_fmodf fmodf
896 #define __ieee754_powf powf
897 #define __ieee754_lgammaf lgammaf
898 #define __ieee754_gammaf gammaf
899 #define __ieee754_lgammaf_r lgammaf_r
900 #define __ieee754_gammaf_r gammaf_r
901 #define __ieee754_log10f log10f
902 #define __ieee754_log2f log2f
903 #define __ieee754_sinhf sinhf
904 #define __ieee754_hypotf hypotf
905 #define __ieee754_j0f j0f
906 #define __ieee754_j1f j1f
907 #define __ieee754_y0f y0f
908 #define __ieee754_y1f y1f
909 #define __ieee754_jnf jnf
910 #define __ieee754_ynf ynf
911 #define __ieee754_remainderf remainderf
912 #define __ieee754_scalbf scalbf
913
914 #define acos fdlibm::acos
915 #define asin fdlibm::asin
916 #define atan fdlibm::atan
917 #define atan2 fdlibm::atan2
918 #define cos fdlibm::cos
919 #define sin fdlibm::sin
920 #define tan fdlibm::tan
921 #define cosh fdlibm::cosh
922 #define sinh fdlibm::sinh
923 #define tanh fdlibm::tanh
924 #define exp fdlibm::exp
925 #define log fdlibm::log
926 #define log10 fdlibm::log10
927 #define pow fdlibm::pow
928 #define ceil fdlibm::ceil
929 #define ceilf fdlibm::ceilf
930 #define fabs fdlibm::fabs
931 #define floor fdlibm::floor
932 #define acosh fdlibm::acosh
933 #define asinh fdlibm::asinh
934 #define atanh fdlibm::atanh
935 #define cbrt fdlibm::cbrt
936 #define expm1 fdlibm::expm1
937 #define hypot fdlibm::hypot
938 #define log1p fdlibm::log1p
939 #define log2 fdlibm::log2
940 #define scalb fdlibm::scalb
941 #define copysign fdlibm::copysign
942 #define scalbn fdlibm::scalbn
943 #define trunc fdlibm::trunc
944 #define truncf fdlibm::truncf
945 #define floorf fdlibm::floorf
946 #define nearbyint fdlibm::nearbyint
947 #define nearbyintf fdlibm::nearbyintf
948 #define rint fdlibm::rint
949 #define rintf fdlibm::rintf
950
951 /* fdlibm kernel function */
952 int __kernel_rem_pio2(double*,double*,int,int,int);
953
954 /* double precision kernel functions */
955 #ifndef INLINE_REM_PIO2
956 int __ieee754_rem_pio2(double,double*);
957 #endif
958 double __kernel_sin(double,double,int);
959 double __kernel_cos(double,double);
960 double __kernel_tan(double,double,int);
961 double __ldexp_exp(double,int);
962 #ifdef _COMPLEX_H
963 double complex __ldexp_cexp(double complex,int);
964 #endif
965
966 /* float precision kernel functions */
967 #ifndef INLINE_REM_PIO2F
968 int __ieee754_rem_pio2f(float,double*);
969 #endif
970 #ifndef INLINE_KERNEL_SINDF
971 float __kernel_sindf(double);
972 #endif
973 #ifndef INLINE_KERNEL_COSDF
974 float __kernel_cosdf(double);
975 #endif
976 #ifndef INLINE_KERNEL_TANDF
977 float __kernel_tandf(double,int);
978 #endif
979 float __ldexp_expf(float,int);
980 #ifdef _COMPLEX_H
981 float complex __ldexp_cexpf(float complex,int);
982 #endif
983
984 /* long double precision kernel functions */
985 long double __kernel_sinl(long double, long double, int);
986 long double __kernel_cosl(long double, long double);
987 long double __kernel_tanl(long double, long double, int);
988
989 #endif /* !_MATH_PRIVATE_H_ */
990