summaryrefslogtreecommitdiffstats
path: root/lib
diff options
context:
space:
mode:
authorbde <bde@FreeBSD.org>1997-02-20 12:37:49 +0000
committerbde <bde@FreeBSD.org>1997-02-20 12:37:49 +0000
commitd9c878777ffccc6ea50f3a143bb94d5a827fb9fd (patch)
treefdffc89d06b18a6465de3e1b9e67ae8a69046ecc /lib
parent301b9cec8851010b8f538d86f85f2a2f2b266297 (diff)
downloadFreeBSD-src-d9c878777ffccc6ea50f3a143bb94d5a827fb9fd.zip
FreeBSD-src-d9c878777ffccc6ea50f3a143bb94d5a827fb9fd.tar.gz
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
Diffstat (limited to 'lib')
-rw-r--r--lib/msun/i387/e_acos.S18
-rw-r--r--lib/msun/i387/e_asin.S16
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
OpenPOWER on IntegriCloud