From d9c878777ffccc6ea50f3a143bb94d5a827fb9fd Mon Sep 17 00:00:00 2001 From: bde Date: Thu, 20 Feb 1997 12:37:49 +0000 Subject: 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 --- lib/msun/i387/e_acos.S | 18 ++++++++++++------ lib/msun/i387/e_asin.S | 16 +++++++++++----- 2 files changed, 23 insertions(+), 11 deletions(-) (limited to 'lib/msun/i387') 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 -- cgit v1.1