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