Update to FreeBSD libm r336665.

This reverts commit 253a8306316cedfd6fd3e3a169fbffe4cac04035 and moves
us forward to a revision that contains fixes for the problem with the
previous attempt.

This also makes sincos(3)/sincosf(3)/sincosl(3) available to `_BSD_SOURCE`
as well as `_GNU_SOURCE`.

The new FreeBSD libm code requires the FreeBSD `__CONCAT` macro, and all
our existing callers are FreeBSD too, so update that.

There's also an assumption that <complex.h> drags in <math.h> which isn't
true for us, so work around that with `-include` in the makefile. This
then causes clang to recognize a bug -- returning from a void function --
in our fake (LP32) sincosl(3), so fix that too.

Bug: http://b/111710419
Change-Id: I84703ad844f8afde6ec6b11604ab3c096ccb62c3
Test: ran tests
diff --git a/libm/upstream-freebsd/lib/msun/src/s_csqrt.c b/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
index ad8ef4c..215b7d7 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csqrt.c
@@ -27,7 +27,7 @@
  */
 
 #include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrt.c 326219 2017-11-26 02:00:33Z pfg $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrt.c 336488 2018-07-19 15:04:10Z bde $");
 
 #include <complex.h>
 #include <float.h>
@@ -35,25 +35,14 @@
 
 #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
 csqrt(double complex z)
 {
 	double complex result;
-	double a, b;
-	double t;
-	int scale;
+	double a, b, rx, ry, scale, t;
 
 	a = creal(z);
 	b = cimag(z);
@@ -65,7 +54,7 @@
 		return (CMPLX(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (CMPLX(a, t));	/* return NaN + NaN i */
+		return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
@@ -79,34 +68,46 @@
 		else
 			return (CMPLX(a, copysign(b - b, b)));
 	}
-	/*
-	 * The remaining special case (b is NaN) is handled just fine by
-	 * the normal code path below.
-	 */
+	if (isnan(b)) {
+		t = (a - a) / (a - a);	/* raise invalid */
+		return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
+	}
 
 	/* Scale to avoid overflow. */
 	if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
-		a *= 0.25;
-		b *= 0.25;
-		scale = 1;
+		/*
+		 * Don't scale a or b if this might give (spurious)
+		 * underflow.  Then the unscaled value is an equivalent
+		 * infinitesmal (or 0).
+		 */
+		if (fabs(a) >= 0x1p-1020)
+			a *= 0.25;
+		if (fabs(b) >= 0x1p-1020)
+			b *= 0.25;
+		scale = 2;
 	} else {
-		scale = 0;
+		scale = 1;
+	}
+
+	/* Scale to reduce inaccuracies when both components are denormal. */
+	if (fabs(a) < 0x1p-1022 && fabs(b) < 0x1p-1022) {
+		a *= 0x1p54;
+		b *= 0x1p54;
+		scale = 0x1p-27;
 	}
 
 	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrt((a + hypot(a, b)) * 0.5);
-		result = CMPLX(t, b / (2 * t));
+		rx = scale * t;
+		ry = scale * b / (2 * t);
 	} else {
 		t = sqrt((-a + hypot(a, b)) * 0.5);
-		result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
+		rx = scale * fabs(b) / (2 * t);
+		ry = copysign(scale * t, b);
 	}
 
-	/* Rescale. */
-	if (scale)
-		return (result * 2);
-	else
-		return (result);
+	return (CMPLX(rx, ry));
 }
 
 #if LDBL_MANT_DIG == 53