Fix x86 acos near 1 (bug 13942).

diff --git a/ChangeLog b/ChangeLog
index f436af4083e307042dac9ebcf54f97cfa216facd..a13a367b29fed4fef8536566f7c1237a20b9af56 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,15 @@
 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.
+
        [BZ #14034]
        * sysdeps/i386/fpu/e_acos.S (__ieee754_acos): Take absolute value
        of square root.

diff --git a/NEWS b/NEWS
index 4481ea43fc007ef4797d67afa297114e2f914416..c9848e68e529e1092a09ca001ae1146e07883082 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -22,7 +22,7 @@ Version 2.16
   13844, 13846, 13851, 13852, 13854, 13871, 13872, 13873, 13879, 13883,
   13886, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915, 13916,
   13917, 13918, 13919, 13920, 13921, 13924, 13926, 13927, 13928, 13938,
-  13941, 13963, 13967, 13970, 13973, 14027, 14033, 14034
+  13941, 13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034
 
 * ISO C11 support:
 

diff --git a/math/libm-test.inc b/math/libm-test.inc
index bedff09fea5de3f58ae3196429945f2abd159f14..d643badfe23c55231a6d74627e21bc09b5f63b6a 100644 (file)
--- 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 d10a054b9cc4b0cc32d285ff705d3cc7c5d0e866..586c7fc406b1867f62379d25c748bb206606aba0 100644 (file)
--- a/sysdeps/i386/fpu/e_acos.S
+++ b/sysdeps/i386/fpu/e_acos.S
@@ -7,13 +7,15 @@
 
 RCSID("$NetBSD: e_acos.S,v 1.4 1995/05/08 23:44:37 jtc Exp $")
 
-/* acos = atan (sqrt(1 - x^2) / x) */
+/* acos = atan (sqrt((1-x) (1+x)) / x) */
 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 d249d5a956041c8d1a2a26800cdb99c9a999e3ab..ab0893192490923b471a83733474b2f89c5ee323 100644 (file)
--- a/sysdeps/i386/fpu/e_acosl.c
+++ b/sysdeps/i386/fpu/e_acosl.c
@@ -12,11 +12,13 @@ __ieee754_acosl (long double x)
 {
   long double res;
 
-  /* acosl = atanl (sqrtl(1 - x^2) / x) */
+  /* acosl = atanl (sqrtl((1-x) (1+x)) / 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 a17e922b6dfa1f14a5fd7671dcff7bdf8dcfced7..9a44cb62ac413d9cee7f12a66aa0b7ab18544dd6 100644 (file)
--- a/sysdeps/i386/fpu/e_asin.S
+++ b/sysdeps/i386/fpu/e_asin.S
@@ -7,13 +7,15 @@
 
 RCSID("$NetBSD: e_asin.S,v 1.4 1995/05/08 23:45:40 jtc Exp $")
 
-/* asin = atan (x / sqrt(1 - x^2)) */
+/* asin = atan (x / sqrt((1-x) (1+x))) */
 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 fdaff3525d6267386f3f7a3a87bc85c56937fedb..a49a69d4e98ce67825c861510429837923f02736 100644 (file)
--- 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 5068fe695266ca7f3b17d2cbaeae4eb5a39096b6..9a1e80ed8da9c2d4848d5fc5c3696513547a828d 100644 (file)
--- 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