From cf551d94ef2f28f9710c71918393a3a0ab3572b2 Mon Sep 17 00:00:00 2001
From: Ryan Libby <rlibby@FreeBSD.org>
Date: Tue, 29 Aug 2017 22:32:29 +0000
Subject: [PATCH] lib/msun: avoid referring to broken LDBL_MAX

LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html

Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'".  Avoid this by
referring to proxy values instead.

Reviewed by:	bde
Approved by:	markj (mentor)
Sponsored by:	Dell EMC Isilon
---
 lib/msun/src/catrig.c   |  9 +++++++--
 lib/msun/src/catrigl.c  |  7 ++++++-
 lib/msun/src/s_csqrtl.c | 12 ++++++++++--
 3 files changed, 23 insertions(+), 5 deletions(-)

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)