d/freebsd-nq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix spurious and extra underflows and resulting inaccuracies for some cases

Browse Source 
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.
This commit is contained in:
Bruce Evans 2018-07-19 15:04:10 +00:00
parent 50c8bd4e53
commit b7092eef4d
2 changed files with 10 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
lib/msun/src/s_csqrt.c
View File

@ -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

10
lib/msun/src/s_csqrtl.c
View File

@ -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));
}