This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] problems in hypotf math function
- From: Bob Wilson <bwilson at tensilica dot com>
- To: libc-alpha at sources dot redhat dot com
- Cc: richard at codesourcery dot com
- Date: Mon, 01 Aug 2005 14:08:11 -0700
- Subject: [PATCH] problems in hypotf math function
The hypotf function in sysdeps/ieee754/flt-32/e_hypotf.c has problems with
scaling small input values.
Values that are not denormalized are scaled by directly incrementing the
exponent. The increment value in this case is not a floating-point value and
must not include the exponent bias. The current code for hypotf has this wrong
-- it includes the exponent bias in the increment. Richard Sandiford fixed this
in newlib (http://sources.redhat.com/ml/newlib-cvs/2002-q1/msg00048.html) when
he changed it to scale by 2^68 instead of 2^60. I don't know which scaling
factor (60 vs. 68) is better, but the exponent bias needs to be removed from the
increment in either case.
For denormal values, the scaling requires a floating-point multiplication to
normalize the mantissa. The 2^126 scale factor in this case is a floating-point
value and must include the exponent bias. It currently does not. I just
submitted a patch for this in newlib. The test program below will demonstrate
the problem when run on a platform using this version of hypotf (e.g., RedHat
Enterprise 3 on x86_64) -- without the patch, the value returned by hypotf will
be zero.
#include <math.h>
unsigned int ta = 0x80000002;
main ()
{
float a = *(float *)&ta;
float c = hypotf(a, a);
float d = hypot(a, a);
printf("a = %e hypotf(a,a) = %e hypot(a,a) = %e\n", a, c, d);
}
The attached patch combines my fix for denormal values and Richard Sandiford's
change from newlib. Presumably he had a good reason for changing the scaling
factor to 2^68. (I'm including Richard in the ChangeLog entry because I'm using
his newlib patch.)
2005-08-01 Bob Wilson <bob.wilson@acm.org>
Richard Sandiford <richard@codesourcery.com>
* sysdeps/ieee754/flt-32/e_hypotf.c (__ieee754_hypotf): Add missing
exponent bias to the value for 2^126. Change 2^60 scaling factor to
2^68 and remove the exponent bias from the increment value.
Index: sysdeps/ieee754/flt-32/e_hypotf.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/flt-32/e_hypotf.c,v
retrieving revision 1.1
diff -u -r1.1 e_hypotf.c
--- sysdeps/ieee754/flt-32/e_hypotf.c 13 Jul 1999 23:58:13 -0000 1.1
+++ sysdeps/ieee754/flt-32/e_hypotf.c 1 Aug 2005 20:40:51 -0000
@@ -46,22 +46,22 @@
if(hb == 0x7f800000) w = b;
return w;
}
- /* scale a and b by 2**-60 */
- ha -= 0x5d800000; hb -= 0x5d800000; k += 60;
+ /* scale a and b by 2**-68 */
+ ha -= 0x22000000; hb -= 0x22000000; k += 68;
SET_FLOAT_WORD(a,ha);
SET_FLOAT_WORD(b,hb);
}
if(hb < 0x26800000) { /* b < 2**-50 */
if(hb <= 0x007fffff) { /* subnormal b or 0 */
if(hb==0) return a;
- SET_FLOAT_WORD(t1,0x3f000000); /* t1=2^126 */
+ SET_FLOAT_WORD(t1,0x7e800000); /* t1=2^126 */
b *= t1;
a *= t1;
k -= 126;
- } else { /* scale a and b by 2^60 */
- ha += 0x5d800000; /* a *= 2^60 */
- hb += 0x5d800000; /* b *= 2^60 */
- k -= 60;
+ } else { /* scale a and b by 2^68 */
+ ha += 0x22000000; /* a *= 2^68 */
+ hb += 0x22000000; /* b *= 2^68 */
+ k -= 68;
SET_FLOAT_WORD(a,ha);
SET_FLOAT_WORD(b,hb);
}