Index: head/lib/msun/src/catrigl.c =================================================================== --- head/lib/msun/src/catrigl.c (revision 322434) +++ head/lib/msun/src/catrigl.c (revision 322435) @@ -1,412 +1,412 @@ /*- * Copyright (c) 2012 Stephen Montgomery-Smith * Copyright (c) 2017 Mahdi Mokhtari * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ /* * The algorithm is very close to that in "Implementing the complex arcsine * and arccosine functions using exception handling" by T. E. Hull, Thomas F. * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335, * http://dl.acm.org/citation.cfm?id=275324. * * See catrig.c for complete comments. * * XXX comments were removed automatically, and even short ones on the right * of statements were removed (all of them), contrary to normal style. Only * a few comments on the right of declarations remain. */ #include __FBSDID("$FreeBSD$"); #include #include #include "invtrig.h" #include "math.h" #include "math_private.h" #undef isinf #define isinf(x) (fabsl(x) == INFINITY) #undef isnan #define isnan(x) ((x) != (x)) #define raise_inexact() do { volatile float junk __unused = 1 + tiny; } while(0) #undef signbit #define signbit(x) (__builtin_signbitl(x)) static const long double A_crossover = 10, B_crossover = 0.6417, FOUR_SQRT_MIN = 0x1p-8189L, QUARTER_SQRT_MAX = 0x1p8189L, RECIP_EPSILON = 1 / LDBL_EPSILON, SQRT_MIN = 0x1p-8191L; #if LDBL_MANT_DIG == 64 static const union IEEEl2bits um_e = LD80C(0xadf85458a2bb4a9b, 1, 2.71828182845904523536e+0L), um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); #define m_e um_e.e #define m_ln2 um_ln2.e static const long double /* The next 2 literals for non-i386. Misrounding them on i386 is harmless. */ SQRT_3_EPSILON = 5.70316273435758915310e-10, /* 0x9cc470a0490973e8.0p-94 */ SQRT_6_EPSILON = 8.06549008734932771664e-10; /* 0xddb3d742c265539e.0p-94 */ #elif LDBL_MANT_DIG == 113 static const long double m_e = 2.71828182845904523536028747135266250e0L, /* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */ m_ln2 = 6.93147180559945309417232121458176568e-1L, /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17, /* 0x1bb67ae8584caa73b25742d7078b8.0p-168 */ SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17; /* 0x13988e1409212e7d0321914321a55.0p-167 */ #else #error "Unsupported long double format" #endif static const volatile float tiny = 0x1p-100; static long double complex clog_for_large_values(long double complex z); static inline long double f(long double a, long double b, long double hypot_a_b) { if (b < 0) return ((hypot_a_b - b) / 2); if (b == 0) return (a / 2); return (a * a / (hypot_a_b + b) / 2); } static inline void do_hard_work(long double x, long double y, long double *rx, int *B_is_usable, long double *B, long double *sqrt_A2my2, long double *new_y) { long double R, S, A; long double Am1, Amy; R = hypotl(x, y + 1); S = hypotl(x, y - 1); A = (R + S) / 2; if (A < 1) A = 1; if (A < A_crossover) { if (y == 1 && x < LDBL_EPSILON * LDBL_EPSILON / 128) { *rx = sqrtl(x); } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); *rx = log1pl(Am1 + sqrtl(Am1 * (A + 1))); } else if (y < 1) { *rx = x / sqrtl((1 - y) * (1 + y)); } else { *rx = log1pl((y - 1) + sqrtl((y - 1) * (y + 1))); } } else { *rx = logl(A + sqrtl(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { *B_is_usable = 0; *sqrt_A2my2 = A * (2 / LDBL_EPSILON); *new_y = y * (2 / LDBL_EPSILON); return; } *B = y / A; *B_is_usable = 1; if (*B > B_crossover) { *B_is_usable = 0; if (y == 1 && x < LDBL_EPSILON / 128) { *sqrt_A2my2 = sqrtl(x) * sqrtl((A + y) / 2); } else if (x >= LDBL_EPSILON * fabsl(y - 1)) { Amy = f(x, y + 1, R) + f(x, y - 1, S); *sqrt_A2my2 = sqrtl(Amy * (A + y)); } else if (y > 1) { *sqrt_A2my2 = x * (4 / LDBL_EPSILON / LDBL_EPSILON) * y / sqrtl((y + 1) * (y - 1)); *new_y = y * (4 / LDBL_EPSILON / LDBL_EPSILON); } else { *sqrt_A2my2 = sqrtl((1 - y) * (1 + y)); } } } long double complex casinhl(long double complex z) { long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; long double complex w; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(x, y + y)); if (isinf(y)) return (CMPLXL(y, x + x)); if (y == 0) return (CMPLXL(x + x, y)); return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { if (signbit(x) == 0) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; return (CMPLXL(copysignl(creall(w), x), copysignl(cimagl(w), y))); } if (x == 0 && y == 0) return (z); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (z); do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y); if (B_is_usable) ry = asinl(B); else ry = atan2l(new_y, sqrt_A2my2); return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } long double complex casinl(long double complex z) { long double complex w; w = casinhl(CMPLXL(cimagl(z), creall(z))); return (CMPLXL(cimagl(w), creall(w))); } long double complex cacosl(long double complex z) { long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; long double complex w; x = creall(z); y = cimagl(z); sx = signbit(x); sy = signbit(y); ax = fabsl(x); ay = fabsl(y); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(y + y, -INFINITY)); if (isinf(y)) return (CMPLXL(x + x, -y)); if (x == 0) return (CMPLXL(pio2_hi + pio2_lo, y + y)); return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { w = clog_for_large_values(z); rx = fabsl(cimagl(w)); ry = creall(w) + m_ln2; if (sy == 0) ry = -ry; return (CMPLXL(rx, ry)); } if (x == 1 && y == 0) return (CMPLXL(0, -y)); raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (CMPLXL(pio2_hi - (x - pio2_lo), -y)); do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x); if (B_is_usable) { if (sx == 0) rx = acosl(B); else rx = acosl(-B); } else { if (sx == 0) rx = atan2l(sqrt_A2mx2, new_x); else rx = atan2l(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; return (CMPLXL(rx, ry)); } long double complex cacoshl(long double complex z) { long double complex w; long double rx, ry; w = cacosl(z); rx = creall(w); ry = cimagl(w); if (isnan(rx) && isnan(ry)) return (CMPLXL(ry, rx)); if (isnan(rx)) return (CMPLXL(fabsl(ry), rx)); if (isnan(ry)) return (CMPLXL(ry, ry)); return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z)))); } static long double complex clog_for_large_values(long double complex z) { long double x, y; long double ax, ay, t; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (ax < ay) { t = ax; ax = ay; ay = t; } - if (ax >= HALF_LDBL_MAX) + if (ax > LDBL_MAX / 2) return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, atan2l(y, x))); if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x))); return (CMPLXL(logl(ax * ax + ay * ay) / 2, atan2l(y, x))); } static inline long double sum_squares(long double x, long double y) { if (y < SQRT_MIN) return (x * x); return (x * x + y * y); } static inline long double real_part_reciprocal(long double x, long double y) { long double scale; uint16_t hx, hy; int16_t ix, iy; GET_LDBL_EXPSIGN(hx, x); ix = hx & 0x7fff; GET_LDBL_EXPSIGN(hy, y); iy = hy & 0x7fff; #define BIAS (LDBL_MAX_EXP - 1) #define CUTOFF (LDBL_MANT_DIG / 2 + 1) if (ix - iy >= CUTOFF || isinf(x)) return (1 / x); if (iy - ix >= CUTOFF) return (x / y / y); if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF) return (x / (x * x + y * y)); scale = 1; SET_LDBL_EXPSIGN(scale, 0x7fff - ix); x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } long double complex catanhl(long double complex z) { long double x, y, ax, ay, rx, ry; x = creall(z); y = cimagl(z); ax = fabsl(x); ay = fabsl(y); if (y == 0 && ax <= 1) return (CMPLXL(atanhl(x), y)); if (x == 0) return (CMPLXL(x, atanl(y))); if (isnan(x) || isnan(y)) { if (isinf(x)) return (CMPLXL(copysignl(0, x), y + y)); if (isinf(y)) return (CMPLXL(copysignl(0, x), copysignl(pio2_hi + pio2_lo, y))); return (CMPLXL(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) return (CMPLXL(real_part_reciprocal(x, y), copysignl(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { raise_inexact(); return (z); } if (ax == 1 && ay < LDBL_EPSILON) rx = (m_ln2 - logl(ay)) / 2; else rx = log1pl(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) ry = atan2l(2, -ay) / 2; else if (ay < LDBL_EPSILON) ry = atan2l(2 * ay, (1 - ax) * (1 + ax)) / 2; else ry = atan2l(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; return (CMPLXL(copysignl(rx, x), copysignl(ry, y))); } long double complex catanl(long double complex z) { long double complex w; w = catanhl(CMPLXL(cimagl(z), creall(z))); return (CMPLXL(cimagl(w), creall(w))); } Index: head/lib/msun/src/math_private.h =================================================================== --- head/lib/msun/src/math_private.h (revision 322434) +++ head/lib/msun/src/math_private.h (revision 322435) @@ -1,797 +1,788 @@ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* * from: @(#)fdlibm.h 5.1 93/09/24 * $FreeBSD$ */ #ifndef _MATH_PRIVATE_H_ #define _MATH_PRIVATE_H_ #include #include /* * The original fdlibm code used statements like: * n0 = ((*(int*)&one)>>29)^1; * index of high word * * ix0 = *(n0+(int*)&x); * high word of x * * ix1 = *((1-n0)+(int*)&x); * low word of x * * to dig two 32 bit words out of the 64 bit IEEE floating point * value. That is non-ANSI, and, moreover, the gcc instruction * scheduler gets it wrong. We instead use the following macros. * Unlike the original code, we determine the endianness at compile * time, not at run time; I don't see much benefit to selecting * endianness at run time. */ /* * A union which permits us to convert between a double and two 32 bit * ints. */ #ifdef __arm__ #if defined(__VFP_FP__) || defined(__ARM_EABI__) #define IEEE_WORD_ORDER BYTE_ORDER #else #define IEEE_WORD_ORDER BIG_ENDIAN #endif #else /* __arm__ */ #define IEEE_WORD_ORDER BYTE_ORDER #endif #if IEEE_WORD_ORDER == BIG_ENDIAN typedef union { double value; struct { u_int32_t msw; u_int32_t lsw; } parts; struct { u_int64_t w; } xparts; } ieee_double_shape_type; #endif #if IEEE_WORD_ORDER == LITTLE_ENDIAN typedef union { double value; struct { u_int32_t lsw; u_int32_t msw; } parts; struct { u_int64_t w; } xparts; } ieee_double_shape_type; #endif /* Get two 32 bit ints from a double. */ #define EXTRACT_WORDS(ix0,ix1,d) \ do { \ ieee_double_shape_type ew_u; \ ew_u.value = (d); \ (ix0) = ew_u.parts.msw; \ (ix1) = ew_u.parts.lsw; \ } while (0) /* Get a 64-bit int from a double. */ #define EXTRACT_WORD64(ix,d) \ do { \ ieee_double_shape_type ew_u; \ ew_u.value = (d); \ (ix) = ew_u.xparts.w; \ } while (0) /* Get the more significant 32 bit int from a double. */ #define GET_HIGH_WORD(i,d) \ do { \ ieee_double_shape_type gh_u; \ gh_u.value = (d); \ (i) = gh_u.parts.msw; \ } while (0) /* Get the less significant 32 bit int from a double. */ #define GET_LOW_WORD(i,d) \ do { \ ieee_double_shape_type gl_u; \ gl_u.value = (d); \ (i) = gl_u.parts.lsw; \ } while (0) /* Set a double from two 32 bit ints. */ #define INSERT_WORDS(d,ix0,ix1) \ do { \ ieee_double_shape_type iw_u; \ iw_u.parts.msw = (ix0); \ iw_u.parts.lsw = (ix1); \ (d) = iw_u.value; \ } while (0) /* Set a double from a 64-bit int. */ #define INSERT_WORD64(d,ix) \ do { \ ieee_double_shape_type iw_u; \ iw_u.xparts.w = (ix); \ (d) = iw_u.value; \ } while (0) /* Set the more significant 32 bits of a double from an int. */ #define SET_HIGH_WORD(d,v) \ do { \ ieee_double_shape_type sh_u; \ sh_u.value = (d); \ sh_u.parts.msw = (v); \ (d) = sh_u.value; \ } while (0) /* Set the less significant 32 bits of a double from an int. */ #define SET_LOW_WORD(d,v) \ do { \ ieee_double_shape_type sl_u; \ sl_u.value = (d); \ sl_u.parts.lsw = (v); \ (d) = sl_u.value; \ } while (0) /* * A union which permits us to convert between a float and a 32 bit * int. */ typedef union { float value; /* FIXME: Assumes 32 bit int. */ unsigned int word; } ieee_float_shape_type; /* Get a 32 bit int from a float. */ #define GET_FLOAT_WORD(i,d) \ do { \ ieee_float_shape_type gf_u; \ gf_u.value = (d); \ (i) = gf_u.word; \ } while (0) /* Set a float from a 32 bit int. */ #define SET_FLOAT_WORD(d,i) \ do { \ ieee_float_shape_type sf_u; \ sf_u.word = (i); \ (d) = sf_u.value; \ } while (0) /* * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long * double. */ #define EXTRACT_LDBL80_WORDS(ix0,ix1,d) \ do { \ union IEEEl2bits ew_u; \ ew_u.e = (d); \ (ix0) = ew_u.xbits.expsign; \ (ix1) = ew_u.xbits.man; \ } while (0) /* * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit * long double. */ #define EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d) \ do { \ union IEEEl2bits ew_u; \ ew_u.e = (d); \ (ix0) = ew_u.xbits.expsign; \ (ix1) = ew_u.xbits.manh; \ (ix2) = ew_u.xbits.manl; \ } while (0) /* Get expsign as a 16 bit int from a long double. */ #define GET_LDBL_EXPSIGN(i,d) \ do { \ union IEEEl2bits ge_u; \ ge_u.e = (d); \ (i) = ge_u.xbits.expsign; \ } while (0) /* * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int * mantissa. */ #define INSERT_LDBL80_WORDS(d,ix0,ix1) \ do { \ union IEEEl2bits iw_u; \ iw_u.xbits.expsign = (ix0); \ iw_u.xbits.man = (ix1); \ (d) = iw_u.e; \ } while (0) /* * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints * comprising the mantissa. */ #define INSERT_LDBL128_WORDS(d,ix0,ix1,ix2) \ do { \ union IEEEl2bits iw_u; \ iw_u.xbits.expsign = (ix0); \ iw_u.xbits.manh = (ix1); \ iw_u.xbits.manl = (ix2); \ (d) = iw_u.e; \ } while (0) /* Set expsign of a long double from a 16 bit int. */ #define SET_LDBL_EXPSIGN(d,v) \ do { \ union IEEEl2bits se_u; \ se_u.e = (d); \ se_u.xbits.expsign = (v); \ (d) = se_u.e; \ } while (0) #ifdef __i386__ /* Long double constants are broken on i386. */ #define LD80C(m, ex, v) { \ .xbits.man = __CONCAT(m, ULL), \ .xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0), \ } #else /* The above works on non-i386 too, but we use this to check v. */ #define LD80C(m, ex, v) { .e = (v), } #endif -/* - * XXX LDBL_MAX is broken on i386. If the precise value of LDBL_MAX is not - * needed, this may be worked around by instead referring to a proxy, such - * as HALF_LDBL_MAX, below. HALF_LDBL_MAX is approximately LDBL_MAX / 2, - * actually just greater than. Note that 2 * HALF_LDBL_MAX will always - * overflow to infinity, regardless of the precision and rounding modes. - */ -#define HALF_LDBL_MAX __CONCAT(__CONCAT(0x0.8p, LDBL_MAX_EXP), L) - #ifdef FLT_EVAL_METHOD /* * Attempt to get strict C99 semantics for assignment with non-C99 compilers. */ #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0 #define STRICT_ASSIGN(type, lval, rval) ((lval) = (rval)) #else #define STRICT_ASSIGN(type, lval, rval) do { \ volatile type __lval; \ \ if (sizeof(type) >= sizeof(long double)) \ (lval) = (rval); \ else { \ __lval = (rval); \ (lval) = __lval; \ } \ } while (0) #endif #endif /* FLT_EVAL_METHOD */ /* Support switching the mode to FP_PE if necessary. */ #if defined(__i386__) && !defined(NO_FPSETPREC) #define ENTERI() \ long double __retval; \ fp_prec_t __oprec; \ \ if ((__oprec = fpgetprec()) != FP_PE) \ fpsetprec(FP_PE) #define RETURNI(x) do { \ __retval = (x); \ if (__oprec != FP_PE) \ fpsetprec(__oprec); \ RETURNF(__retval); \ } while (0) #define ENTERV() \ fp_prec_t __oprec; \ \ if ((__oprec = fpgetprec()) != FP_PE) \ fpsetprec(FP_PE) #define RETURNV() do { \ if (__oprec != FP_PE) \ fpsetprec(__oprec); \ return; \ } while (0) #else #define ENTERI() #define RETURNI(x) RETURNF(x) #define ENTERV() #define RETURNV() return #endif /* Default return statement if hack*_t() is not used. */ #define RETURNF(v) return (v) /* * 2sum gives the same result as 2sumF without requiring |a| >= |b| or * a == 0, but is slower. */ #define _2sum(a, b) do { \ __typeof(a) __s, __w; \ \ __w = (a) + (b); \ __s = __w - (a); \ (b) = ((a) - (__w - __s)) + ((b) - __s); \ (a) = __w; \ } while (0) /* * 2sumF algorithm. * * "Normalize" the terms in the infinite-precision expression a + b for * the sum of 2 floating point values so that b is as small as possible * relative to 'a'. (The resulting 'a' is the value of the expression in * the same precision as 'a' and the resulting b is the rounding error.) * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and * exponent overflow or underflow must not occur. This uses a Theorem of * Dekker (1971). See Knuth (1981) 4.2.2 Theorem C. The name "TwoSum" * is apparently due to Skewchuk (1997). * * For this to always work, assignment of a + b to 'a' must not retain any * extra precision in a + b. This is required by C standards but broken * in many compilers. The brokenness cannot be worked around using * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this * algorithm would be destroyed by non-null strict assignments. (The * compilers are correct to be broken -- the efficiency of all floating * point code calculations would be destroyed similarly if they forced the * conversions.) * * Fortunately, a case that works well can usually be arranged by building * any extra precision into the type of 'a' -- 'a' should have type float_t, * double_t or long double. b's type should be no larger than 'a's type. * Callers should use these types with scopes as large as possible, to * reduce their own extra-precision and efficiciency problems. In * particular, they shouldn't convert back and forth just to call here. */ #ifdef DEBUG #define _2sumF(a, b) do { \ __typeof(a) __w; \ volatile __typeof(a) __ia, __ib, __r, __vw; \ \ __ia = (a); \ __ib = (b); \ assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib)); \ \ __w = (a) + (b); \ (b) = ((a) - __w) + (b); \ (a) = __w; \ \ /* The next 2 assertions are weak if (a) is already long double. */ \ assert((long double)__ia + __ib == (long double)(a) + (b)); \ __vw = __ia + __ib; \ __r = __ia - __vw; \ __r += __ib; \ assert(__vw == (a) && __r == (b)); \ } while (0) #else /* !DEBUG */ #define _2sumF(a, b) do { \ __typeof(a) __w; \ \ __w = (a) + (b); \ (b) = ((a) - __w) + (b); \ (a) = __w; \ } while (0) #endif /* DEBUG */ /* * Set x += c, where x is represented in extra precision as a + b. * x must be sufficiently normalized and sufficiently larger than c, * and the result is then sufficiently normalized. * * The details of ordering are that |a| must be >= |c| (so that (a, c) * can be normalized without extra work to swap 'a' with c). The details of * the normalization are that b must be small relative to the normalized 'a'. * Normalization of (a, c) makes the normalized c tiny relative to the * normalized a, so b remains small relative to 'a' in the result. However, * b need not ever be tiny relative to 'a'. For example, b might be about * 2**20 times smaller than 'a' to give about 20 extra bits of precision. * That is usually enough, and adding c (which by normalization is about * 2**53 times smaller than a) cannot change b significantly. However, * cancellation of 'a' with c in normalization of (a, c) may reduce 'a' * significantly relative to b. The caller must ensure that significant * cancellation doesn't occur, either by having c of the same sign as 'a', * or by having |c| a few percent smaller than |a|. Pre-normalization of * (a, b) may help. * * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2 * exercise 19). We gain considerable efficiency by requiring the terms to * be sufficiently normalized and sufficiently increasing. */ #define _3sumF(a, b, c) do { \ __typeof(a) __tmp; \ \ __tmp = (c); \ _2sumF(__tmp, (a)); \ (b) += (a); \ (a) = __tmp; \ } while (0) /* * Common routine to process the arguments to nan(), nanf(), and nanl(). */ void _scan_nan(uint32_t *__words, int __num_words, const char *__s); #ifdef _COMPLEX_H /* * C99 specifies that complex numbers have the same representation as * an array of two elements, where the first element is the real part * and the second element is the imaginary part. */ typedef union { float complex f; float a[2]; } float_complex; typedef union { double complex f; double a[2]; } double_complex; typedef union { long double complex f; long double a[2]; } long_double_complex; #define REALPART(z) ((z).a[0]) #define IMAGPART(z) ((z).a[1]) /* * Inline functions that can be used to construct complex values. * * The C99 standard intends x+I*y to be used for this, but x+I*y is * currently unusable in general since gcc introduces many overflow, * underflow, sign and efficiency bugs by rewriting I*y as * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product. * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted * to -0.0+I*0.0. * * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL() * to construct complex values. Compilers that conform to the C99 * standard require the following functions to avoid the above issues. */ #ifndef CMPLXF static __inline float complex CMPLXF(float x, float y) { float_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #ifndef CMPLX static __inline double complex CMPLX(double x, double y) { double_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #ifndef CMPLXL static __inline long double complex CMPLXL(long double x, long double y) { long_double_complex z; REALPART(z) = x; IMAGPART(z) = y; return (z.f); } #endif #endif /* _COMPLEX_H */ #ifdef __GNUCLIKE_ASM /* Asm versions of some functions. */ #ifdef __amd64__ static __inline int irint(double x) { int n; asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); return (n); } #define HAVE_EFFICIENT_IRINT #endif #ifdef __i386__ static __inline int irint(double x) { int n; asm("fistl %0" : "=m" (n) : "t" (x)); return (n); } #define HAVE_EFFICIENT_IRINT #endif #if defined(__amd64__) || defined(__i386__) static __inline int irintl(long double x) { int n; asm("fistl %0" : "=m" (n) : "t" (x)); return (n); } #define HAVE_EFFICIENT_IRINTL #endif #endif /* __GNUCLIKE_ASM */ #ifdef DEBUG #if defined(__amd64__) || defined(__i386__) #define breakpoint() asm("int $3") #else #include #define breakpoint() raise(SIGTRAP) #endif #endif /* Write a pari script to test things externally. */ #ifdef DOPRINT #include #ifndef DOPRINT_SWIZZLE #define DOPRINT_SWIZZLE 0 #endif #ifdef DOPRINT_LD80 #define DOPRINT_START(xp) do { \ uint64_t __lx; \ uint16_t __hx; \ \ /* Hack to give more-problematic args. */ \ EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \ __lx ^= DOPRINT_SWIZZLE; \ INSERT_LDBL80_WORDS(*xp, __hx, __lx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #elif defined(DOPRINT_D64) #define DOPRINT_START(xp) do { \ uint32_t __hx, __lx; \ \ EXTRACT_WORDS(__hx, __lx, *xp); \ __lx ^= DOPRINT_SWIZZLE; \ INSERT_WORDS(*xp, __hx, __lx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #elif defined(DOPRINT_F32) #define DOPRINT_START(xp) do { \ uint32_t __hx; \ \ GET_FLOAT_WORD(__hx, *xp); \ __hx ^= DOPRINT_SWIZZLE; \ SET_FLOAT_WORD(*xp, __hx); \ printf("x = %.21Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */ #ifndef DOPRINT_SWIZZLE_HIGH #define DOPRINT_SWIZZLE_HIGH 0 #endif #define DOPRINT_START(xp) do { \ uint64_t __lx, __llx; \ uint16_t __hx; \ \ EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \ __llx ^= DOPRINT_SWIZZLE; \ __lx ^= DOPRINT_SWIZZLE_HIGH; \ INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \ printf("x = %.36Lg; ", (long double)*xp); \ } while (0) #define DOPRINT_END1(v) \ printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v)) #define DOPRINT_END2(hi, lo) \ printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \ (long double)(hi), (long double)(lo)) #endif /* DOPRINT_LD80 */ #else /* !DOPRINT */ #define DOPRINT_START(xp) #define DOPRINT_END1(v) #define DOPRINT_END2(hi, lo) #endif /* DOPRINT */ #define RETURNP(x) do { \ DOPRINT_END1(x); \ RETURNF(x); \ } while (0) #define RETURNPI(x) do { \ DOPRINT_END1(x); \ RETURNI(x); \ } while (0) #define RETURN2P(x, y) do { \ DOPRINT_END2((x), (y)); \ RETURNF((x) + (y)); \ } while (0) #define RETURN2PI(x, y) do { \ DOPRINT_END2((x), (y)); \ RETURNI((x) + (y)); \ } while (0) #ifdef STRUCT_RETURN #define RETURNSP(rp) do { \ if (!(rp)->lo_set) \ RETURNP((rp)->hi); \ RETURN2P((rp)->hi, (rp)->lo); \ } while (0) #define RETURNSPI(rp) do { \ if (!(rp)->lo_set) \ RETURNPI((rp)->hi); \ RETURN2PI((rp)->hi, (rp)->lo); \ } while (0) #endif #define SUM2P(x, y) ({ \ const __typeof (x) __x = (x); \ const __typeof (y) __y = (y); \ \ DOPRINT_END2(__x, __y); \ __x + __y; \ }) /* * ieee style elementary functions * * We rename functions here to improve other sources' diffability * against fdlibm. */ #define __ieee754_sqrt sqrt #define __ieee754_acos acos #define __ieee754_acosh acosh #define __ieee754_log log #define __ieee754_log2 log2 #define __ieee754_atanh atanh #define __ieee754_asin asin #define __ieee754_atan2 atan2 #define __ieee754_exp exp #define __ieee754_cosh cosh #define __ieee754_fmod fmod #define __ieee754_pow pow #define __ieee754_lgamma lgamma #define __ieee754_gamma gamma #define __ieee754_lgamma_r lgamma_r #define __ieee754_gamma_r gamma_r #define __ieee754_log10 log10 #define __ieee754_sinh sinh #define __ieee754_hypot hypot #define __ieee754_j0 j0 #define __ieee754_j1 j1 #define __ieee754_y0 y0 #define __ieee754_y1 y1 #define __ieee754_jn jn #define __ieee754_yn yn #define __ieee754_remainder remainder #define __ieee754_scalb scalb #define __ieee754_sqrtf sqrtf #define __ieee754_acosf acosf #define __ieee754_acoshf acoshf #define __ieee754_logf logf #define __ieee754_atanhf atanhf #define __ieee754_asinf asinf #define __ieee754_atan2f atan2f #define __ieee754_expf expf #define __ieee754_coshf coshf #define __ieee754_fmodf fmodf #define __ieee754_powf powf #define __ieee754_lgammaf lgammaf #define __ieee754_gammaf gammaf #define __ieee754_lgammaf_r lgammaf_r #define __ieee754_gammaf_r gammaf_r #define __ieee754_log10f log10f #define __ieee754_log2f log2f #define __ieee754_sinhf sinhf #define __ieee754_hypotf hypotf #define __ieee754_j0f j0f #define __ieee754_j1f j1f #define __ieee754_y0f y0f #define __ieee754_y1f y1f #define __ieee754_jnf jnf #define __ieee754_ynf ynf #define __ieee754_remainderf remainderf #define __ieee754_scalbf scalbf /* fdlibm kernel function */ int __kernel_rem_pio2(double*,double*,int,int,int); /* double precision kernel functions */ #ifndef INLINE_REM_PIO2 int __ieee754_rem_pio2(double,double*); #endif double __kernel_sin(double,double,int); double __kernel_cos(double,double); double __kernel_tan(double,double,int); double __ldexp_exp(double,int); #ifdef _COMPLEX_H double complex __ldexp_cexp(double complex,int); #endif /* float precision kernel functions */ #ifndef INLINE_REM_PIO2F int __ieee754_rem_pio2f(float,double*); #endif #ifndef INLINE_KERNEL_SINDF float __kernel_sindf(double); #endif #ifndef INLINE_KERNEL_COSDF float __kernel_cosdf(double); #endif #ifndef INLINE_KERNEL_TANDF float __kernel_tandf(double,int); #endif float __ldexp_expf(float,int); #ifdef _COMPLEX_H float complex __ldexp_cexpf(float complex,int); #endif /* long double precision kernel functions */ long double __kernel_sinl(long double, long double, int); long double __kernel_cosl(long double, long double); long double __kernel_tanl(long double, long double, int); #endif /* !_MATH_PRIVATE_H_ */ Index: head/lib/msun/src/s_csqrtl.c =================================================================== --- head/lib/msun/src/s_csqrtl.c (revision 322434) +++ head/lib/msun/src/s_csqrtl.c (revision 322435) @@ -1,108 +1,108 @@ /*- * Copyright (c) 2007-2008 David Schultz * All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include __FBSDID("$FreeBSD$"); #include #include #include #include "math_private.h" /* * gcc doesn't implement complex multiplication or division correctly, * so we need to handle infinities specially. We turn on this pragma to * notify conforming c99 compilers that the fast-but-incorrect code that * gcc generates is acceptable, since the special cases have already been * handled. */ #pragma STDC CX_LIMITED_RANGE ON /* We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). */ -#define THRESH (HALF_LDBL_MAX / 1.207106781186547524400844362104849L) +#define THRESH (LDBL_MAX / 2.414213562373095048801688724209698L) long double complex csqrtl(long double complex z) { long double complex result; long double a, b; long double t; int scale; a = creall(z); b = cimagl(z); /* Handle special cases. */ if (z == 0) return (CMPLXL(0, b)); if (isinf(b)) return (CMPLXL(INFINITY, b)); if (isnan(a)) { t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ return (CMPLXL(a, t)); /* return NaN + NaN i */ } if (isinf(a)) { /* * csqrt(inf + NaN i) = inf + NaN i * csqrt(inf + y i) = inf + 0 i * csqrt(-inf + NaN i) = NaN +- inf i * csqrt(-inf + y i) = 0 + inf i */ if (signbit(a)) return (CMPLXL(fabsl(b - b), copysignl(a, b))); else return (CMPLXL(a, copysignl(b - b, b))); } /* * The remaining special case (b is NaN) is handled just fine by * the normal code path below. */ /* Scale to avoid overflow. */ if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) { a *= 0.25; b *= 0.25; scale = 1; } else { scale = 0; } /* Algorithm 312, CACM vol 10, Oct 1967. */ if (a >= 0) { t = sqrtl((a + hypotl(a, b)) * 0.5); result = CMPLXL(t, b / (2 * t)); } else { t = sqrtl((-a + hypotl(a, b)) * 0.5); result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b)); } /* Rescale. */ if (scale) return (result * 2); else return (result); }