rjw | 1f88458 | 2022-01-06 17:20:42 +0800 | [diff] [blame^] | 1 | /* |
| 2 | * ==================================================== |
| 3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 4 | * |
| 5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
| 6 | * Permission to use, copy, modify, and distribute this |
| 7 | * software is freely granted, provided that this notice |
| 8 | * is preserved. |
| 9 | * ==================================================== |
| 10 | */ |
| 11 | |
| 12 | /* |
| 13 | * from: @(#)fdlibm.h 5.1 93/09/24 |
| 14 | * $FreeBSD$ |
| 15 | */ |
| 16 | |
| 17 | #ifndef _MATH_PRIVATE_H_ |
| 18 | #define _MATH_PRIVATE_H_ |
| 19 | |
| 20 | #include <sys/types.h> |
| 21 | #include <endian.h> |
| 22 | |
| 23 | /* |
| 24 | * The original fdlibm code used statements like: |
| 25 | * n0 = ((*(int*)&one)>>29)^1; * index of high word * |
| 26 | * ix0 = *(n0+(int*)&x); * high word of x * |
| 27 | * ix1 = *((1-n0)+(int*)&x); * low word of x * |
| 28 | * to dig two 32 bit words out of the 64 bit IEEE floating point |
| 29 | * value. That is non-ANSI, and, moreover, the gcc instruction |
| 30 | * scheduler gets it wrong. We instead use the following macros. |
| 31 | * Unlike the original code, we determine the endianness at compile |
| 32 | * time, not at run time; I don't see much benefit to selecting |
| 33 | * endianness at run time. |
| 34 | */ |
| 35 | |
| 36 | /* |
| 37 | * A union which permits us to convert between a double and two 32 bit |
| 38 | * ints. |
| 39 | */ |
| 40 | |
| 41 | #ifdef __arm__ |
| 42 | #if defined(__VFP_FP__) || defined(__ARM_EABI__) |
| 43 | #define IEEE_WORD_ORDER BYTE_ORDER |
| 44 | #else |
| 45 | #define IEEE_WORD_ORDER BIG_ENDIAN |
| 46 | #endif |
| 47 | #else /* __arm__ */ |
| 48 | #define IEEE_WORD_ORDER BYTE_ORDER |
| 49 | #endif |
| 50 | |
| 51 | #if IEEE_WORD_ORDER == BIG_ENDIAN |
| 52 | |
| 53 | typedef union { |
| 54 | double value; |
| 55 | struct { |
| 56 | u_int32_t msw; |
| 57 | u_int32_t lsw; |
| 58 | } parts; |
| 59 | struct { |
| 60 | u_int64_t w; |
| 61 | } xparts; |
| 62 | } ieee_double_shape_type; |
| 63 | |
| 64 | #endif |
| 65 | |
| 66 | #if IEEE_WORD_ORDER == LITTLE_ENDIAN |
| 67 | |
| 68 | typedef union { |
| 69 | double value; |
| 70 | struct { |
| 71 | u_int32_t lsw; |
| 72 | u_int32_t msw; |
| 73 | } parts; |
| 74 | struct { |
| 75 | u_int64_t w; |
| 76 | } xparts; |
| 77 | } ieee_double_shape_type; |
| 78 | |
| 79 | #endif |
| 80 | |
| 81 | /* Get two 32 bit ints from a double. */ |
| 82 | |
| 83 | #define EXTRACT_WORDS(ix0,ix1,d) \ |
| 84 | do { \ |
| 85 | ieee_double_shape_type ew_u; \ |
| 86 | ew_u.value = (d); \ |
| 87 | (ix0) = ew_u.parts.msw; \ |
| 88 | (ix1) = ew_u.parts.lsw; \ |
| 89 | } while (0) |
| 90 | |
| 91 | /* Get a 64-bit int from a double. */ |
| 92 | #define EXTRACT_WORD64(ix,d) \ |
| 93 | do { \ |
| 94 | ieee_double_shape_type ew_u; \ |
| 95 | ew_u.value = (d); \ |
| 96 | (ix) = ew_u.xparts.w; \ |
| 97 | } while (0) |
| 98 | |
| 99 | /* Get the more significant 32 bit int from a double. */ |
| 100 | |
| 101 | #define GET_HIGH_WORD(i,d) \ |
| 102 | do { \ |
| 103 | ieee_double_shape_type gh_u; \ |
| 104 | gh_u.value = (d); \ |
| 105 | (i) = gh_u.parts.msw; \ |
| 106 | } while (0) |
| 107 | |
| 108 | /* Get the less significant 32 bit int from a double. */ |
| 109 | |
| 110 | #define GET_LOW_WORD(i,d) \ |
| 111 | do { \ |
| 112 | ieee_double_shape_type gl_u; \ |
| 113 | gl_u.value = (d); \ |
| 114 | (i) = gl_u.parts.lsw; \ |
| 115 | } while (0) |
| 116 | |
| 117 | /* Set a double from two 32 bit ints. */ |
| 118 | |
| 119 | #define INSERT_WORDS(d,ix0,ix1) \ |
| 120 | do { \ |
| 121 | ieee_double_shape_type iw_u; \ |
| 122 | iw_u.parts.msw = (ix0); \ |
| 123 | iw_u.parts.lsw = (ix1); \ |
| 124 | (d) = iw_u.value; \ |
| 125 | } while (0) |
| 126 | |
| 127 | /* Set a double from a 64-bit int. */ |
| 128 | #define INSERT_WORD64(d,ix) \ |
| 129 | do { \ |
| 130 | ieee_double_shape_type iw_u; \ |
| 131 | iw_u.xparts.w = (ix); \ |
| 132 | (d) = iw_u.value; \ |
| 133 | } while (0) |
| 134 | |
| 135 | /* Set the more significant 32 bits of a double from an int. */ |
| 136 | |
| 137 | #define SET_HIGH_WORD(d,v) \ |
| 138 | do { \ |
| 139 | ieee_double_shape_type sh_u; \ |
| 140 | sh_u.value = (d); \ |
| 141 | sh_u.parts.msw = (v); \ |
| 142 | (d) = sh_u.value; \ |
| 143 | } while (0) |
| 144 | |
| 145 | /* Set the less significant 32 bits of a double from an int. */ |
| 146 | |
| 147 | #define SET_LOW_WORD(d,v) \ |
| 148 | do { \ |
| 149 | ieee_double_shape_type sl_u; \ |
| 150 | sl_u.value = (d); \ |
| 151 | sl_u.parts.lsw = (v); \ |
| 152 | (d) = sl_u.value; \ |
| 153 | } while (0) |
| 154 | |
| 155 | /* |
| 156 | * A union which permits us to convert between a float and a 32 bit |
| 157 | * int. |
| 158 | */ |
| 159 | |
| 160 | typedef union { |
| 161 | float value; |
| 162 | /* FIXME: Assumes 32 bit int. */ |
| 163 | unsigned int word; |
| 164 | } ieee_float_shape_type; |
| 165 | |
| 166 | /* Get a 32 bit int from a float. */ |
| 167 | |
| 168 | #define GET_FLOAT_WORD(i,d) \ |
| 169 | do { \ |
| 170 | ieee_float_shape_type gf_u; \ |
| 171 | gf_u.value = (d); \ |
| 172 | (i) = gf_u.word; \ |
| 173 | } while (0) |
| 174 | |
| 175 | /* Set a float from a 32 bit int. */ |
| 176 | |
| 177 | #define SET_FLOAT_WORD(d,i) \ |
| 178 | do { \ |
| 179 | ieee_float_shape_type sf_u; \ |
| 180 | sf_u.word = (i); \ |
| 181 | (d) = sf_u.value; \ |
| 182 | } while (0) |
| 183 | |
| 184 | /* |
| 185 | * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long |
| 186 | * double. |
| 187 | */ |
| 188 | |
| 189 | #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \ |
| 190 | do { \ |
| 191 | union IEEEl2bits ew_u; \ |
| 192 | ew_u.e = (d); \ |
| 193 | (ix0) = ew_u.xbits.expsign; \ |
| 194 | (ix1) = ew_u.xbits.man; \ |
| 195 | } while (0) |
| 196 | |
| 197 | /* |
| 198 | * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit |
| 199 | * long double. |
| 200 | */ |
| 201 | |
| 202 | #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ |
| 203 | do { \ |
| 204 | union IEEEl2bits ew_u; \ |
| 205 | ew_u.e = (d); \ |
| 206 | (ix0) = ew_u.xbits.expsign; \ |
| 207 | (ix1) = ew_u.xbits.manh; \ |
| 208 | (ix2) = ew_u.xbits.manl; \ |
| 209 | } while (0) |
| 210 | |
| 211 | /* Get expsign as a 16 bit int from a long double. */ |
| 212 | |
| 213 | #define GET_LDBL_EXPSIGN(i,d) \ |
| 214 | do { \ |
| 215 | union IEEEl2bits ge_u; \ |
| 216 | ge_u.e = (d); \ |
| 217 | (i) = ge_u.xbits.expsign; \ |
| 218 | } while (0) |
| 219 | |
| 220 | /* |
| 221 | * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int |
| 222 | * mantissa. |
| 223 | */ |
| 224 | |
| 225 | #define INSERT_LDBL80_WORDS(d,ix0,ix1) \ |
| 226 | do { \ |
| 227 | union IEEEl2bits iw_u; \ |
| 228 | iw_u.xbits.expsign = (ix0); \ |
| 229 | iw_u.xbits.man = (ix1); \ |
| 230 | (d) = iw_u.e; \ |
| 231 | } while (0) |
| 232 | |
| 233 | /* |
| 234 | * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints |
| 235 | * comprising the mantissa. |
| 236 | */ |
| 237 | |
| 238 | #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \ |
| 239 | do { \ |
| 240 | union IEEEl2bits iw_u; \ |
| 241 | iw_u.xbits.expsign = (ix0); \ |
| 242 | iw_u.xbits.manh = (ix1); \ |
| 243 | iw_u.xbits.manl = (ix2); \ |
| 244 | (d) = iw_u.e; \ |
| 245 | } while (0) |
| 246 | |
| 247 | /* Set expsign of a long double from a 16 bit int. */ |
| 248 | |
| 249 | #define SET_LDBL_EXPSIGN(d,v) \ |
| 250 | do { \ |
| 251 | union IEEEl2bits se_u; \ |
| 252 | se_u.e = (d); \ |
| 253 | se_u.xbits.expsign = (v); \ |
| 254 | (d) = se_u.e; \ |
| 255 | } while (0) |
| 256 | |
| 257 | #ifdef __i386__ |
| 258 | /* Long double constants are broken on i386. */ |
| 259 | #define LD80C(m, ex, v) { \ |
| 260 | .xbits.man = __CONCAT(m, ULL), \ |
| 261 | .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \ |
| 262 | } |
| 263 | #else |
| 264 | /* The above works on non-i386 too, but we use this to check v. */ |
| 265 | #define LD80C(m, ex, v) { .e = (v), } |
| 266 | #endif |
| 267 | |
| 268 | #ifdef FLT_EVAL_METHOD |
| 269 | /* |
| 270 | * Attempt to get strict C99 semantics for assignment with non-C99 compilers. |
| 271 | */ |
| 272 | #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 |
| 273 | #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) |
| 274 | #else |
| 275 | #define STRICT_ASSIGN(type, lval, rval) do { \ |
| 276 | volatile type __lval; \ |
| 277 | \ |
| 278 | if (sizeof(type) >= sizeof(long double)) \ |
| 279 | (lval) = (rval); \ |
| 280 | else { \ |
| 281 | __lval = (rval); \ |
| 282 | (lval) = __lval; \ |
| 283 | } \ |
| 284 | } while (0) |
| 285 | #endif |
| 286 | #endif /* FLT_EVAL_METHOD */ |
| 287 | |
| 288 | /* Support switching the mode to FP_PE if necessary. */ |
| 289 | #if defined(__i386__) && !defined(NO_FPSETPREC) |
| 290 | #define ENTERI() \ |
| 291 | long double __retval; \ |
| 292 | fp_prec_t __oprec; \ |
| 293 | \ |
| 294 | if ((__oprec = fpgetprec()) != FP_PE) \ |
| 295 | fpsetprec(FP_PE) |
| 296 | #define RETURNI(x) do { \ |
| 297 | __retval = (x); \ |
| 298 | if (__oprec != FP_PE) \ |
| 299 | fpsetprec(__oprec); \ |
| 300 | RETURNF(__retval); \ |
| 301 | } while (0) |
| 302 | #else |
| 303 | #define ENTERI(x) |
| 304 | #define RETURNI(x) RETURNF(x) |
| 305 | #endif |
| 306 | |
| 307 | /* Default return statement if hack*_t() is not used. */ |
| 308 | #define RETURNF(v) return (v) |
| 309 | |
| 310 | /* |
| 311 | * 2sum gives the same result as 2sumF without requiring |a| >= |b| or |
| 312 | * a == 0, but is slower. |
| 313 | */ |
| 314 | #define _2sum(a, b) do { \ |
| 315 | __typeof(a) __s, __w; \ |
| 316 | \ |
| 317 | __w = (a) + (b); \ |
| 318 | __s = __w - (a); \ |
| 319 | (b) = ((a) - (__w - __s)) + ((b) - __s); \ |
| 320 | (a) = __w; \ |
| 321 | } while (0) |
| 322 | |
| 323 | /* |
| 324 | * 2sumF algorithm. |
| 325 | * |
| 326 | * "Normalize" the terms in the infinite-precision expression a + b for |
| 327 | * the sum of 2 floating point values so that b is as small as possible |
| 328 | * relative to 'a'. (The resulting 'a' is the value of the expression in |
| 329 | * the same precision as 'a' and the resulting b is the rounding error.) |
| 330 | * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and |
| 331 | * exponent overflow or underflow must not occur. This uses a Theorem of |
| 332 | * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum" |
| 333 | * is apparently due to Skewchuk (1997). |
| 334 | * |
| 335 | * For this to always work, assignment of a + b to 'a' must not retain any |
| 336 | * extra precision in a + b. This is required by C standards but broken |
| 337 | * in many compilers. The brokenness cannot be worked around using |
| 338 | * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this |
| 339 | * algorithm would be destroyed by non-null strict assignments. (The |
| 340 | * compilers are correct to be broken -- the efficiency of all floating |
| 341 | * point code calculations would be destroyed similarly if they forced the |
| 342 | * conversions.) |
| 343 | * |
| 344 | * Fortunately, a case that works well can usually be arranged by building |
| 345 | * any extra precision into the type of 'a' -- 'a' should have type float_t, |
| 346 | * double_t or long double. b's type should be no larger than 'a's type. |
| 347 | * Callers should use these types with scopes as large as possible, to |
| 348 | * reduce their own extra-precision and efficiciency problems. In |
| 349 | * particular, they shouldn't convert back and forth just to call here. |
| 350 | */ |
| 351 | #ifdef DEBUG |
| 352 | #define _2sumF(a, b) do { \ |
| 353 | __typeof(a) __w; \ |
| 354 | volatile __typeof(a) __ia, __ib, __r, __vw; \ |
| 355 | \ |
| 356 | __ia = (a); \ |
| 357 | __ib = (b); \ |
| 358 | assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \ |
| 359 | \ |
| 360 | __w = (a) + (b); \ |
| 361 | (b) = ((a) - __w) + (b); \ |
| 362 | (a) = __w; \ |
| 363 | \ |
| 364 | /* The next 2 assertions are weak if (a) is already long double. */ \ |
| 365 | assert((long double)__ia + __ib == (long double)(a) + (b)); \ |
| 366 | __vw = __ia + __ib; \ |
| 367 | __r = __ia - __vw; \ |
| 368 | __r += __ib; \ |
| 369 | assert(__vw == (a) && __r == (b)); \ |
| 370 | } while (0) |
| 371 | #else /* !DEBUG */ |
| 372 | #define _2sumF(a, b) do { \ |
| 373 | __typeof(a) __w; \ |
| 374 | \ |
| 375 | __w = (a) + (b); \ |
| 376 | (b) = ((a) - __w) + (b); \ |
| 377 | (a) = __w; \ |
| 378 | } while (0) |
| 379 | #endif /* DEBUG */ |
| 380 | |
| 381 | /* |
| 382 | * Set x += c, where x is represented in extra precision as a + b. |
| 383 | * x must be sufficiently normalized and sufficiently larger than c, |
| 384 | * and the result is then sufficiently normalized. |
| 385 | * |
| 386 | * The details of ordering are that |a| must be >= |c| (so that (a, c) |
| 387 | * can be normalized without extra work to swap 'a' with c). The details of |
| 388 | * the normalization are that b must be small relative to the normalized 'a'. |
| 389 | * Normalization of (a, c) makes the normalized c tiny relative to the |
| 390 | * normalized a, so b remains small relative to 'a' in the result. However, |
| 391 | * b need not ever be tiny relative to 'a'. For example, b might be about |
| 392 | * 2**20 times smaller than 'a' to give about 20 extra bits of precision. |
| 393 | * That is usually enough, and adding c (which by normalization is about |
| 394 | * 2**53 times smaller than a) cannot change b significantly. However, |
| 395 | * cancellation of 'a' with c in normalization of (a, c) may reduce 'a' |
| 396 | * significantly relative to b. The caller must ensure that significant |
| 397 | * cancellation doesn't occur, either by having c of the same sign as 'a', |
| 398 | * or by having |c| a few percent smaller than |a|. Pre-normalization of |
| 399 | * (a, b) may help. |
| 400 | * |
| 401 | * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2 |
| 402 | * exercise 19). We gain considerable efficiency by requiring the terms to |
| 403 | * be sufficiently normalized and sufficiently increasing. |
| 404 | */ |
| 405 | #define _3sumF(a, b, c) do { \ |
| 406 | __typeof(a) __tmp; \ |
| 407 | \ |
| 408 | __tmp = (c); \ |
| 409 | _2sumF(__tmp, (a)); \ |
| 410 | (b) += (a); \ |
| 411 | (a) = __tmp; \ |
| 412 | } while (0) |
| 413 | |
| 414 | /* |
| 415 | * Common routine to process the arguments to nan(), nanf(), and nanl(). |
| 416 | */ |
| 417 | void _scan_nan(uint32_t *__words, int __num_words, const char *__s); |
| 418 | |
| 419 | #ifdef _COMPLEX_H |
| 420 | |
| 421 | /* |
| 422 | * C99 specifies that complex numbers have the same representation as |
| 423 | * an array of two elements, where the first element is the real part |
| 424 | * and the second element is the imaginary part. |
| 425 | */ |
| 426 | typedef union { |
| 427 | float complex f; |
| 428 | float a[2]; |
| 429 | } float_complex; |
| 430 | typedef union { |
| 431 | double complex f; |
| 432 | double a[2]; |
| 433 | } double_complex; |
| 434 | typedef union { |
| 435 | long double complex f; |
| 436 | long double a[2]; |
| 437 | } long_double_complex; |
| 438 | #define REALPART(z) ((z).a[0]) |
| 439 | #define IMAGPART(z) ((z).a[1]) |
| 440 | |
| 441 | /* |
| 442 | * Inline functions that can be used to construct complex values. |
| 443 | * |
| 444 | * The C99 standard intends x+I*y to be used for this, but x+I*y is |
| 445 | * currently unusable in general since gcc introduces many overflow, |
| 446 | * underflow, sign and efficiency bugs by rewriting I*y as |
| 447 | * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. |
| 448 | * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted |
| 449 | * to -0.0+I*0.0. |
| 450 | */ |
| 451 | static __inline float complex |
| 452 | cpackf(float x, float y) |
| 453 | { |
| 454 | float_complex z; |
| 455 | |
| 456 | REALPART(z) = x; |
| 457 | IMAGPART(z) = y; |
| 458 | return (z.f); |
| 459 | } |
| 460 | |
| 461 | static __inline double complex |
| 462 | cpack(double x, double y) |
| 463 | { |
| 464 | double_complex z; |
| 465 | |
| 466 | REALPART(z) = x; |
| 467 | IMAGPART(z) = y; |
| 468 | return (z.f); |
| 469 | } |
| 470 | |
| 471 | static __inline long double complex |
| 472 | cpackl(long double x, long double y) |
| 473 | { |
| 474 | long_double_complex z; |
| 475 | |
| 476 | REALPART(z) = x; |
| 477 | IMAGPART(z) = y; |
| 478 | return (z.f); |
| 479 | } |
| 480 | #endif /* _COMPLEX_H */ |
| 481 | |
| 482 | #ifdef __GNUCLIKE_ASM |
| 483 | |
| 484 | /* Asm versions of some functions. */ |
| 485 | |
| 486 | #ifdef __amd64__ |
| 487 | static __inline int |
| 488 | irint(double x) |
| 489 | { |
| 490 | int n; |
| 491 | |
| 492 | asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); |
| 493 | return (n); |
| 494 | } |
| 495 | #define HAVE_EFFICIENT_IRINT |
| 496 | #endif |
| 497 | |
| 498 | #ifdef __i386__ |
| 499 | static __inline int |
| 500 | irint(double x) |
| 501 | { |
| 502 | int n; |
| 503 | |
| 504 | asm("fistl %0" : "=m" (n) : "t" (x)); |
| 505 | return (n); |
| 506 | } |
| 507 | #define HAVE_EFFICIENT_IRINT |
| 508 | #endif |
| 509 | |
| 510 | #if defined(__amd64__) || defined(__i386__) |
| 511 | static __inline int |
| 512 | irintl(long double x) |
| 513 | { |
| 514 | int n; |
| 515 | |
| 516 | asm("fistl %0" : "=m" (n) : "t" (x)); |
| 517 | return (n); |
| 518 | } |
| 519 | #define HAVE_EFFICIENT_IRINTL |
| 520 | #endif |
| 521 | |
| 522 | #endif /* __GNUCLIKE_ASM */ |
| 523 | |
| 524 | #ifdef DEBUG |
| 525 | #if defined(__amd64__) || defined(__i386__) |
| 526 | #define breakpoint() asm("int $3") |
| 527 | #else |
| 528 | #include <signal.h> |
| 529 | |
| 530 | #define breakpoint() raise(SIGTRAP) |
| 531 | #endif |
| 532 | #endif |
| 533 | |
| 534 | /* Write a pari script to test things externally. */ |
| 535 | #ifdef DOPRINT |
| 536 | #include <stdio.h> |
| 537 | |
| 538 | #ifndef DOPRINT_SWIZZLE |
| 539 | #define DOPRINT_SWIZZLE 0 |
| 540 | #endif |
| 541 | |
| 542 | #ifdef DOPRINT_LD80 |
| 543 | |
| 544 | #define DOPRINT_START(xp) do { \ |
| 545 | uint64_t __lx; \ |
| 546 | uint16_t __hx; \ |
| 547 | \ |
| 548 | /* Hack to give more-problematic args. */ \ |
| 549 | EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \ |
| 550 | __lx ^= DOPRINT_SWIZZLE; \ |
| 551 | INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ |
| 552 | printf("x = %.21Lg; ", (long double)*xp); \ |
| 553 | } while (0) |
| 554 | #define DOPRINT_END1(v) \ |
| 555 | printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) |
| 556 | #define DOPRINT_END2(hi, lo) \ |
| 557 | printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ |
| 558 | (long double)(hi), (long double)(lo)) |
| 559 | |
| 560 | #elif defined(DOPRINT_D64) |
| 561 | |
| 562 | #define DOPRINT_START(xp) do { \ |
| 563 | uint32_t __hx, __lx; \ |
| 564 | \ |
| 565 | EXTRACT_WORDS(__hx, __lx, *xp); \ |
| 566 | __lx ^= DOPRINT_SWIZZLE; \ |
| 567 | INSERT_WORDS(*xp, __hx, __lx); \ |
| 568 | printf("x = %.21Lg; ", (long double)*xp); \ |
| 569 | } while (0) |
| 570 | #define DOPRINT_END1(v) \ |
| 571 | printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) |
| 572 | #define DOPRINT_END2(hi, lo) \ |
| 573 | printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ |
| 574 | (long double)(hi), (long double)(lo)) |
| 575 | |
| 576 | #elif defined(DOPRINT_F32) |
| 577 | |
| 578 | #define DOPRINT_START(xp) do { \ |
| 579 | uint32_t __hx; \ |
| 580 | \ |
| 581 | GET_FLOAT_WORD(__hx, *xp); \ |
| 582 | __hx ^= DOPRINT_SWIZZLE; \ |
| 583 | SET_FLOAT_WORD(*xp, __hx); \ |
| 584 | printf("x = %.21Lg; ", (long double)*xp); \ |
| 585 | } while (0) |
| 586 | #define DOPRINT_END1(v) \ |
| 587 | printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) |
| 588 | #define DOPRINT_END2(hi, lo) \ |
| 589 | printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ |
| 590 | (long double)(hi), (long double)(lo)) |
| 591 | |
| 592 | #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */ |
| 593 | |
| 594 | #ifndef DOPRINT_SWIZZLE_HIGH |
| 595 | #define DOPRINT_SWIZZLE_HIGH 0 |
| 596 | #endif |
| 597 | |
| 598 | #define DOPRINT_START(xp) do { \ |
| 599 | uint64_t __lx, __llx; \ |
| 600 | uint16_t __hx; \ |
| 601 | \ |
| 602 | EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \ |
| 603 | __llx ^= DOPRINT_SWIZZLE; \ |
| 604 | __lx ^= DOPRINT_SWIZZLE_HIGH; \ |
| 605 | INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ |
| 606 | printf("x = %.36Lg; ", (long double)*xp); \ |
| 607 | } while (0) |
| 608 | #define DOPRINT_END1(v) \ |
| 609 | printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v)) |
| 610 | #define DOPRINT_END2(hi, lo) \ |
| 611 | printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \ |
| 612 | (long double)(hi), (long double)(lo)) |
| 613 | |
| 614 | #endif /* DOPRINT_LD80 */ |
| 615 | |
| 616 | #else /* !DOPRINT */ |
| 617 | #define DOPRINT_START(xp) |
| 618 | #define DOPRINT_END1(v) |
| 619 | #define DOPRINT_END2(hi, lo) |
| 620 | #endif /* DOPRINT */ |
| 621 | |
| 622 | #define RETURNP(x) do { \ |
| 623 | DOPRINT_END1(x); \ |
| 624 | RETURNF(x); \ |
| 625 | } while (0) |
| 626 | #define RETURNPI(x) do { \ |
| 627 | DOPRINT_END1(x); \ |
| 628 | RETURNI(x); \ |
| 629 | } while (0) |
| 630 | #define RETURN2P(x, y) do { \ |
| 631 | DOPRINT_END2((x), (y)); \ |
| 632 | RETURNF((x) + (y)); \ |
| 633 | } while (0) |
| 634 | #define RETURN2PI(x, y) do { \ |
| 635 | DOPRINT_END2((x), (y)); \ |
| 636 | RETURNI((x) + (y)); \ |
| 637 | } while (0) |
| 638 | #ifdef STRUCT_RETURN |
| 639 | #define RETURNSP(rp) do { \ |
| 640 | if (!(rp)->lo_set) \ |
| 641 | RETURNP((rp)->hi); \ |
| 642 | RETURN2P((rp)->hi, (rp)->lo); \ |
| 643 | } while (0) |
| 644 | #define RETURNSPI(rp) do { \ |
| 645 | if (!(rp)->lo_set) \ |
| 646 | RETURNPI((rp)->hi); \ |
| 647 | RETURN2PI((rp)->hi, (rp)->lo); \ |
| 648 | } while (0) |
| 649 | #endif |
| 650 | #define SUM2P(x, y) ({ \ |
| 651 | const __typeof (x) __x = (x); \ |
| 652 | const __typeof (y) __y = (y); \ |
| 653 | \ |
| 654 | DOPRINT_END2(__x, __y); \ |
| 655 | __x + __y; \ |
| 656 | }) |
| 657 | |
| 658 | /* |
| 659 | * ieee style elementary functions |
| 660 | * |
| 661 | * We rename functions here to improve other sources' diffability |
| 662 | * against fdlibm. |
| 663 | */ |
| 664 | #define __ieee754_sqrt sqrt |
| 665 | #define __ieee754_acos acos |
| 666 | #define __ieee754_acosh acosh |
| 667 | #define __ieee754_log log |
| 668 | #define __ieee754_log2 log2 |
| 669 | #define __ieee754_atanh atanh |
| 670 | #define __ieee754_asin asin |
| 671 | #define __ieee754_atan2 atan2 |
| 672 | #define __ieee754_exp exp |
| 673 | #define __ieee754_cosh cosh |
| 674 | #define __ieee754_fmod fmod |
| 675 | #define __ieee754_pow pow |
| 676 | #define __ieee754_lgamma lgamma |
| 677 | #define __ieee754_gamma gamma |
| 678 | #define __ieee754_lgamma_r lgamma_r |
| 679 | #define __ieee754_gamma_r gamma_r |
| 680 | #define __ieee754_log10 log10 |
| 681 | #define __ieee754_sinh sinh |
| 682 | #define __ieee754_hypot hypot |
| 683 | #define __ieee754_j0 j0 |
| 684 | #define __ieee754_j1 j1 |
| 685 | #define __ieee754_y0 y0 |
| 686 | #define __ieee754_y1 y1 |
| 687 | #define __ieee754_jn jn |
| 688 | #define __ieee754_yn yn |
| 689 | #define __ieee754_remainder remainder |
| 690 | #define __ieee754_scalb scalb |
| 691 | #define __ieee754_sqrtf sqrtf |
| 692 | #define __ieee754_acosf acosf |
| 693 | #define __ieee754_acoshf acoshf |
| 694 | #define __ieee754_logf logf |
| 695 | #define __ieee754_atanhf atanhf |
| 696 | #define __ieee754_asinf asinf |
| 697 | #define __ieee754_atan2f atan2f |
| 698 | #define __ieee754_expf expf |
| 699 | #define __ieee754_coshf coshf |
| 700 | #define __ieee754_fmodf fmodf |
| 701 | #define __ieee754_powf powf |
| 702 | #define __ieee754_lgammaf lgammaf |
| 703 | #define __ieee754_gammaf gammaf |
| 704 | #define __ieee754_lgammaf_r lgammaf_r |
| 705 | #define __ieee754_gammaf_r gammaf_r |
| 706 | #define __ieee754_log10f log10f |
| 707 | #define __ieee754_log2f log2f |
| 708 | #define __ieee754_sinhf sinhf |
| 709 | #define __ieee754_hypotf hypotf |
| 710 | #define __ieee754_j0f j0f |
| 711 | #define __ieee754_j1f j1f |
| 712 | #define __ieee754_y0f y0f |
| 713 | #define __ieee754_y1f y1f |
| 714 | #define __ieee754_jnf jnf |
| 715 | #define __ieee754_ynf ynf |
| 716 | #define __ieee754_remainderf remainderf |
| 717 | #define __ieee754_scalbf scalbf |
| 718 | |
| 719 | /* fdlibm kernel function */ |
| 720 | int __kernel_rem_pio2(double*,double*,int,int,int); |
| 721 | |
| 722 | /* double precision kernel functions */ |
| 723 | #ifndef INLINE_REM_PIO2 |
| 724 | int __ieee754_rem_pio2(double,double*); |
| 725 | #endif |
| 726 | double __kernel_sin(double,double,int); |
| 727 | double __kernel_cos(double,double); |
| 728 | double __kernel_tan(double,double,int); |
| 729 | double __ldexp_exp(double,int); |
| 730 | #ifdef _COMPLEX_H |
| 731 | double complex __ldexp_cexp(double complex,int); |
| 732 | #endif |
| 733 | |
| 734 | /* float precision kernel functions */ |
| 735 | #ifndef INLINE_REM_PIO2F |
| 736 | int __ieee754_rem_pio2f(float,double*); |
| 737 | #endif |
| 738 | #ifndef INLINE_KERNEL_SINDF |
| 739 | float __kernel_sindf(double); |
| 740 | #endif |
| 741 | #ifndef INLINE_KERNEL_COSDF |
| 742 | float __kernel_cosdf(double); |
| 743 | #endif |
| 744 | #ifndef INLINE_KERNEL_TANDF |
| 745 | float __kernel_tandf(double,int); |
| 746 | #endif |
| 747 | float __ldexp_expf(float,int); |
| 748 | #ifdef _COMPLEX_H |
| 749 | float complex __ldexp_cexpf(float complex,int); |
| 750 | #endif |
| 751 | |
| 752 | /* long double precision kernel functions */ |
| 753 | long double __kernel_sinl(long double, long double, int); |
| 754 | long double __kernel_cosl(long double, long double); |
| 755 | long double __kernel_tanl(long double, long double, int); |
| 756 | |
| 757 | #endif /* !_MATH_PRIVATE_H_ */ |