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 * $NetBSD: math_private.h,v 1.34 2024/07/17 12:00:13 riastradh Exp $
15 */
16
17 #ifndef _MATH_PRIVATE_H_
18 #define _MATH_PRIVATE_H_
19
20 #include <assert.h>
21 #include <sys/types.h>
22
23 /*
24 * The original fdlibm code used statements like:
25 *
26 * n0 = ((*(int*)&one)>>29)^1; // index of high word
27 * ix0 = *(n0+(int*)&x); // high word of x
28 * ix1 = *((1-n0)+(int*)&x); // low word of x
29 *
30 * to dig two 32-bit words out of the-64 bit IEEE floating point value.
31 * That is non-ANSI, and, moreover, the gcc instruction scheduler gets
32 * it wrong. We instead use the following macros. Unlike the original
33 * code, we determine the endianness at compile time, not at run time;
34 * I don't see much benefit to selecting endianness at run time.
35 */
36
37 #ifdef __arm__
38 #if defined(__VFP_FP__) || defined(__ARM_EABI__)
39 #define IEEE_WORD_ORDER BYTE_ORDER
40 #else
41 #define IEEE_WORD_ORDER BIG_ENDIAN
42 #endif
43 #else /* __arm__ */
44 #define IEEE_WORD_ORDER BYTE_ORDER
45 #endif
46
47 /*
48 * A union which permits us to convert between a long double and
49 * four 32-bit integers.
50 */
51
52 #if IEEE_WORD_ORDER == BIG_ENDIAN
53
54 typedef union
55 {
56 long double value;
57 struct {
58 u_int32_t mswhi;
59 u_int32_t mswlo;
60 u_int32_t lswhi;
61 u_int32_t lswlo;
62 } parts32;
63 struct {
64 u_int64_t msw;
65 u_int64_t lsw;
66 } parts64;
67 } ieee_quad_shape_type;
68
69 #endif
70
71 #if IEEE_WORD_ORDER == LITTLE_ENDIAN
72
73 typedef union
74 {
75 long double value;
76 struct {
77 u_int32_t lswlo;
78 u_int32_t lswhi;
79 u_int32_t mswlo;
80 u_int32_t mswhi;
81 } parts32;
82 struct {
83 u_int64_t lsw;
84 u_int64_t msw;
85 } parts64;
86 } ieee_quad_shape_type;
87
88 #endif
89
90 /*
91 * A union which permits us to convert between a double and two 32-bit
92 * integers.
93 */
94
95 #if IEEE_WORD_ORDER == BIG_ENDIAN
96
97 typedef union
98 {
99 double value;
100 struct
101 {
102 u_int32_t msw;
103 u_int32_t lsw;
104 } parts;
105 struct
106 {
107 u_int64_t w;
108 } xparts;
109 } ieee_double_shape_type;
110
111 #endif
112
113 #if IEEE_WORD_ORDER == LITTLE_ENDIAN
114
115 typedef union
116 {
117 double value;
118 struct
119 {
120 u_int32_t lsw;
121 u_int32_t msw;
122 } parts;
123 struct
124 {
125 u_int64_t w;
126 } xparts;
127 } ieee_double_shape_type;
128
129 #endif
130
131 /* Get two 32-bit integers from a double. */
132
133 #define EXTRACT_WORDS(ix0,ix1,d) \
134 do { \
135 ieee_double_shape_type ew_u; \
136 ew_u.value = (d); \
137 (ix0) = ew_u.parts.msw; \
138 (ix1) = ew_u.parts.lsw; \
139 } while (0)
140
141 /* Get a 64-bit integer from a double. */
142 #define EXTRACT_WORD64(ix,d) \
143 do { \
144 ieee_double_shape_type ew_u; \
145 ew_u.value = (d); \
146 (ix) = ew_u.xparts.w; \
147 } while (0)
148
149
150 /* Get the more significant 32-bit integer from a double. */
151
152 #define GET_HIGH_WORD(i,d) \
153 do { \
154 ieee_double_shape_type gh_u; \
155 gh_u.value = (d); \
156 (i) = gh_u.parts.msw; \
157 } while (0)
158
159 /* Get the less significant 32-bit integer from a double. */
160
161 #define GET_LOW_WORD(i,d) \
162 do { \
163 ieee_double_shape_type gl_u; \
164 gl_u.value = (d); \
165 (i) = gl_u.parts.lsw; \
166 } while (0)
167
168 /* Set a double from two 32-bit integers. */
169
170 #define INSERT_WORDS(d,ix0,ix1) \
171 do { \
172 ieee_double_shape_type iw_u; \
173 iw_u.parts.msw = (ix0); \
174 iw_u.parts.lsw = (ix1); \
175 (d) = iw_u.value; \
176 } while (0)
177
178 /* Set a double from a 64-bit integer. */
179
180 #define INSERT_WORD64(d,ix) \
181 do { \
182 ieee_double_shape_type iw_u; \
183 iw_u.xparts.w = (ix); \
184 (d) = iw_u.value; \
185 } while (0)
186
187 /* Set the more significant 32 bits of a double from an integer. */
188
189 #define SET_HIGH_WORD(d,v) \
190 do { \
191 ieee_double_shape_type sh_u; \
192 sh_u.value = (d); \
193 sh_u.parts.msw = (v); \
194 (d) = sh_u.value; \
195 } while (0)
196
197 /* Set the less significant 32 bits of a double from an integer. */
198
199 #define SET_LOW_WORD(d,v) \
200 do { \
201 ieee_double_shape_type sl_u; \
202 sl_u.value = (d); \
203 sl_u.parts.lsw = (v); \
204 (d) = sl_u.value; \
205 } while (0)
206
207 /*
208 * A union which permits us to convert between a float and a 32-bit
209 * integer.
210 */
211
212 typedef union
213 {
214 float value;
215 u_int32_t word;
216 } ieee_float_shape_type;
217
218 /* Get a 32-bit integer from a float. */
219
220 #define GET_FLOAT_WORD(i,d) \
221 do { \
222 ieee_float_shape_type gf_u; \
223 gf_u.value = (d); \
224 (i) = gf_u.word; \
225 } while (0)
226
227 /* Set a float from a 32-bit integer. */
228
229 #define SET_FLOAT_WORD(d,i) \
230 do { \
231 ieee_float_shape_type sf_u; \
232 sf_u.word = (i); \
233 (d) = sf_u.value; \
234 } while (0)
235
236 #define GET_EXPSIGN(u) \
237 (((u)->extu_sign << EXT_EXPBITS) | (u)->extu_exp)
238 #define SET_EXPSIGN(u, x) \
239 (u)->extu_exp = (x), \
240 (u)->extu_sign = ((x) >> EXT_EXPBITS)
241 #define GET_LDBL80_MAN(u) \
242 (((uint64_t)(u)->extu_frach << EXT_FRACLBITS) | (u)->extu_fracl)
243 #define SET_LDBL80_MAN(u, v) \
244 ((u)->extu_fracl = (v) & ((1ULL << EXT_FRACLBITS) - 1), \
245 (u)->extu_frach = (v) >> EXT_FRACLBITS)
246
247 /*
248 * Get expsign as 16-bit integer ix0 and significand as 64-bit integer
249 * ix1 from an 80-bit long double d.
250 */
251
252 #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \
253 do { \
254 union ieee_ext_u ew_u; \
255 ew_u.extu_ld = (d); \
256 (ix0) = GET_EXPSIGN(&ew_u); \
257 (ix1) = GET_LDBL80_MAN(&ew_u); \
258 } while (0)
259
260 /*
261 * Get expsign as 16-bit integer ix0 and significand as two 64-bit
262 * integers, ix1 high-order and ix2 low-order, from a 128-bit long
263 * double d.
264 */
265
266 #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \
267 do { \
268 union ieee_ext_u ew_u; \
269 ew_u.extu_ld = (d); \
270 (ix0) = GET_EXPSIGN(&ew_u); \
271 (ix1) = ew_u.extu_frach; \
272 (ix2) = ew_u.extu_fracl; \
273 } while (0)
274
275 /* Get expsign as a 16-bit integer i from a long double d. */
276
277 #define GET_LDBL_EXPSIGN(i,d) \
278 do { \
279 union ieee_ext_u ge_u; \
280 ge_u.extu_ld = (d); \
281 (i) = GET_EXPSIGN(&ge_u); \
282 } while (0)
283
284 /*
285 * Set an 80-bit long double d from a 16-bit integer expsign ix0 and a
286 * 64-bit integer significand ix1.
287 */
288
289 #define INSERT_LDBL80_WORDS(d,ix0,ix1) \
290 do { \
291 union ieee_ext_u iw_u; \
292 SET_EXPSIGN(&iw_u, ix0); \
293 SET_LDBL80_MAN(&iw_u, ix1); \
294 (d) = iw_u.extu_ld; \
295 } while (0)
296
297 /*
298 * Set a 128-bit long double d from a 16-bit integer expsign ix0 and
299 * two 64-bit integers composing the significand, ix1 high-order and
300 * ix2 low-order.
301 */
302
303 #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \
304 do { \
305 union ieee_ext_u iw_u; \
306 SET_EXPSIGN(&iw_u, ix0); \
307 iw_u.extu_frach = (ix1); \
308 iw_u.extu_fracl = (ix2); \
309 (d) = iw_u.extu_ld; \
310 } while (0)
311
312 /* Set expsign of a long double from a 16-bit integer. */
313
314 #define SET_LDBL_EXPSIGN(d,v) \
315 do { \
316 union ieee_ext_u se_u; \
317 se_u.extu_ld = (d); \
318 SET_EXPSIGN(&se_u, v); \
319 (d) = se_u.extu_ld; \
320 } while (0)
321
322 #ifdef __i386__
323 /* Long double constants are broken on i386. */
324 #define LD80C(m, ex, v) { \
325 .extu_fracl = (uint32_t)(__CONCAT(m, ULL)), \
326 .extu_frach = __CONCAT(m, ULL) >> EXT_FRACLBITS, \
327 .extu_exp = (0x3fff + (ex)), \
328 .extu_sign = ((v) < 0), \
329 }
330 #else
331 /**XXX: the following comment may no longer be true: kre 20240122 **/
332 /* The above works on non-i386 too, but we use this to check v. */
333 #define LD80C(m, ex, v) { .extu_ld = (v), }
334 #endif
335
336 /*
337 * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
338 */
339 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
340 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
341 #else
342 #define STRICT_ASSIGN(type, lval, rval) do { \
343 volatile type __lval; \
344 \
345 if (sizeof(type) >= sizeof(long double)) \
346 (lval) = (rval); \
347 else { \
348 __lval = (rval); \
349 (lval) = __lval; \
350 } \
351 } while (0)
352 #endif
353
354 /* Support switching the mode to FP_PE if necessary. */
355 #if defined(__i386__) && !defined(NO_FPSETPREC)
356
357 #include <ieeefp.h>
358
359 #define ENTERI() ENTERIT(long double)
360 #define ENTERIT(returntype) \
361 returntype __retval; \
362 fp_prec_t __oprec; \
363 \
364 if ((__oprec = fpgetprec()) != FP_PE) \
365 fpsetprec(FP_PE)
366 #define RETURNI(x) do { \
367 __retval = (x); \
368 if (__oprec != FP_PE) \
369 fpsetprec(__oprec); \
370 RETURNF(__retval); \
371 } while (0)
372 #define ENTERV() \
373 fp_prec_t __oprec; \
374 \
375 if ((__oprec = fpgetprec()) != FP_PE) \
376 fpsetprec(FP_PE)
377 #define RETURNV() do { \
378 if (__oprec != FP_PE) \
379 fpsetprec(__oprec); \
380 return; \
381 } while (0)
382 #else
383 #define ENTERI()
384 #define ENTERIT(x)
385 #define RETURNI(x) RETURNF(x)
386 #define ENTERV()
387 #define RETURNV() return
388 #endif
389
390 /* Default return statement if hack*_t() is not used. */
391 #define RETURNF(v) return (v)
392
393 /*
394 * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
395 * a == 0, but is slower.
396 */
397 #define _2sum(a, b) do { \
398 __typeof(a) __s, __w; \
399 \
400 __w = (a) + (b); \
401 __s = __w - (a); \
402 (b) = ((a) - (__w - __s)) + ((b) - __s); \
403 (a) = __w; \
404 } while (0)
405
406 /*
407 * 2sumF algorithm.
408 *
409 * "Normalize" the terms in the infinite-precision expression a + b for
410 * the sum of 2 floating point values so that b is as small as possible
411 * relative to 'a'. (The resulting 'a' is the value of the expression in
412 * the same precision as 'a' and the resulting b is the rounding error.)
413 * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
414 * exponent overflow or underflow must not occur. This uses a Theorem of
415 * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum"
416 * is apparently due to Skewchuk (1997).
417 *
418 * For this to always work, assignment of a + b to 'a' must not retain any
419 * extra precision in a + b. This is required by C standards but broken
420 * in many compilers. The brokenness cannot be worked around using
421 * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
422 * algorithm would be destroyed by non-null strict assignments. (The
423 * compilers are correct to be broken -- the efficiency of all floating
424 * point code calculations would be destroyed similarly if they forced the
425 * conversions.)
426 *
427 * Fortunately, a case that works well can usually be arranged by building
428 * any extra precision into the type of 'a' -- 'a' should have type float_t,
429 * double_t or long double. b's type should be no larger than 'a's type.
430 * Callers should use these types with scopes as large as possible, to
431 * reduce their own extra-precision and efficiciency problems. In
432 * particular, they shouldn't convert back and forth just to call here.
433 */
434 #ifdef DEBUG
435 #define _2sumF(a, b) do { \
436 __typeof(a) __w; \
437 volatile __typeof(a) __ia, __ib, __r, __vw; \
438 \
439 __ia = (a); \
440 __ib = (b); \
441 assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \
442 \
443 __w = (a) + (b); \
444 (b) = ((a) - __w) + (b); \
445 (a) = __w; \
446 \
447 /* The next 2 assertions are weak if (a) is already long double. */ \
448 assert((long double)__ia + __ib == (long double)(a) + (b)); \
449 __vw = __ia + __ib; \
450 __r = __ia - __vw; \
451 __r += __ib; \
452 assert(__vw == (a) && __r == (b)); \
453 } while (0)
454 #else /* !DEBUG */
455 #define _2sumF(a, b) do { \
456 __typeof(a) __w; \
457 \
458 __w = (a) + (b); \
459 (b) = ((a) - __w) + (b); \
460 (a) = __w; \
461 } while (0)
462 #endif /* DEBUG */
463
464 /*
465 * Set x += c, where x is represented in extra precision as a + b.
466 * x must be sufficiently normalized and sufficiently larger than c,
467 * and the result is then sufficiently normalized.
468 *
469 * The details of ordering are that |a| must be >= |c| (so that (a, c)
470 * can be normalized without extra work to swap 'a' with c). The details of
471 * the normalization are that b must be small relative to the normalized 'a'.
472 * Normalization of (a, c) makes the normalized c tiny relative to the
473 * normalized a, so b remains small relative to 'a' in the result. However,
474 * b need not ever be tiny relative to 'a'. For example, b might be about
475 * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
476 * That is usually enough, and adding c (which by normalization is about
477 * 2**53 times smaller than a) cannot change b significantly. However,
478 * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
479 * significantly relative to b. The caller must ensure that significant
480 * cancellation doesn't occur, either by having c of the same sign as 'a',
481 * or by having |c| a few percent smaller than |a|. Pre-normalization of
482 * (a, b) may help.
483 *
484 * This is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
485 * exercise 19). We gain considerable efficiency by requiring the terms to
486 * be sufficiently normalized and sufficiently increasing.
487 */
488 #define _3sumF(a, b, c) do { \
489 __typeof(a) __tmp; \
490 \
491 __tmp = (c); \
492 _2sumF(__tmp, (a)); \
493 (b) += (a); \
494 (a) = __tmp; \
495 } while (0)
496
497 /*
498 * Common routine to process the arguments to nan(), nanf(), and nanl().
499 */
500 void _scan_nan(uint32_t *__words, int __num_words, const char *__s);
501
502 /*
503 * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns
504 * signaling NaNs into quiet NaNs by setting a quiet bit. We do this
505 * because we want to never return a signaling NaN, and also because we
506 * don't want the quiet bit to affect the result. Then mix the converted
507 * args using the specified operation.
508 *
509 * When one arg is NaN, the result is typically that arg quieted. When both
510 * args are NaNs, the result is typically the quietening of the arg whose
511 * significand is largest after quietening. When neither arg is NaN, the
512 * result may be NaN because it is indeterminate, or finite for subsequent
513 * construction of a NaN as the indeterminate 0.0L/0.0L.
514 *
515 * Technical complications: the result in bits after rounding to the final
516 * precision might depend on the runtime precision and/or on compiler
517 * optimizations, especially when different register sets are used for
518 * different precisions. Try to make the result not depend on at least the
519 * runtime precision by always doing the main mixing step in long double
520 * precision. Try to reduce dependencies on optimizations by adding the
521 * the 0's in different precisions (unless everything is in long double
522 * precision).
523 */
524 #define nan_mix(x, y) (nan_mix_op((x), (y), +))
525 #define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0))
526
527 #ifdef _COMPLEX_H
528
529 /*
530 * Quoting from ISO/IEC 9899:TC2:
531 *
532 * 6.2.5.13 Types
533 * Each complex type has the same representation and alignment requirements as
534 * an array type containing exactly two elements of the corresponding real type;
535 * the first element is equal to the real part, and the second element to the
536 * imaginary part, of the complex number.
537 */
538 typedef union {
539 float complex z;
540 float parts[2];
541 } float_complex;
542
543 typedef union {
544 double complex z;
545 double parts[2];
546 } double_complex;
547
548 typedef union {
549 long double complex z;
550 long double parts[2];
551 } long_double_complex;
552
553 #define REAL_PART(z) ((z).parts[0])
554 #define IMAG_PART(z) ((z).parts[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 REAL_PART(z) = x;
578 IMAG_PART(z) = y;
579 return (z.z);
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 REAL_PART(z) = x;
590 IMAG_PART(z) = y;
591 return (z.z);
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 REAL_PART(z) = x;
602 IMAG_PART(z) = y;
603 return (z.z);
604 }
605 #endif
606
607 #endif /* _COMPLEX_H */
608
609 /* ieee style elementary functions */
610 extern double __ieee754_sqrt __P((double));
611 extern double __ieee754_acos __P((double));
612 extern double __ieee754_acosh __P((double));
613 extern double __ieee754_log __P((double));
614 extern double __ieee754_atanh __P((double));
615 extern double __ieee754_asin __P((double));
616 extern double __ieee754_atan2 __P((double,double));
617 extern double __ieee754_exp __P((double));
618 extern double __ieee754_cosh __P((double));
619 extern double __ieee754_fmod __P((double,double));
620 extern double __ieee754_pow __P((double,double));
621 extern double __ieee754_lgamma_r __P((double,int *));
622 extern double __ieee754_gamma_r __P((double,int *));
623 extern double __ieee754_lgamma __P((double));
624 extern double __ieee754_gamma __P((double));
625 extern double __ieee754_log10 __P((double));
626 extern double __ieee754_log2 __P((double));
627 extern double __ieee754_sinh __P((double));
628 extern double __ieee754_hypot __P((double,double));
629 extern double __ieee754_j0 __P((double));
630 extern double __ieee754_j1 __P((double));
631 extern double __ieee754_y0 __P((double));
632 extern double __ieee754_y1 __P((double));
633 extern double __ieee754_jn __P((int,double));
634 extern double __ieee754_yn __P((int,double));
635 extern double __ieee754_remainder __P((double,double));
636 extern int32_t __ieee754_rem_pio2 __P((double,double*));
637 extern double __ieee754_scalb __P((double,double));
638
639 /* fdlibm kernel function */
640 extern double __kernel_standard __P((double,double,int));
641 extern double __kernel_sin __P((double,double,int));
642 extern double __kernel_cos __P((double,double));
643 extern double __kernel_tan __P((double,double,int));
644 extern int __kernel_rem_pio2 __P((double*,double*,int,int,int));
645
646
647 /* ieee style elementary float functions */
648 extern float __ieee754_sqrtf __P((float));
649 extern float __ieee754_acosf __P((float));
650 extern float __ieee754_acoshf __P((float));
651 extern float __ieee754_logf __P((float));
652 extern float __ieee754_atanhf __P((float));
653 extern float __ieee754_asinf __P((float));
654 extern float __ieee754_atan2f __P((float,float));
655 extern float __ieee754_expf __P((float));
656 extern float __ieee754_coshf __P((float));
657 extern float __ieee754_fmodf __P((float,float));
658 extern float __ieee754_powf __P((float,float));
659 extern float __ieee754_lgammaf_r __P((float,int *));
660 extern float __ieee754_gammaf_r __P((float,int *));
661 extern float __ieee754_lgammaf __P((float));
662 extern float __ieee754_gammaf __P((float));
663 extern float __ieee754_log10f __P((float));
664 extern float __ieee754_log2f __P((float));
665 extern float __ieee754_sinhf __P((float));
666 extern float __ieee754_hypotf __P((float,float));
667 extern float __ieee754_j0f __P((float));
668 extern float __ieee754_j1f __P((float));
669 extern float __ieee754_y0f __P((float));
670 extern float __ieee754_y1f __P((float));
671 extern float __ieee754_jnf __P((int,float));
672 extern float __ieee754_ynf __P((int,float));
673 extern float __ieee754_remainderf __P((float,float));
674 extern int32_t __ieee754_rem_pio2f __P((float,float*));
675 extern float __ieee754_scalbf __P((float,float));
676
677 /* float versions of fdlibm kernel functions */
678 extern float __kernel_sinf __P((float,float,int));
679 extern float __kernel_cosf __P((float,float));
680 extern float __kernel_tanf __P((float,float,int));
681 extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const int32_t*));
682
683 /* ieee style elementary long double functions */
684 extern long double __ieee754_fmodl(long double, long double);
685 extern long double __ieee754_sqrtl(long double);
686
687 /*
688 * TRUNC() is a macro that sets the trailing 27 bits in the significand of an
689 * IEEE double variable to zero. It must be expression-like for syntactic
690 * reasons, and we implement this expression using an inline function
691 * instead of a pure macro to avoid depending on the gcc feature of
692 * statement-expressions.
693 */
694 #define TRUNC(d) (_b_trunc(&(d)))
695
696 static __inline void
_b_trunc(volatile double * _dp)697 _b_trunc(volatile double *_dp)
698 {
699 uint32_t _lw;
700
701 GET_LOW_WORD(_lw, *_dp);
702 SET_LOW_WORD(*_dp, _lw & 0xf8000000);
703 }
704
705 struct Double {
706 double a;
707 double b;
708 };
709
710 /*
711 * Functions internal to the math package, yet not static.
712 */
713 double __exp__D(double, double);
714 struct Double __log__D(double);
715
716 /*
717 * The rnint() family rounds to the nearest integer for a restricted range
718 * range of args (up to about 2**MANT_DIG). We assume that the current
719 * rounding mode is FE_TONEAREST so that this can be done efficiently.
720 * Extra precision causes more problems in practice, and we only centralize
721 * this here to reduce those problems, and have not solved the efficiency
722 * problems. The exp2() family uses a more delicate version of this that
723 * requires extracting bits from the intermediate value, so it is not
724 * centralized here and should copy any solution of the efficiency problems.
725 */
726
727 static inline double
rnint(double x)728 rnint(double x)
729 {
730 /*
731 * This casts to double to kill any extra precision. This depends
732 * on the cast being applied to a double_t to avoid compiler bugs
733 * (this is a cleaner version of STRICT_ASSIGN()). This is
734 * inefficient if there actually is extra precision, but is hard
735 * to improve on. We use double_t in the API to minimise conversions
736 * for just calling here. Note that we cannot easily change the
737 * magic number to the one that works directly with double_t, since
738 * the rounding precision is variable at runtime on x86 so the
739 * magic number would need to be variable. Assuming that the
740 * rounding precision is always the default is too fragile. This
741 * and many other complications will move when the default is
742 * changed to FP_PE.
743 */
744 return ((double)(x + 0x1.8p52) - 0x1.8p52);
745 }
746
747 static inline float
rnintf(float x)748 rnintf(float x)
749 {
750 /*
751 * As for rnint(), except we could just call that to handle the
752 * extra precision case, usually without losing efficiency.
753 */
754 return ((float)(x + 0x1.8p23F) - 0x1.8p23F);
755 }
756
757 #ifdef LDBL_MANT_DIG
758 /*
759 * The complications for extra precision are smaller for rnintl() since it
760 * can safely assume that the rounding precision has been increased from
761 * its default to FP_PE on x86. We don't exploit that here to get small
762 * optimizations from limiting the range to double. We just need it for
763 * the magic number to work with long doubles. ld128 callers should use
764 * rnint() instead of this if possible. ld80 callers should prefer
765 * rnintl() since for amd64 this avoids swapping the register set, while
766 * for i386 it makes no difference (assuming FP_PE), and for other arches
767 * it makes little difference.
768 */
769 static inline long double
rnintl(long double x)770 rnintl(long double x)
771 {
772 return (x + ___CONCAT(0x1.8p,LDBL_MANT_DIG) / 2 -
773 ___CONCAT(0x1.8p,LDBL_MANT_DIG) / 2);
774 }
775 #endif /* LDBL_MANT_DIG */
776
777 /*
778 * irint() and i64rint() give the same result as casting to their integer
779 * return type provided their arg is a floating point integer. They can
780 * sometimes be more efficient because no rounding is required.
781 */
782 #if (defined(amd64) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
783 #define irint(x) \
784 (sizeof(x) == sizeof(float) && \
785 sizeof(__float_t) == sizeof(long double) ? irintf(x) : \
786 sizeof(x) == sizeof(double) && \
787 sizeof(__double_t) == sizeof(long double) ? irintd(x) : \
788 sizeof(x) == sizeof(long double) ? irintl(x) : (int)(x))
789 #else
790 #define irint(x) ((int)(x))
791 #endif
792
793 #define i64rint(x) ((int64_t)(x)) /* only needed for ld128 so not opt. */
794
795 #if defined(__i386__) && defined(__GNUCLIKE_ASM)
796 static __inline int
irintf(float x)797 irintf(float x)
798 {
799 int n;
800
801 __asm("fistl %0" : "=m" (n) : "t" (x));
802 return (n);
803 }
804
805 static __inline int
irintd(double x)806 irintd(double x)
807 {
808 int n;
809
810 __asm("fistl %0" : "=m" (n) : "t" (x));
811 return (n);
812 }
813 #endif
814
815 #if (defined(__amd64__) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
816 static __inline int
irintl(long double x)817 irintl(long double x)
818 {
819 int n;
820
821 __asm("fistl %0" : "=m" (n) : "t" (x));
822 return (n);
823 }
824 #endif
825
826 /*
827 * The following are fast floor macros for 0 <= |x| < 0x1p(N-1), where
828 * N is the precision of the type of x. These macros are used in the
829 * half-cycle trignometric functions (e.g., sinpi(x)).
830 */
831 #define FFLOORF(x, j0, ix) do { \
832 (j0) = (((ix) >> 23) & 0xff) - 0x7f; \
833 (ix) &= ~(0x007fffff >> (j0)); \
834 SET_FLOAT_WORD((x), (ix)); \
835 } while (0)
836
837 #define FFLOOR(x, j0, ix, lx) do { \
838 (j0) = (((ix) >> 20) & 0x7ff) - 0x3ff; \
839 if ((j0) < 20) { \
840 (ix) &= ~(0x000fffff >> (j0)); \
841 (lx) = 0; \
842 } else { \
843 (lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20)); \
844 } \
845 INSERT_WORDS((x), (ix), (lx)); \
846 } while (0)
847
848 #define FFLOORL80(x, j0, ix, lx) do { \
849 j0 = ix - 0x3fff + 1; \
850 if ((j0) < 32) { \
851 (lx) = ((lx) >> 32) << 32; \
852 (lx) &= ~((((lx) << 32)-1) >> (j0)); \
853 } else { \
854 uint64_t _m; \
855 _m = (uint64_t)-1 >> (j0); \
856 if ((lx) & _m) (lx) &= ~_m; \
857 } \
858 INSERT_LDBL80_WORDS((x), (ix), (lx)); \
859 } while (0)
860
861 #define FFLOORL128(x, ai, ar) do { \
862 union ieee_ext_u u; \
863 uint64_t m; \
864 int e; \
865 u.extu_ld = (x); \
866 e = u.extu_exp - 16383; \
867 if (e < 48) { \
868 m = ((1llu << 49) - 1) >> (e + 1); \
869 u.extu_frach &= ~m; \
870 u.extu_fracl = 0; \
871 } else { \
872 m = (uint64_t)-1 >> (e - 48); \
873 u.extu_fracl &= ~m; \
874 } \
875 (ai) = u.extu_ld; \
876 (ar) = (x) - (ai); \
877 } while (0)
878
879 #ifdef DEBUG
880 #if defined(__amd64__) || defined(__i386__)
881 #define breakpoint() asm("int $3")
882 #else
883 #include <signal.h>
884
885 #define breakpoint() raise(SIGTRAP)
886 #endif
887 #endif
888
889 #ifdef STRUCT_RETURN
890 #define RETURNSP(rp) do { \
891 if (!(rp)->lo_set) \
892 RETURNF((rp)->hi); \
893 RETURNF((rp)->hi + (rp)->lo); \
894 } while (0)
895 #define RETURNSPI(rp) do { \
896 if (!(rp)->lo_set) \
897 RETURNI((rp)->hi); \
898 RETURNI((rp)->hi + (rp)->lo); \
899 } while (0)
900 #endif
901
902 #define SUM2P(x, y) ({ \
903 const __typeof (x) __x = (x); \
904 const __typeof (y) __y = (y); \
905 __x + __y; \
906 })
907
908 #ifndef INLINE_KERNEL_SINDF
909 float __kernel_sindf(double);
910 #endif
911 #ifndef INLINE_KERNEL_COSDF
912 float __kernel_cosdf(double);
913 #endif
914 #ifndef INLINE_KERNEL_TANDF
915 float __kernel_tandf(double,int);
916 #endif
917
918 /* long double precision kernel functions */
919 long double __kernel_sinl(long double, long double, int);
920 long double __kernel_cosl(long double, long double);
921 long double __kernel_tanl(long double, long double, int);
922
923 #endif /* _MATH_PRIVATE_H_ */
924