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_csqrtl.c b/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
index a975833..e58e4ae 100644
--- a/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
+++ b/libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
@@ -27,7 +27,7 @@
*/
#include <sys/cdefs.h>
-__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 326219 2017-11-26 02:00:33Z pfg $");
+__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 336488 2018-07-19 15:04:10Z bde $");
#include <complex.h>
#include <float.h>
@@ -36,32 +36,28 @@
#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.
+ * Several thresholds require a 15-bit exponent and also 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
+/*
+ * Overflow must be avoided for components >= LDBL_MAX / (1 + sqrt(2)).
+ * The precise threshold is nontrivial to determine and spell, so use a
+ * lower threshold of approximaely LDBL_MAX / 4, and don't use LDBL_MAX
+ * to spell this since LDBL_MAX is broken on i386 (it overflows in 53-bit
+ * precision).
+ */
+#define THRESH 0x1p16382L
+
long double complex
csqrtl(long double complex z)
{
long double complex result;
- long double a, b;
- long double t;
- int scale;
+ long double a, b, rx, ry, scale, t;
a = creall(z);
b = cimagl(z);
@@ -73,7 +69,7 @@
return (CMPLXL(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return (CMPLXL(a, t)); /* return NaN + NaN i */
+ return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -87,32 +83,44 @@
else
return (CMPLXL(a, copysignl(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 (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */
+ }
/* Scale to avoid overflow. */
if (fabsl(a) >= THRESH || fabsl(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 (fabsl(a) >= 0x1p-16380L)
+ a *= 0.25;
+ if (fabsl(b) >= 0x1p-16380L)
+ b *= 0.25;
+ scale = 2;
} else {
- scale = 0;
+ scale = 1;
+ }
+
+ /* Scale to reduce inaccuracies when both components are denormal. */
+ if (fabsl(a) < 0x1p-16382L && fabsl(b) < 0x1p-16382L) {
+ a *= 0x1p64;
+ b *= 0x1p64;
+ scale = 0x1p-32;
}
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrtl((a + hypotl(a, b)) * 0.5);
- result = CMPLXL(t, b / (2 * t));
+ rx = scale * t;
+ ry = scale * b / (2 * t);
} else {
t = sqrtl((-a + hypotl(a, b)) * 0.5);
- result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b));
+ rx = scale * fabsl(b) / (2 * t);
+ ry = copysignl(scale * t, b);
}
- /* Rescale. */
- if (scale)
- return (result * 2);
- else
- return (result);
+ return (CMPLXL(rx, ry));
}