diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c index cde8792c36f6..ca4924fc2bdc 100644 --- a/lib/msun/src/s_csqrt.c +++ b/lib/msun/src/s_csqrt.c @@ -35,16 +35,7 @@ __FBSDID("$FreeBSD$"); #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 >= DBL_MAX / (1 + sqrt(2)). */ +/* For avoiding overflow for components >= DBL_MAX / (1 + sqrt(2)). */ #define THRESH 0x1.a827999fcef32p+1022 double complex diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c index 0f36aa0ffb51..6836a3903c25 100644 --- a/lib/msun/src/s_csqrtf.c +++ b/lib/msun/src/s_csqrtf.c @@ -34,20 +34,14 @@ __FBSDID("$FreeBSD$"); #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 - float complex csqrtf(float complex z) { - float a = crealf(z), b = cimagf(z); double t; + float a, b; + + a = creal(z); + b = cimag(z); /* Handle special cases. */ if (z == 0) @@ -82,9 +76,9 @@ csqrtf(float complex z) */ if (a >= 0) { t = sqrt((a + hypot(a, b)) * 0.5); - return (CMPLXF(t, b / (2.0 * t))); + return (CMPLXF(t, b / (2 * t))); } else { t = sqrt((-a + hypot(a, b)) * 0.5); - return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b))); + return (CMPLXF(fabsf(b) / (2 * t), copysignf(t, b))); } } diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c index 8aefbeab1a0b..5ab9406558af 100644 --- a/lib/msun/src/s_csqrtl.c +++ b/lib/msun/src/s_csqrtl.c @@ -36,25 +36,18 @@ __FBSDID("$FreeBSD$"); #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)). - * Rather than determining the fully precise value at which we might - * overflow, just use a threshold of approximately LDBL_MAX / 4. + * THRESH is now calculated portably (up to 113-bit precision). However, + * the denormal threshold is hard-coded for a 15-bit exponent with the usual + * bias. s_logl.c and e_hypotl have less hard-coding but end up requiring + * the same for the exponent and more for the mantissa. */ #if LDBL_MAX_EXP != 0x4000 #error "Unsupported long double format" -#else -#define THRESH 0x1p16382L #endif +/* For avoiding overflow for components >= LDBL_MAX / (1 + sqrt(2)). */ +#define THRESH (LDBL_MAX / 2.414213562373095048801688724209698L) + long double complex csqrtl(long double complex z) {