Fix x86 acos near 1 (bug 13942)

Mon Apr 30 19:20:00 GMT 2012

On 04/30/2012 02:11 PM, Joseph S. Myers wrote:
> Bug 13942 is the problem of acos on x86, and acosl on x86/x86_64,
> being very inaccurate for arguments near 1 because of a naive
> calculation of 1 - x*x, which can be calculated more accurately as (1
> - x) * (1 + x).
>
> I propose this patch to fix this.  I also fixed the corresponding
> calculation in asin, although that only makes a small accuracy
> improvement there (2ulp error changes to 0ulp).  acosf and asinf are
> unaffected, since the squaring is exact there, and asinl uses a
> different (C) implementation, not x87 assembly.  Note that the tests
> added would suffice to catch any such problems for the wider long
> double implementations, since squaring a 64-bit value is substantially
> inaccurate even for 113-bit mantissa, so there is no need for further
> tests specifically for wider long double formats.
>
> Tested x86 and x86_64 and ulps updated accordingly.
>
> 2012-04-30  Joseph Myers<joseph@codesourcery.com>
>
> 	[BZ #13942]
> 	* sysdeps/i386/fpu/e_acos.S (__ieee754_acos): Calculate 1 - x^2 as
> 	(1 - x) * (1 + x).
> 	* sysdeps/i386/fpu/e_acosl.c (__ieee754_acosl): Likewise.
> 	* sysdeps/i386/fpu/e_asin.S (__ieee754_asin): Likewise.
> 	* math/libm-test.inc (acos_test): Add more tests.
> 	(asin_test): Likewise.
> 	* sysdeps/i386/fpu/libm-test-ulps: Update.
> 	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
>
> diff --git a/math/libm-test.inc b/math/libm-test.inc
> index bedff09..d643bad 100644
> --- a/math/libm-test.inc
> +++ b/math/libm-test.inc
> @@ -758,6 +758,18 @@ acos_test (void)
>     TEST_f_f (acos, 0.75L, 0.722734247813415611178377352641333362L);
>     TEST_f_f (acos, 2e-17L, 1.57079632679489659923132169163975144L);
>     TEST_f_f (acos, 0.0625L, 1.50825556499840522843072005474337068L);
> +  TEST_f_f (acos, 0x0.ffffffp0L, 3.4526698471620358760324948263873649728491e-4L);
> +  TEST_f_f (acos, -0x0.ffffffp0L, 3.1412473866050770348750401337968641476999L);
> +#ifndef TEST_FLOAT
> +  TEST_f_f (acos, 0x0.ffffffff8p0L, 1.5258789062648029736620564947844627548516e-5L);
> +  TEST_f_f (acos, -0x0.ffffffff8p0L, 3.1415773948007305904329067627145550395696L);
> +  TEST_f_f (acos, 0x0.ffffffffffffp0L, 8.4293697021788088529885473244391795127130e-8L);
> +  TEST_f_f (acos, -0x0.ffffffffffffp0L, 3.1415925692960962166745548533940296398054L);
> +#endif
> +#if defined TEST_LDOUBLE&&  LDBL_MANT_DIG>= 64
> +  TEST_f_f (acos, 0x0.ffffffffffffffffp0L, 3.2927225399135962333718255320079907245059e-10L);
> +  TEST_f_f (acos, -0x0.ffffffffffffffffp0L, 3.1415926532605209844712837599423203309964L);
> +#endif
>     END (acos);
>   }
>
> @@ -933,6 +945,18 @@ asin_test (void)
>     TEST_f_f (asin, 1.0, M_PI_2l);
>     TEST_f_f (asin, -1.0, -M_PI_2l);
>     TEST_f_f (asin, 0.75L, 0.848062078981481008052944338998418080L);
> +  TEST_f_f (asin, 0x0.ffffffp0L, 1.5704510598101804156437184421571127056013L);
> +  TEST_f_f (asin, -0x0.ffffffp0L, -1.5704510598101804156437184421571127056013L);
> +#ifndef TEST_FLOAT
> +  TEST_f_f (asin, 0x0.ffffffff8p0L, 1.5707810680058339712015850710748035974710L);
> +  TEST_f_f (asin, -0x0.ffffffff8p0L, -1.5707810680058339712015850710748035974710L);
> +  TEST_f_f (asin, 0x0.ffffffffffffp0L, 1.5707962425011995974432331617542781977068L);
> +  TEST_f_f (asin, -0x0.ffffffffffffp0L, -1.5707962425011995974432331617542781977068L);
> +#endif
> +#if defined TEST_LDOUBLE&&  LDBL_MANT_DIG>= 64
> +  TEST_f_f (asin, 0x0.ffffffffffffffffp0L, 1.5707963264656243652399620683025688888978L);
> +  TEST_f_f (asin, -0x0.ffffffffffffffffp0L, -1.5707963264656243652399620683025688888978L);
> +#endif
>
>     END (asin);
>   }
> diff --git a/sysdeps/i386/fpu/e_acos.S b/sysdeps/i386/fpu/e_acos.S
> index d10a054..fc680ee 100644
> --- a/sysdeps/i386/fpu/e_acos.S
> +++ b/sysdeps/i386/fpu/e_acos.S
> @@ -11,9 +11,11 @@ RCSID("$NetBSD: e_acos.S,v 1.4 1995/05/08 23:44:37 jtc Exp $")
>   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 d249d5a..42c5f01 100644
> --- a/sysdeps/i386/fpu/e_acosl.c
> +++ b/sysdeps/i386/fpu/e_acosl.c
> @@ -14,9 +14,11 @@ __ieee754_acosl (long double x)
>
>     /* acosl = atanl (sqrtl(1 - x^2) / x) */

I suggest to enhance the comment above with
= atanl (sqrtl (((1-x) (1+x))/ x)

Other files have the same comment, please add it there as well.

Ok with those changes,

Andreas

>     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 a17e922..b4d17b0 100644
> --- a/sysdeps/i386/fpu/e_asin.S
> +++ b/sysdeps/i386/fpu/e_asin.S
> @@ -11,9 +11,11 @@ RCSID("$NetBSD: e_asin.S,v 1.4 1995/05/08 23:45:40 jtc Exp $")
>   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

-- 
  Andreas Jaeger aj@{suse.com,opensuse.org} Twitter/Identica: jaegerandi
   SUSE LINUX Products GmbH, Maxfeldstr. 5, 90409 NÃ¼rnberg, Germany
    GF: Jeff Hawn,Jennifer Guild,Felix ImendÃ¶rffer,HRB16746 (AG NÃ¼rnberg)
     GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126