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));
}
_______________________________________________
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"
