From adfbc8ac9e192b6e3007f7a47852df937afa2145 Mon Sep 17 00:00:00 2001
From: Joseph Myers <joseph@codesourcery.com>
Date: Mon, 30 Apr 2012 18:56:39 +0000
Subject: Fix x86 acos near 1 (bug 13942).

---
 sysdeps/i386/fpu/e_acos.S       | 10 ++++++----
 sysdeps/i386/fpu/e_acosl.c      |  8 +++++---
 sysdeps/i386/fpu/e_asin.S       | 10 ++++++----
 sysdeps/i386/fpu/libm-test-ulps | 12 ++++++++++++
 4 files changed, 29 insertions(+), 11 deletions(-)

(limited to 'sysdeps/i386')

diff --git a/sysdeps/i386/fpu/e_acos.S b/sysdeps/i386/fpu/e_acos.S
index d10a054b9c..586c7fc406 100644
--- 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 d249d5a956..ab08931924 100644
--- 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 a17e922b6d..9a44cb62ac 100644
--- 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 fdaff3525d..a49a69d4e9 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
-- 
cgit 1.4.1