Compute (1 - x^2) as ((1 - x) * (1 + x)) instead of as (1 - x * x) to

avoid easily avoidable loss of precision when |x| is nearly 1.

Extended (64-bit) precision only moves the meaning of "nearly" here.

This probably could be done better by splitting up the range into
|x| <= 0.5 and |x| > 0.5 like the C version.  However, ucbtest
does't report any errors in this version.  Perhaps the C version
should be used anyway.  It's only 25% slower now on a P5, provided
the C version of sqrt() isn't used, and the C version could be
optimized better.

Errors checked by:	ucbtest

2 files changed, 23 insertions, 11 deletions

diff --git a/lib/msun/i387/e_acos.S b/lib/msun/i387/e_acos.S
index 84cb542..774c659 100644
--- a/lib/msun/i387/e_acos.S
+++ b/lib/msun/i387/e_acos.S
@@ -37,14 +37,20 @@
 
 RCSID("$FreeBSD$")
 
-/* acos = atan (sqrt(1 - x^2) / x) */
+/*
+ * acos(x) = atan2(sqrt(1  - x^2, x).
+ * Actually evaluate (1 - x^2) as (1 - x) * (1 + x) to avoid loss of
+ * precision when |x| is nearly 1.
+ */
 ENTRY(__ieee754_acos)
 	fldl	4(%esp)			/* x */
-	fst	%st(1)
-	fmul	%st(0)			/* x^2 */
-	fld1				
-	fsubp				/* 1 - x^2 */
-	fsqrt				/* sqrt (1 - x^2) */
+	fld1
+	fld	%st(0)
+	fsub	%st(2)			/* 1 - x */
+	fxch	%st(1)
+	fadd	%st(2)			/* 1 + x */
+	fmulp	%st(1)			/* (1 - x) * (1 + x) */
+	fsqrt
 	fxch	%st(1)
 	fpatan
 	ret
diff --git a/lib/msun/i387/e_asin.S b/lib/msun/i387/e_asin.S
index 5b2a1bd..de031cf 100644
--- a/lib/msun/i387/e_asin.S
+++ b/lib/msun/i387/e_asin.S
@@ -37,13 +37,19 @@
 
 RCSID("$FreeBSD$")
 
-/* asin = atan (x / sqrt(1 - x^2)) */
+/*
+ * asin(x) = atan2(x, sqrt(1  - x^2).
+ * Actually evaluate (1 - x^2) as (1 - x) * (1 + x) to avoid loss of
+ * precision when |x| is nearly 1.
+ */
 ENTRY(__ieee754_asin)
 	fldl	4(%esp)			/* x */
-	fst	%st(1)
-	fmul	%st(0)			/* x^2 */
 	fld1
-	fsubp				/* 1 - x^2 */
-	fsqrt				/* sqrt (1 - x^2) */
+	fld	%st(0)
+	fsub	%st(2)			/* 1 - x */
+	fxch	%st(1)
+	fadd	%st(2)			/* 1 + x */
+	fmulp	%st(1)			/* (1 - x) * (1 + x) */
+	fsqrt
 	fpatan
 	ret