Skip site navigation (1)Skip section navigation (2) 
Date:      Thu, 19 Jul 2018 15:04:10 +0000 (UTC)
From:      Bruce Evans <bde@FreeBSD.org>
To:        src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org
Subject:   svn commit: r336488 - head/lib/msun/src
Message-ID:  <201807191504.w6JF4AFD080232@repo.freebsd.org>
next in thread | raw e-mail | index | archive | help 
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));
 }

Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201807191504.w6JF4AFD080232>