svn commit: r336488 - head/lib/msun/src

[Bruce Evans]

Author: bde
Date: Thu Jul 19 15:04:10 2018
New Revision: 336488
URL: https://svnweb.freebsd.org/changeset/base/336488

Log:
  Fix spurious and extra underflows and resulting inaccuracies for some cases
  with 1 huge component and 1 tiny (but nowhere near denormal) component.
  Rescale earlier so that a scale factor of 2 can be combined with a non-
  scale divisor of 2, so that the division doesn't shift out a bit.  In the
  usual case where the scale factor is just 1, the division may shift out a
  bit, but then the underflow is not spurious and the inaccuracies are harder
  to fix.

Modified:
  head/lib/msun/src/s_csqrt.c
  head/lib/msun/src/s_csqrtl.c

Modified: head/lib/msun/src/s_csqrt.c
==============================================================================
--- head/lib/msun/src/s_csqrt.c	Thu Jul 19 14:41:46 2018	(r336487)
+++ head/lib/msun/src/s_csqrt.c	Thu Jul 19 15:04:10 2018	(r336488)
@@ -99,15 +99,15 @@ csqrt(double complex z)
 	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrt((a + hypot(a, b)) * 0.5);
-		rx = t;
-		ry = b / (2 * t);
+		rx = scale * t;
+		ry = scale * b / (2 * t);
 	} else {
 		t = sqrt((-a + hypot(a, b)) * 0.5);
-		rx = fabs(b) / (2 * t);
-		ry = copysign(t, b);
+		rx = scale * fabs(b) / (2 * t);
+		ry = copysign(scale * t, b);
 	}
 
-	return (CMPLX(rx * scale, ry * scale));
+	return (CMPLX(rx, ry));
 }
 
 #if LDBL_MANT_DIG == 53

Modified: head/lib/msun/src/s_csqrtl.c
==============================================================================
--- head/lib/msun/src/s_csqrtl.c	Thu Jul 19 14:41:46 2018	(r336487)
+++ head/lib/msun/src/s_csqrtl.c	Thu Jul 19 15:04:10 2018	(r336488)
@@ -114,13 +114,13 @@ csqrtl(long double complex z)
 	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrtl((a + hypotl(a, b)) * 0.5);
-		rx = t;
-		ry = b / (2 * t);
+		rx = scale * t;
+		ry = scale * b / (2 * t);
 	} else {
 		t = sqrtl((-a + hypotl(a, b)) * 0.5);
-		rx = fabsl(b) / (2 * t);
-		ry = copysignl(t, b);
+		rx = scale * fabsl(b) / (2 * t);
+		ry = copysignl(scale * t, b);
 	}
 
-	return (CMPLXL(rx * scale, ry * scale));
+	return (CMPLXL(rx, ry));
 }