Index: stable/11/lib/msun/src/catrig.c =================================================================== --- stable/11/lib/msun/src/catrig.c (revision 323474) +++ stable/11/lib/msun/src/catrig.c (revision 323475) @@ -1,639 +1,644 @@ /*- * Copyright (c) 2012 Stephen Montgomery-Smith * 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 "math.h" #include "math_private.h" #undef isinf #define isinf(x) (fabs(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_signbit(x)) /* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */ static const double A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */ B_crossover = 0.6417, /* suggested by Hull et al */ FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */ QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */ m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */ m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */ pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */ RECIP_EPSILON = 1 / DBL_EPSILON, SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */ SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */ SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */ static const volatile double pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */ static const volatile float tiny = 0x1p-100; static double complex clog_for_large_values(double complex z); /* * Testing indicates that all these functions are accurate up to 4 ULP. * The functions casin(h) and cacos(h) are about 2.5 times slower than asinh. * The functions catan(h) are a little under 2 times slower than atanh. * * The code for casinh, casin, cacos, and cacosh comes first. The code is * rather complicated, and the four functions are highly interdependent. * * The code for catanh and catan comes at the end. It is much simpler than * the other functions, and the code for these can be disconnected from the * rest of the code. */ /* * ================================ * | casinh, casin, cacos, cacosh | * ================================ */ /* * 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. * * Throughout we use the convention z = x + I*y. * * casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B) * where * A = (|z+I| + |z-I|) / 2 * B = (|z+I| - |z-I|) / 2 = y/A * * These formulas become numerically unstable: * (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that * is, Re(casinh(z)) is close to 0); * (b) for Im(casinh(z)) when z is close to either of the intervals * [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is * close to PI/2). * * These numerical problems are overcome by defining * f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2 * Then if A < A_crossover, we use * log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1))) * A-1 = f(x, 1+y) + f(x, 1-y) * and if B > B_crossover, we use * asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y))) * A-y = f(x, y+1) + f(x, y-1) * where without loss of generality we have assumed that x and y are * non-negative. * * Much of the difficulty comes because the intermediate computations may * produce overflows or underflows. This is dealt with in the paper by Hull * et al by using exception handling. We do this by detecting when * computations risk underflow or overflow. The hardest part is handling the * underflows when computing f(a, b). * * Note that the function f(a, b) does not appear explicitly in the paper by * Hull et al, but the idea may be found on pages 308 and 309. Introducing the * function f(a, b) allows us to concentrate many of the clever tricks in this * paper into one function. */ /* * Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2. * Pass hypot(a, b) as the third argument. */ static inline double f(double a, double b, 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); } /* * All the hard work is contained in this function. * x and y are assumed positive or zero, and less than RECIP_EPSILON. * Upon return: * rx = Re(casinh(z)) = -Im(cacos(y + I*x)). * B_is_usable is set to 1 if the value of B is usable. * If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y. * If returning sqrt_A2my2 has potential to result in an underflow, it is * rescaled, and new_y is similarly rescaled. */ static inline void do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B, double *sqrt_A2my2, double *new_y) { double R, S, A; /* A, B, R, and S are as in Hull et al. */ double Am1, Amy; /* A-1, A-y. */ R = hypot(x, y + 1); /* |z+I| */ S = hypot(x, y - 1); /* |z-I| */ /* A = (|z+I| + |z-I|) / 2 */ A = (R + S) / 2; /* * Mathematically A >= 1. There is a small chance that this will not * be so because of rounding errors. So we will make certain it is * so. */ if (A < 1) A = 1; if (A < A_crossover) { /* * Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y). * rx = log1p(Am1 + sqrt(Am1*(A+1))) */ if (y == 1 && x < DBL_EPSILON * DBL_EPSILON / 128) { /* * fp is of order x^2, and fm = x/2. * A = 1 (inexactly). */ *rx = sqrt(x); } else if (x >= DBL_EPSILON * fabs(y - 1)) { /* * Underflow will not occur because * x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN */ Am1 = f(x, 1 + y, R) + f(x, 1 - y, S); *rx = log1p(Am1 + sqrt(Am1 * (A + 1))); } else if (y < 1) { /* * fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and * A = 1 (inexactly). */ *rx = x / sqrt((1 - y) * (1 + y)); } else { /* if (y > 1) */ /* * A-1 = y-1 (inexactly). */ *rx = log1p((y - 1) + sqrt((y - 1) * (y + 1))); } } else { *rx = log(A + sqrt(A * A - 1)); } *new_y = y; if (y < FOUR_SQRT_MIN) { /* * Avoid a possible underflow caused by y/A. For casinh this * would be legitimate, but will be picked up by invoking atan2 * later on. For cacos this would not be legitimate. */ *B_is_usable = 0; *sqrt_A2my2 = A * (2 / DBL_EPSILON); *new_y = y * (2 / DBL_EPSILON); return; } /* B = (|z+I| - |z-I|) / 2 = y/A */ *B = y / A; *B_is_usable = 1; if (*B > B_crossover) { *B_is_usable = 0; /* * Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1). * sqrt_A2my2 = sqrt(Amy*(A+y)) */ if (y == 1 && x < DBL_EPSILON / 128) { /* * fp is of order x^2, and fm = x/2. * A = 1 (inexactly). */ *sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2); } else if (x >= DBL_EPSILON * fabs(y - 1)) { /* * Underflow will not occur because * x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN * and * x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN */ Amy = f(x, y + 1, R) + f(x, y - 1, S); *sqrt_A2my2 = sqrt(Amy * (A + y)); } else if (y > 1) { /* * fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and * A = y (inexactly). * * y < RECIP_EPSILON. So the following * scaling should avoid any underflow problems. */ *sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y / sqrt((y + 1) * (y - 1)); *new_y = y * (4 / DBL_EPSILON / DBL_EPSILON); } else { /* if (y < 1) */ /* * fm = 1-y >= DBL_EPSILON, fp is of order x^2, and * A = 1 (inexactly). */ *sqrt_A2my2 = sqrt((1 - y) * (1 + y)); } } } /* * casinh(z) = z + O(z^3) as z -> 0 * * casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity * The above formula works for the imaginary part as well, because * Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3) * as z -> infinity, uniformly in y */ double complex casinh(double complex z) { double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y; int B_is_usable; double complex w; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); if (isnan(x) || isnan(y)) { /* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */ if (isinf(x)) return (CMPLX(x, y + y)); /* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */ if (isinf(y)) return (CMPLX(y, x + x)); /* casinh(NaN + I*0) = NaN + I*0 */ if (y == 0) return (CMPLX(x + x, y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { /* clog...() will raise inexact unless x or y is infinite. */ if (signbit(x) == 0) w = clog_for_large_values(z) + m_ln2; else w = clog_for_large_values(-z) + m_ln2; return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y))); } /* Avoid spuriously raising inexact for z = 0. */ if (x == 0 && y == 0) return (z); /* All remaining cases are inexact. */ 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 = asin(B); else ry = atan2(new_y, sqrt_A2my2); return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* * casin(z) = reverse(casinh(reverse(z))) * where reverse(x + I*y) = y + I*x = I*conj(z). */ double complex casin(double complex z) { double complex w = casinh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(w), creal(w))); } /* * cacos(z) = PI/2 - casin(z) * but do the computation carefully so cacos(z) is accurate when z is * close to 1. * * cacos(z) = PI/2 - z + O(z^3) as z -> 0 * * cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity * The above formula works for the real part as well, because * Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3) * as z -> infinity, uniformly in y */ double complex cacos(double complex z) { double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x; int sx, sy; int B_is_usable; double complex w; x = creal(z); y = cimag(z); sx = signbit(x); sy = signbit(y); ax = fabs(x); ay = fabs(y); if (isnan(x) || isnan(y)) { /* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */ if (isinf(x)) return (CMPLX(y + y, -INFINITY)); /* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */ if (isinf(y)) return (CMPLX(x + x, -y)); /* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */ if (x == 0) return (CMPLX(pio2_hi + pio2_lo, y + y)); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) { /* clog...() will raise inexact unless x or y is infinite. */ w = clog_for_large_values(z); rx = fabs(cimag(w)); ry = creal(w) + m_ln2; if (sy == 0) ry = -ry; return (CMPLX(rx, ry)); } /* Avoid spuriously raising inexact for z = 1. */ if (x == 1 && y == 0) return (CMPLX(0, -y)); /* All remaining cases are inexact. */ raise_inexact(); if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4) return (CMPLX(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 = acos(B); else rx = acos(-B); } else { if (sx == 0) rx = atan2(sqrt_A2mx2, new_x); else rx = atan2(sqrt_A2mx2, -new_x); } if (sy == 0) ry = -ry; return (CMPLX(rx, ry)); } /* * cacosh(z) = I*cacos(z) or -I*cacos(z) * where the sign is chosen so Re(cacosh(z)) >= 0. */ double complex cacosh(double complex z) { double complex w; double rx, ry; w = cacos(z); rx = creal(w); ry = cimag(w); /* cacosh(NaN + I*NaN) = NaN + I*NaN */ if (isnan(rx) && isnan(ry)) return (CMPLX(ry, rx)); /* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */ /* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */ if (isnan(rx)) return (CMPLX(fabs(ry), rx)); /* cacosh(0 + I*NaN) = NaN + I*NaN */ if (isnan(ry)) return (CMPLX(ry, ry)); return (CMPLX(fabs(ry), copysign(rx, cimag(z)))); } /* * Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON. */ static double complex clog_for_large_values(double complex z) { double x, y; double ax, ay, t; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); if (ax < ay) { t = ax; ax = ay; ay = t; } /* * Avoid overflow in hypot() when x and y are both very large. - * Divide x and y by E, and then add 1 to the logarithm. This depends - * on E being larger than sqrt(2). + * Divide x and y by E, and then add 1 to the logarithm. This + * depends on E being larger than sqrt(2), since the return value of + * hypot cannot overflow if neither argument is greater in magnitude + * than 1/sqrt(2) of the maximum value of the return type. Likewise + * this determines the necessary threshold for using this method + * (however, actually use 1/2 instead as it is simpler). + * * Dividing by E causes an insignificant loss of accuracy; however * this method is still poor since it is uneccessarily slow. */ if (ax > DBL_MAX / 2) return (CMPLX(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x))); /* * Avoid overflow when x or y is large. Avoid underflow when x or * y is small. */ if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN) return (CMPLX(log(hypot(x, y)), atan2(y, x))); return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x))); } /* * ================= * | catanh, catan | * ================= */ /* * sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow). * Assumes x*x and y*y will not overflow. * Assumes x and y are finite. * Assumes y is non-negative. * Assumes fabs(x) >= DBL_EPSILON. */ static inline double sum_squares(double x, double y) { /* Avoid underflow when y is small. */ if (y < SQRT_MIN) return (x * x); return (x * x + y * y); } /* * real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y). * Assumes x and y are not NaN, and one of x and y is larger than * RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use * the code creal(1/z), because the imaginary part may produce an unwanted * underflow. * This is only called in a context where inexact is always raised before * the call, so no effort is made to avoid or force inexact. */ static inline double real_part_reciprocal(double x, double y) { double scale; uint32_t hx, hy; int32_t ix, iy; /* * This code is inspired by the C99 document n1124.pdf, Section G.5.1, * example 2. */ GET_HIGH_WORD(hx, x); ix = hx & 0x7ff00000; GET_HIGH_WORD(hy, y); iy = hy & 0x7ff00000; #define BIAS (DBL_MAX_EXP - 1) /* XXX more guard digits are useful iff there is extra precision. */ #define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */ if (ix - iy >= CUTOFF << 20 || isinf(x)) return (1 / x); /* +-Inf -> +-0 is special */ if (iy - ix >= CUTOFF << 20) return (x / y / y); /* should avoid double div, but hard */ if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20) return (x / (x * x + y * y)); scale = 1; SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */ x *= scale; y *= scale; return (x / (x * x + y * y) * scale); } /* * catanh(z) = log((1+z)/(1-z)) / 2 * = log1p(4*x / |z-1|^2) / 4 * + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2 * * catanh(z) = z + O(z^3) as z -> 0 * * catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity * The above formula works for the real part as well, because * Re(catanh(z)) = x/|z|^2 + O(x/z^4) * as z -> infinity, uniformly in x */ double complex catanh(double complex z) { double x, y, ax, ay, rx, ry; x = creal(z); y = cimag(z); ax = fabs(x); ay = fabs(y); /* This helps handle many cases. */ if (y == 0 && ax <= 1) return (CMPLX(atanh(x), y)); /* To ensure the same accuracy as atan(), and to filter out z = 0. */ if (x == 0) return (CMPLX(x, atan(y))); if (isnan(x) || isnan(y)) { /* catanh(+-Inf + I*NaN) = +-0 + I*NaN */ if (isinf(x)) return (CMPLX(copysign(0, x), y + y)); /* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */ if (isinf(y)) return (CMPLX(copysign(0, x), copysign(pio2_hi + pio2_lo, y))); /* * All other cases involving NaN return NaN + I*NaN. * C99 leaves it optional whether to raise invalid if one of * the arguments is not NaN, so we opt not to raise it. */ return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0))); } if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) return (CMPLX(real_part_reciprocal(x, y), copysign(pio2_hi + pio2_lo, y))); if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) { /* * z = 0 was filtered out above. All other cases must raise * inexact, but this is the only only that needs to do it * explicitly. */ raise_inexact(); return (z); } if (ax == 1 && ay < DBL_EPSILON) rx = (m_ln2 - log(ay)) / 2; else rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4; if (ax == 1) ry = atan2(2, -ay) / 2; else if (ay < DBL_EPSILON) ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2; else ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2; return (CMPLX(copysign(rx, x), copysign(ry, y))); } /* * catan(z) = reverse(catanh(reverse(z))) * where reverse(x + I*y) = y + I*x = I*conj(z). */ double complex catan(double complex z) { double complex w = catanh(CMPLX(cimag(z), creal(z))); return (CMPLX(cimag(w), creal(w))); } Index: stable/11/lib/msun/src/s_csqrtl.c =================================================================== --- stable/11/lib/msun/src/s_csqrtl.c (revision 323474) +++ stable/11/lib/msun/src/s_csqrtl.c (revision 323475) @@ -1,108 +1,116 @@ /*- * 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 (LDBL_MAX / 2.414213562373095048801688724209698L) +/* + * We risk spurious overflow for components >= LDBL_MAX / (1 + sqrt(2)). + * Rather than determining the fully precise value at which we might + * overflow, just use a threshold of approximately LDBL_MAX / 4. + */ +#if LDBL_MAX_EXP != 0x4000 +#error "Unsupported long double format" +#else +#define THRESH 0x1p16382L +#endif 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); } Index: stable/11/lib/msun/tests/csqrt_test.c =================================================================== --- stable/11/lib/msun/tests/csqrt_test.c (revision 323474) +++ stable/11/lib/msun/tests/csqrt_test.c (revision 323475) @@ -1,295 +1,376 @@ /*- * Copyright (c) 2007 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. */ /* * Tests for csqrt{,f}() */ #include __FBSDID("$FreeBSD$"); #include #include #include #include #include #include #include "test-utils.h" /* * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf(). * The latter two convert to float or double, respectively, and test csqrtf() * and csqrt() with the same arguments. */ static long double complex (*t_csqrt)(long double complex); static long double complex _csqrtf(long double complex d) { return (csqrtf((float complex)d)); } static long double complex _csqrt(long double complex d) { return (csqrt((double complex)d)); } #pragma STDC CX_LIMITED_RANGE OFF /* * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. * Fail an assertion if they differ. */ static void assert_equal(long double complex d1, long double complex d2) { assert(cfpequal(d1, d2)); } /* * Test csqrt for some finite arguments where the answer is exact. * (We do not test if it produces correctly rounded answers when the * result is inexact, nor do we check whether it throws spurious * exceptions.) */ static void test_finite(void) { static const double tests[] = { /* csqrt(a + bI) = x + yI */ /* a b x y */ 0, 8, 2, 2, 0, -8, 2, -2, 4, 0, 2, 0, -4, 0, 0, 2, 3, 4, 2, 1, 3, -4, 2, -1, -3, 4, 1, 2, -3, -4, 1, -2, 5, 12, 3, 2, 7, 24, 4, 3, 9, 40, 5, 4, 11, 60, 6, 5, 13, 84, 7, 6, 33, 56, 7, 4, 39, 80, 8, 5, 65, 72, 9, 4, 987, 9916, 74, 67, 5289, 6640, 83, 40, 460766389075.0, 16762287900.0, 678910, 12345 }; /* * We also test some multiples of the above arguments. This * array defines which multiples we use. Note that these have * to be small enough to not cause overflow for float precision * with all of the constants in the above table. */ static const double mults[] = { 1, 2, 3, 13, 16, 0x1.p30, 0x1.p-30, }; double a, b; double x, y; unsigned i, j; for (i = 0; i < nitems(tests); i += 4) { for (j = 0; j < nitems(mults); j++) { a = tests[i] * mults[j] * mults[j]; b = tests[i + 1] * mults[j] * mults[j]; x = tests[i + 2] * mults[j]; y = tests[i + 3] * mults[j]; assert(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y)); } } } /* * Test the handling of +/- 0. */ static void test_zeros(void) { assert_equal(t_csqrt(CMPLXL(0.0, 0.0)), CMPLXL(0.0, 0.0)); assert_equal(t_csqrt(CMPLXL(-0.0, 0.0)), CMPLXL(0.0, 0.0)); assert_equal(t_csqrt(CMPLXL(0.0, -0.0)), CMPLXL(0.0, -0.0)); assert_equal(t_csqrt(CMPLXL(-0.0, -0.0)), CMPLXL(0.0, -0.0)); } /* * Test the handling of infinities when the other argument is not NaN. */ static void test_infinities(void) { static const double vals[] = { 0.0, -0.0, 42.0, -42.0, INFINITY, -INFINITY, }; unsigned i; for (i = 0; i < nitems(vals); i++) { if (isfinite(vals[i])) { assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])), CMPLXL(0.0, copysignl(INFINITY, vals[i]))); assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])), CMPLXL(INFINITY, copysignl(0.0, vals[i]))); } assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)), CMPLXL(INFINITY, INFINITY)); assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)), CMPLXL(INFINITY, -INFINITY)); } } /* * Test the handling of NaNs. */ static void test_nans(void) { assert(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY); assert(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN))))); assert(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN))))); assert(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN))))); assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)), CMPLXL(INFINITY, INFINITY)); assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)), CMPLXL(INFINITY, -INFINITY)); assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN)); assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN)); } /* * Test whether csqrt(a + bi) works for inputs that are large enough to - * cause overflow in hypot(a, b) + a. In this case we are using - * csqrt(115 + 252*I) == 14 + 9*I - * scaled up to near MAX_EXP. + * cause overflow in hypot(a, b) + a. Each of the tests is scaled up to + * near MAX_EXP. */ static void test_overflow(int maxexp) { long double a, b; long double complex result; + int exp, i; - a = ldexpl(115 * 0x1p-8, maxexp); - b = ldexpl(252 * 0x1p-8, maxexp); - result = t_csqrt(CMPLXL(a, b)); - assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2)); - assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2)); + assert(maxexp > 0 && maxexp % 2 == 0); + + for (i = 0; i < 4; i++) { + exp = maxexp - 2 * i; + + /* csqrt(115 + 252*I) == 14 + 9*I */ + a = ldexpl(115 * 0x1p-8, exp); + b = ldexpl(252 * 0x1p-8, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(14 * 0x1p-4, exp / 2)); + assert(cimagl(result) == ldexpl(9 * 0x1p-4, exp / 2)); + + /* csqrt(-11 + 60*I) = 5 + 6*I */ + a = ldexpl(-11 * 0x1p-6, exp); + b = ldexpl(60 * 0x1p-6, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(5 * 0x1p-3, exp / 2)); + assert(cimagl(result) == ldexpl(6 * 0x1p-3, exp / 2)); + + /* csqrt(225 + 0*I) == 15 + 0*I */ + a = ldexpl(225 * 0x1p-8, exp); + b = 0; + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(15 * 0x1p-4, exp / 2)); + assert(cimagl(result) == 0); + } } +/* + * Test that precision is maintained for some large squares. Set all or + * some bits in the lower mantdig/2 bits, square the number, and try to + * recover the sqrt. Note: + * (x + xI)**2 = 2xxI + */ +static void +test_precision(int maxexp, int mantdig) +{ + long double b, x; + long double complex result; + uint64_t mantbits, sq_mantbits; + int exp, i; + + assert(maxexp > 0 && maxexp % 2 == 0); + assert(mantdig <= 64); + mantdig = rounddown(mantdig, 2); + + for (exp = 0; exp <= maxexp; exp += 2) { + mantbits = ((uint64_t)1 << (mantdig / 2 )) - 1; + for (i = 0; + i < 100 && mantbits > ((uint64_t)1 << (mantdig / 2 - 1)); + i++, mantbits--) { + sq_mantbits = mantbits * mantbits; + /* + * sq_mantibts is a mantdig-bit number. Divide by + * 2**mantdig to normalize it to [0.5, 1), where, + * note, the binary power will be -1. Raise it by + * 2**exp for the test. exp is even. Lower it by + * one to reach a final binary power which is also + * even. The result should be exactly + * representable, given that mantdig is less than or + * equal to the available precision. + */ + b = ldexpl((long double)sq_mantbits, + exp - 1 - mantdig); + x = ldexpl(mantbits, (exp - 2 - mantdig) / 2); + assert(b == x * x * 2); + result = t_csqrt(CMPLXL(0, b)); + assert(creall(result) == x); + assert(cimagl(result) == x); + } + } +} + int main(void) { - printf("1..15\n"); + printf("1..18\n"); /* Test csqrt() */ t_csqrt = _csqrt; test_finite(); printf("ok 1 - csqrt\n"); test_zeros(); printf("ok 2 - csqrt\n"); test_infinities(); printf("ok 3 - csqrt\n"); test_nans(); printf("ok 4 - csqrt\n"); test_overflow(DBL_MAX_EXP); printf("ok 5 - csqrt\n"); + test_precision(DBL_MAX_EXP, DBL_MANT_DIG); + printf("ok 6 - csqrt\n"); + /* Now test csqrtf() */ t_csqrt = _csqrtf; test_finite(); - printf("ok 6 - csqrt\n"); + printf("ok 7 - csqrt\n"); test_zeros(); - printf("ok 7 - csqrt\n"); + printf("ok 8 - csqrt\n"); test_infinities(); - printf("ok 8 - csqrt\n"); + printf("ok 9 - csqrt\n"); test_nans(); - printf("ok 9 - csqrt\n"); + printf("ok 10 - csqrt\n"); test_overflow(FLT_MAX_EXP); - printf("ok 10 - csqrt\n"); + printf("ok 11 - csqrt\n"); + test_precision(FLT_MAX_EXP, FLT_MANT_DIG); + printf("ok 12 - csqrt\n"); + /* Now test csqrtl() */ t_csqrt = csqrtl; test_finite(); - printf("ok 11 - csqrt\n"); + printf("ok 13 - csqrt\n"); test_zeros(); - printf("ok 12 - csqrt\n"); + printf("ok 14 - csqrt\n"); test_infinities(); - printf("ok 13 - csqrt\n"); + printf("ok 15 - csqrt\n"); test_nans(); - printf("ok 14 - csqrt\n"); + printf("ok 16 - csqrt\n"); test_overflow(LDBL_MAX_EXP); - printf("ok 15 - csqrt\n"); + printf("ok 17 - csqrt\n"); + + test_precision(LDBL_MAX_EXP, +#ifndef __i386__ + LDBL_MANT_DIG +#else + DBL_MANT_DIG +#endif + ); + printf("ok 18 - csqrt\n"); return (0); } Index: stable/11 =================================================================== --- stable/11 (revision 323474) +++ stable/11 (revision 323475) Property changes on: stable/11 ___________________________________________________________________ Modified: svn:mergeinfo ## -0,0 +0,1 ## Merged /head:r323003-323004