MFC r323003,r323004:

  r323003:
    lib/msun: avoid referring to broken LDBL_MAX
  r323004:
    lib/msun: add more csqrt unit tests for precision and overflow

3 files changed, 117 insertions, 23 deletions

diff --git a/lib/msun/src/catrig.c b/lib/msun/src/catrig.c
index f392862..a57cc9e 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/s_csqrtl.c b/lib/msun/src/s_csqrtl.c
index ea2ef27..01cc409 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)
diff --git a/lib/msun/tests/csqrt_test.c b/lib/msun/tests/csqrt_test.c
index 89de14b..9f596ed 100644
--- a/lib/msun/tests/csqrt_test.c
+++ b/lib/msun/tests/csqrt_test.c
@@ -214,28 +214,94 @@ test_nans(void)
 
 /*
  * Test whether csqrt(a + bi) works for inputs that are large enough to
- * cause overflow in hypot(a, b) + a. In this case we are using
- *	csqrt(115 + 252*I) == 14 + 9*I
- * scaled up to near MAX_EXP.
+ * cause overflow in hypot(a, b) + a.  Each of the tests is scaled up to
+ * near MAX_EXP.
  */
 static void
 test_overflow(int maxexp)
 {
 	long double a, b;
 	long double complex result;
+	int exp, i;
+
+	assert(maxexp > 0 && maxexp % 2 == 0);
+
+	for (i = 0; i < 4; i++) {
+		exp = maxexp - 2 * i;
+
+		/* csqrt(115 + 252*I) == 14 + 9*I */
+		a = ldexpl(115 * 0x1p-8, exp);
+		b = ldexpl(252 * 0x1p-8, exp);
+		result = t_csqrt(CMPLXL(a, b));
+		assert(creall(result) == ldexpl(14 * 0x1p-4, exp / 2));
+		assert(cimagl(result) == ldexpl(9 * 0x1p-4, exp / 2));
+
+		/* csqrt(-11 + 60*I) = 5 + 6*I */
+		a = ldexpl(-11 * 0x1p-6, exp);
+		b = ldexpl(60 * 0x1p-6, exp);
+		result = t_csqrt(CMPLXL(a, b));
+		assert(creall(result) == ldexpl(5 * 0x1p-3, exp / 2));
+		assert(cimagl(result) == ldexpl(6 * 0x1p-3, exp / 2));
+
+		/* csqrt(225 + 0*I) == 15 + 0*I */
+		a = ldexpl(225 * 0x1p-8, exp);
+		b = 0;
+		result = t_csqrt(CMPLXL(a, b));
+		assert(creall(result) == ldexpl(15 * 0x1p-4, exp / 2));
+		assert(cimagl(result) == 0);
+	}
+}
 
-	a = ldexpl(115 * 0x1p-8, maxexp);
-	b = ldexpl(252 * 0x1p-8, maxexp);
-	result = t_csqrt(CMPLXL(a, b));
-	assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2));
-	assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2));
+/*
+ * Test that precision is maintained for some large squares.  Set all or
+ * some bits in the lower mantdig/2 bits, square the number, and try to
+ * recover the sqrt.  Note:
+ * 	(x + xI)**2 = 2xxI
+ */
+static void
+test_precision(int maxexp, int mantdig)
+{
+	long double b, x;
+	long double complex result;
+	uint64_t mantbits, sq_mantbits;
+	int exp, i;
+
+	assert(maxexp > 0 && maxexp % 2 == 0);
+	assert(mantdig <= 64);
+	mantdig = rounddown(mantdig, 2);
+
+	for (exp = 0; exp <= maxexp; exp += 2) {
+		mantbits = ((uint64_t)1 << (mantdig / 2 )) - 1;
+		for (i = 0;
+		     i < 100 && mantbits > ((uint64_t)1 << (mantdig / 2 - 1));
+		     i++, mantbits--) {
+			sq_mantbits = mantbits * mantbits;
+			/*
+			 * sq_mantibts is a mantdig-bit number.  Divide by
+			 * 2**mantdig to normalize it to [0.5, 1), where,
+			 * note, the binary power will be -1.  Raise it by
+			 * 2**exp for the test.  exp is even.  Lower it by
+			 * one to reach a final binary power which is also
+			 * even.  The result should be exactly
+			 * representable, given that mantdig is less than or
+			 * equal to the available precision.
+			 */
+			b = ldexpl((long double)sq_mantbits,
+			    exp - 1 - mantdig);
+			x = ldexpl(mantbits, (exp - 2 - mantdig) / 2);
+			assert(b == x * x * 2);
+			result = t_csqrt(CMPLXL(0, b));
+			assert(creall(result) == x);
+			assert(cimagl(result) == x);
+		}
+	}
 }
 
 int
 main(void)
 {
 
-	printf("1..15\n");
+	printf("1..18\n");
 
 	/* Test csqrt() */
 	t_csqrt = _csqrt;
@@ -255,41 +321,56 @@ main(void)
 	test_overflow(DBL_MAX_EXP);
 	printf("ok 5 - csqrt\n");
 
+	test_precision(DBL_MAX_EXP, DBL_MANT_DIG);
+	printf("ok 6 - csqrt\n");
+
 	/* Now test csqrtf() */
 	t_csqrt = _csqrtf;
 
 	test_finite();
-	printf("ok 6 - csqrt\n");
+	printf("ok 7 - csqrt\n");
 
 	test_zeros();
-	printf("ok 7 - csqrt\n");
+	printf("ok 8 - csqrt\n");
 
 	test_infinities();
-	printf("ok 8 - csqrt\n");
+	printf("ok 9 - csqrt\n");
 
 	test_nans();
-	printf("ok 9 - csqrt\n");
+	printf("ok 10 - csqrt\n");
 
 	test_overflow(FLT_MAX_EXP);
-	printf("ok 10 - csqrt\n");
+	printf("ok 11 - csqrt\n");
+
+	test_precision(FLT_MAX_EXP, FLT_MANT_DIG);
+	printf("ok 12 - csqrt\n");
 
 	/* Now test csqrtl() */
 	t_csqrt = csqrtl;
 
 	test_finite();
-	printf("ok 11 - csqrt\n");
+	printf("ok 13 - csqrt\n");
 
 	test_zeros();
-	printf("ok 12 - csqrt\n");
+	printf("ok 14 - csqrt\n");
 
 	test_infinities();
-	printf("ok 13 - csqrt\n");
+	printf("ok 15 - csqrt\n");
 
 	test_nans();
-	printf("ok 14 - csqrt\n");
+	printf("ok 16 - csqrt\n");
 
 	test_overflow(LDBL_MAX_EXP);
-	printf("ok 15 - csqrt\n");
+	printf("ok 17 - csqrt\n");
+
+	test_precision(LDBL_MAX_EXP,
+#ifndef __i386__
+	    LDBL_MANT_DIG
+#else
+	    DBL_MANT_DIG
+#endif
+	    );
+	printf("ok 18 - csqrt\n");
 
 	return (0);
 }