This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix hypot handling of subnormals (bug 16316, bug 16330)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Mon, 16 Dec 2013 22:16:05 +0000
- Subject: Fix hypot handling of subnormals (bug 16316, bug 16330)
- Authentication-results: sourceware.org; auth=none
This patch fixes bugs 16316 and 16330, problems with hypot with
subnormal arguments.
The general problem was that the hypot implementations for dbl-64,
ldbl-96 and ldbl-128 compute exponents / high parts of their
arguments, and sort their arguments so a >= b >= 0 (roughly; only the
exponents or high parts may get compared), before doing the main
computation that uses the exponents / high parts and the ordering to
get more accurate results - but the implementations also need to
adjust very large or small arguments to avoid internal overflows and
underflows, and the adjustments in the case of subnormals failed to
adjust the previously computed exponents / high parts or to redo the
ordering (which, when just comparing exponents, won't actually have
done anything useful for subnormals before the adjustment). This
could result in spurious underflows - or in invalid numbers being used
in subsequent computations in the ldbl-96 case - because of the high
parts / exponents of subnormals being used in code that only expected
normal values. (I don't think it would have resulted in significant
errors in results, although they'd have been slightly less accurate in
some cases.)
This patch ensures the relevant values are recomputed, and sorted by
the same logic used in sorting unadjusted values. (The logic that
checked for the second value being zero doesn't need to check for the
case where, before sorting, the first value was zero but the second
wasn't, as having just one zero, in the second place after sorting,
doesn't cause problems for the subsequent code.)
Tested x86_64; no ulps updates needed. (32-bit x86 builds are
currently broken.) Also spot-checked for mips64-linux-gnu to test the
ldbl-128 changes.
(auto-libm-test-out diffs omitted.)
2013-12-16 Joseph Myers <joseph@codesourcery.com>
[BZ #16316]
[BZ #16330]
* sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Update
values of ha and hb and sort them after adjusting subnormal
arguments.
* sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl):
Likewise.
* sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Update
values of ea and eb and sort them after adjusting subnormal
arguments.
* math/auto-libm-test-in: Do not expect some hypot tests of
subnormals to fail. Add more hypot tests.
* math/auto-libm-test-out: Regenerated.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index f468676..980b6b9 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -282,6 +282,10 @@ expm1 100000.0
expm1 max
expm1 -max
+hypot 0 0
+hypot 0 -0
+hypot -0 0
+hypot -0 -0
# hypot (x,y) == hypot (+-x, +-y).
hypot 0.7 12.4
hypot -0.7 12.4
@@ -307,9 +311,15 @@ hypot 0x1p+0 0x0.3ep-1022 no-test-inline:dbl-64
hypot 0x3p16381 0x4p16381 no-test-inline
hypot 0x1p-149 0x1p-149
hypot 0x1p-1074 0x1p-1074
-# Bug 16316: hypot broken for some subnormal arguments.
-hypot 0x1p-16445 0x1p-16445 no-test-inline xfail:ldbl-96-intel xfail:ldbl-96-m68k
-hypot 0x1p-16494 0x1p-16494 no-test-inline xfail:ldbl-96-intel xfail:ldbl-96-m68k
+hypot 0x1p-16445 0x1p-16445 no-test-inline
+hypot 0x1p-16494 0x1p-16494 no-test-inline
+hypot 0x0.fffffep-126 0x0.fp-127
+hypot 0x0.fffffep-126 0x0.fp-130
+hypot 0x0.fffffffffffffp-1022 0x0.fp-1023
+hypot 0x0.fffffffffffffp-1022 0x0.fp-1026
+hypot 0x0.ffffffp-16382 0x0.fp-16383
+hypot 0x0.ffffffp-16382 0x0.fp-16386
+hypot 0 min_subnorm
j0 -1.0
j0 0.0
diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c
index 9f4321e..500658d 100644
--- a/sysdeps/ieee754/dbl-64/e_hypot.c
+++ b/sysdeps/ieee754/dbl-64/e_hypot.c
@@ -102,6 +102,17 @@ __ieee754_hypot (double x, double y)
b *= t1;
a *= t1;
k -= 1022;
+ GET_HIGH_WORD (ha, a);
+ GET_HIGH_WORD (hb, b);
+ if (hb > ha)
+ {
+ t1 = a;
+ a = b;
+ b = t1;
+ j = ha;
+ ha = hb;
+ hb = j;
+ }
}
else /* scale a and b by 2^600 */
{
diff --git a/sysdeps/ieee754/ldbl-128/e_hypotl.c b/sysdeps/ieee754/ldbl-128/e_hypotl.c
index f5ab901..01444cf 100644
--- a/sysdeps/ieee754/ldbl-128/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-128/e_hypotl.c
@@ -89,6 +89,17 @@ __ieee754_hypotl(long double x, long double y)
b *= t1;
a *= t1;
k -= 16382;
+ GET_LDOUBLE_MSW64 (ha, a);
+ GET_LDOUBLE_MSW64 (hb, b);
+ if (hb > ha)
+ {
+ t1 = a;
+ a = b;
+ b = t1;
+ j = ha;
+ ha = hb;
+ hb = j;
+ }
} else { /* scale a and b by 2^9600 */
ha += 0x2580000000000000LL; /* a *= 2^9600 */
hb += 0x2580000000000000LL; /* b *= 2^9600 */
diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c
index 7895488..d3152f9 100644
--- a/sysdeps/ieee754/ldbl-96/e_hypotl.c
+++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c
@@ -89,6 +89,17 @@ long double __ieee754_hypotl(long double x, long double y)
b *= t1;
a *= t1;
k -= 16382;
+ GET_LDOUBLE_EXP (ea, a);
+ GET_LDOUBLE_EXP (eb, b);
+ if (eb > ea)
+ {
+ t1 = a;
+ a = b;
+ b = t1;
+ j = ea;
+ ea = eb;
+ eb = j;
+ }
} else { /* scale a and b by 2^9600 */
ea += 0x2580; /* a *= 2^9600 */
eb += 0x2580; /* b *= 2^9600 */
--
Joseph S. Myers
joseph@codesourcery.com