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
This commit is contained in:
Bruce Evans 1997-02-20 12:37:49 +00:00
parent d5e4d7e154
commit 487ddb8fa4
2 changed files with 23 additions and 11 deletions

View File

@ -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

View File

@ -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