diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c index ebfe3567c168..b304f213fa6d 100644 --- a/lib/msun/src/catrig.c +++ b/lib/msun/src/catrig.c @@ -469,8 +469,13 @@ clog_for_large_values(double complex z) /* * 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. */ diff --git a/lib/msun/src/catrigl.c b/lib/msun/src/catrigl.c index 41357e864456..b79857082700 100644 --- a/lib/msun/src/catrigl.c +++ b/lib/msun/src/catrigl.c @@ -57,10 +57,15 @@ __FBSDID("$FreeBSD$"); #undef signbit #define signbit(x) (__builtin_signbitl(x)) +#if LDBL_MAX_EXP != 0x4000 +#error "Unsupported long double format" +#endif + static const long double A_crossover = 10, B_crossover = 0.6417, FOUR_SQRT_MIN = 0x1p-8189L, +HALF_MAX = 0x1p16383L, QUARTER_SQRT_MAX = 0x1p8189L, RECIP_EPSILON = 1 / LDBL_EPSILON, SQRT_MIN = 0x1p-8191L; @@ -307,7 +312,7 @@ clog_for_large_values(long double complex z) ay = t; } - if (ax > LDBL_MAX / 2) + if (ax > HALF_MAX) return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, atan2l(y, x))); diff --git a/lib/msun/src/s_csqrtl.c b/lib/msun/src/s_csqrtl.c index ea2ef27b5088..01cc4090e964 100644 --- a/lib/msun/src/s_csqrtl.c +++ b/lib/msun/src/s_csqrtl.c @@ -42,8 +42,16 @@ __FBSDID("$FreeBSD$"); */ #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)