This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix x86 acos near 1 (bug 13942)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Mon, 30 Apr 2012 12:11:41 +0000 (UTC)
- Subject: Fix x86 acos near 1 (bug 13942)
Bug 13942 is the problem of acos on x86, and acosl on x86/x86_64,
being very inaccurate for arguments near 1 because of a naive
calculation of 1 - x*x, which can be calculated more accurately as (1
- x) * (1 + x).
I propose this patch to fix this. I also fixed the corresponding
calculation in asin, although that only makes a small accuracy
improvement there (2ulp error changes to 0ulp). acosf and asinf are
unaffected, since the squaring is exact there, and asinl uses a
different (C) implementation, not x87 assembly. Note that the tests
added would suffice to catch any such problems for the wider long
double implementations, since squaring a 64-bit value is substantially
inaccurate even for 113-bit mantissa, so there is no need for further
tests specifically for wider long double formats.
Tested x86 and x86_64 and ulps updated accordingly.
2012-04-30 Joseph Myers <joseph@codesourcery.com>
[BZ #13942]
* sysdeps/i386/fpu/e_acos.S (__ieee754_acos): Calculate 1 - x^2 as
(1 - x) * (1 + x).
* sysdeps/i386/fpu/e_acosl.c (__ieee754_acosl): Likewise.
* sysdeps/i386/fpu/e_asin.S (__ieee754_asin): Likewise.
* math/libm-test.inc (acos_test): Add more tests.
(asin_test): Likewise.
* sysdeps/i386/fpu/libm-test-ulps: Update.
* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index bedff09..d643bad 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -758,6 +758,18 @@ acos_test (void)
TEST_f_f (acos, 0.75L, 0.722734247813415611178377352641333362L);
TEST_f_f (acos, 2e-17L, 1.57079632679489659923132169163975144L);
TEST_f_f (acos, 0.0625L, 1.50825556499840522843072005474337068L);
+ TEST_f_f (acos, 0x0.ffffffp0L, 3.4526698471620358760324948263873649728491e-4L);
+ TEST_f_f (acos, -0x0.ffffffp0L, 3.1412473866050770348750401337968641476999L);
+#ifndef TEST_FLOAT
+ TEST_f_f (acos, 0x0.ffffffff8p0L, 1.5258789062648029736620564947844627548516e-5L);
+ TEST_f_f (acos, -0x0.ffffffff8p0L, 3.1415773948007305904329067627145550395696L);
+ TEST_f_f (acos, 0x0.ffffffffffffp0L, 8.4293697021788088529885473244391795127130e-8L);
+ TEST_f_f (acos, -0x0.ffffffffffffp0L, 3.1415925692960962166745548533940296398054L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 64
+ TEST_f_f (acos, 0x0.ffffffffffffffffp0L, 3.2927225399135962333718255320079907245059e-10L);
+ TEST_f_f (acos, -0x0.ffffffffffffffffp0L, 3.1415926532605209844712837599423203309964L);
+#endif
END (acos);
}
@@ -933,6 +945,18 @@ asin_test (void)
TEST_f_f (asin, 1.0, M_PI_2l);
TEST_f_f (asin, -1.0, -M_PI_2l);
TEST_f_f (asin, 0.75L, 0.848062078981481008052944338998418080L);
+ TEST_f_f (asin, 0x0.ffffffp0L, 1.5704510598101804156437184421571127056013L);
+ TEST_f_f (asin, -0x0.ffffffp0L, -1.5704510598101804156437184421571127056013L);
+#ifndef TEST_FLOAT
+ TEST_f_f (asin, 0x0.ffffffff8p0L, 1.5707810680058339712015850710748035974710L);
+ TEST_f_f (asin, -0x0.ffffffff8p0L, -1.5707810680058339712015850710748035974710L);
+ TEST_f_f (asin, 0x0.ffffffffffffp0L, 1.5707962425011995974432331617542781977068L);
+ TEST_f_f (asin, -0x0.ffffffffffffp0L, -1.5707962425011995974432331617542781977068L);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MANT_DIG >= 64
+ TEST_f_f (asin, 0x0.ffffffffffffffffp0L, 1.5707963264656243652399620683025688888978L);
+ TEST_f_f (asin, -0x0.ffffffffffffffffp0L, -1.5707963264656243652399620683025688888978L);
+#endif
END (asin);
}
diff --git a/sysdeps/i386/fpu/e_acos.S b/sysdeps/i386/fpu/e_acos.S
index d10a054..fc680ee 100644
--- a/sysdeps/i386/fpu/e_acos.S
+++ b/sysdeps/i386/fpu/e_acos.S
@@ -11,9 +11,11 @@ RCSID("$NetBSD: e_acos.S,v 1.4 1995/05/08 23:44:37 jtc Exp $")
ENTRY(__ieee754_acos)
fldl 4(%esp) /* x */
fld %st /* x : x */
- fmul %st(0) /* x^2 : x */
- fld1 /* 1 : x^2 : x */
- fsubp /* 1 - x^2 : x */
+ fld1 /* 1 : x : x */
+ fsubp /* 1 - x : x */
+ fld1 /* 1 : 1 - x : x */
+ fadd %st(2) /* 1 + x : 1 - x : x */
+ fmulp /* 1 - x^2 : x */
fsqrt /* sqrt (1 - x^2) : x */
fabs
fxch %st(1) /* x : sqrt (1 - x^2) */
diff --git a/sysdeps/i386/fpu/e_acosl.c b/sysdeps/i386/fpu/e_acosl.c
index d249d5a..42c5f01 100644
--- a/sysdeps/i386/fpu/e_acosl.c
+++ b/sysdeps/i386/fpu/e_acosl.c
@@ -14,9 +14,11 @@ __ieee754_acosl (long double x)
/* acosl = atanl (sqrtl(1 - x^2) / x) */
asm ( "fld %%st\n"
- "fmul %%st(0)\n" /* x^2 */
"fld1\n"
- "fsubp\n" /* 1 - x^2 */
+ "fsubp\n"
+ "fld1\n"
+ "fadd %%st(2)\n"
+ "fmulp\n" /* 1 - x^2 */
"fsqrt\n" /* sqrtl (1 - x^2) */
"fabs\n"
"fxch %%st(1)\n"
diff --git a/sysdeps/i386/fpu/e_asin.S b/sysdeps/i386/fpu/e_asin.S
index a17e922..b4d17b0 100644
--- a/sysdeps/i386/fpu/e_asin.S
+++ b/sysdeps/i386/fpu/e_asin.S
@@ -11,9 +11,11 @@ RCSID("$NetBSD: e_asin.S,v 1.4 1995/05/08 23:45:40 jtc Exp $")
ENTRY(__ieee754_asin)
fldl 4(%esp) /* x */
fld %st
- fmul %st(0) /* x^2 */
- fld1
- fsubp /* 1 - x^2 */
+ fld1 /* 1 : x : x */
+ fsubp /* 1 - x : x */
+ fld1 /* 1 : 1 - x : x */
+ fadd %st(2) /* 1 + x : 1 - x : x */
+ fmulp /* 1 - x^2 */
fsqrt /* sqrt (1 - x^2) */
fpatan
ret
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index fdaff35..a49a69d 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -80,6 +80,12 @@ ldouble: 1
Test "asin (-0.5) == -pi/6":
ildouble: 1
ldouble: 1
+Test "asin (-0x0.ffffffff8p0) == -1.5707810680058339712015850710748035974710":
+ildouble: 1
+ldouble: 1
+Test "asin (-0x0.ffffffp0) == -1.5704510598101804156437184421571127056013":
+ildouble: 1
+ldouble: 1
Test "asin (-1.0) == -pi/2":
ildouble: 1
ldouble: 1
@@ -89,6 +95,12 @@ ldouble: 1
Test "asin (0.75) == 0.848062078981481008052944338998418080":
ildouble: 1
ldouble: 1
+Test "asin (0x0.ffffffff8p0) == 1.5707810680058339712015850710748035974710":
+ildouble: 1
+ldouble: 1
+Test "asin (0x0.ffffffp0) == 1.5704510598101804156437184421571127056013":
+ildouble: 1
+ldouble: 1
Test "asin (1.0) == pi/2":
ildouble: 1
ldouble: 1
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 5068fe6..9a1e80e 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -55,6 +55,12 @@ ldouble: 1
Test "asin (-0.5) == -pi/6":
ildouble: 1
ldouble: 1
+Test "asin (-0x0.ffffffffffffffffp0) == -1.5707963264656243652399620683025688888978":
+ildouble: 1
+ldouble: 1
+Test "asin (-0x0.ffffffffffffp0) == -1.5707962425011995974432331617542781977068":
+ildouble: 1
+ldouble: 1
Test "asin (-1.0) == -pi/2":
ildouble: 1
ldouble: 1
@@ -64,6 +70,12 @@ ldouble: 1
Test "asin (0.75) == 0.848062078981481008052944338998418080":
ildouble: 1
ldouble: 1
+Test "asin (0x0.ffffffffffffffffp0) == 1.5707963264656243652399620683025688888978":
+ildouble: 1
+ldouble: 1
+Test "asin (0x0.ffffffffffffp0) == 1.5707962425011995974432331617542781977068":
+ildouble: 1
+ldouble: 1
Test "asin (1.0) == pi/2":
ildouble: 1
ldouble: 1
--
Joseph S. Myers
joseph@codesourcery.com