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