This is the mail archive of the
libc-alpha@sources.redhat.com
mailing list for the glibc project.
[PATCH] fix for accuracy problem in atan2f
- From: Bob Wilson <bwilson at tensilica dot com>
- To: libc-alpha at sources dot redhat dot com
- Date: Wed, 20 Jul 2005 08:24:20 -0700
- Subject: [PATCH] fix for accuracy problem in atan2f
There is a problem with the version of atan2f in
sysdeps/ieee754/flt-32/e_atan2f.c that appears to be a typo in the value
of pi_lo, and that causes the results to be slightly inaccurate for some
input values. The attached test program will demonstrate the problem on
a platform that uses this version of atan2f. (I don't see the problem
on i386, presumably because it's using a machine-specific version of
atan2f, but I do see it on x86_64 with RedHat Enterprise 3.)
Back in July 1997, the value of pi in this file (back when it was in
sysdeps/libm-ieee754/) was changed to be the correct round-to-nearest
value, and pi_lo was adjusted accordingly. It appears that the exponent
of pi_lo was mistakenly left unchanged from its previous value. There
are several indications of this:
- The hex value of pi_lo shown in comments is 0xb3bbbd2e, which
corresponds to -8.7422776573e-08, but the float value actually used has
the same mantissa but an exponent of "e-07".
- The values of pi and pi_lo should add up to the real value of pi (as
much as can be represented). The exponent should be "e-08" for this to
be true.
- The test case below gets the wrong answer (atan2f != (float)atan2)
when pi_lo has an exponent of "e-07" but gets the correct answer when
the exponent is changed to "e-08".
The testcase and patch are attached below. Here is a ChangeLog entry:
2005-07-20 Bob Wilson <bob.wilson@acm.org>
Darin Petkov <darin@tensilica.com>
* sysdeps/ieee754/flt-32/e_atan2f.c (pi_lo): Correct exponent value.
#include <math.h>
unsigned int ta = 0x00800000;
unsigned int tb = 0x80800000;
unsigned int tc = 0x4016CBE0;
unsigned int td = 0x4016CBE4;
main ()
{
float a = *(float *)&ta;
float b = *(float *)&tb;
printf("arguments: a = %e b = %e\n", a, b);
printf("atan2f(a, b) = %.13e\n", atan2f(a, b));
printf("atan2(a, b) = %.13e\n", (float) atan2(a, b));
printf("wrong value = %.13e\n", *(float *)&tc);
printf("right value = %.13e\n", *(float *)&td);
}
Index: e_atan2f.c
===================================================================
RCS file: /cvs/glibc/libc/sysdeps/ieee754/flt-32/e_atan2f.c,v
retrieving revision 1.1
diff -u -p -r1.1 e_atan2f.c
--- e_atan2f.c 13 Jul 1999 23:57:45 -0000 1.1
+++ e_atan2f.c 20 Jul 2005 13:38:44 -0000
@@ -30,7 +30,7 @@ zero = 0.0,
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415927410e+00, /* 0x40490fdb */
-pi_lo = -8.7422776573e-07; /* 0xb3bbbd2e */
+pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
#ifdef __STDC__
float __ieee754_atan2f(float y, float x)