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