Author: rlibby
Date: Tue Aug 29 22:37:24 2017
New Revision: 323004
URL: https://svnweb.freebsd.org/changeset/base/323004
Log:
lib/msun: add more csqrt unit tests for precision and overflow
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
Modified:
head/lib/msun/tests/csqrt_test.c
Modified: head/lib/msun/tests/csqrt_test.c
==============================================================================
--- head/lib/msun/tests/csqrt_test.c Tue Aug 29 22:32:29 2017
(r323003)
+++ head/lib/msun/tests/csqrt_test.c Tue Aug 29 22:37:24 2017
(r323004)
@@ -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;
- 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));
+ 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);
+ }
}
+/*
+ * 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);
}
_______________________________________________
svn-src-head@freebsd.org mailing list
https://lists.freebsd.org/mailman/listinfo/svn-src-head
To unsubscribe, send any mail to "svn-src-head-unsubscr...@freebsd.org"
