[Bug math/13942] New: x86 acos inaccurate near 1

[Bug math/13942] New: x86 acos inaccurate near 1

[jsm28 at gcc dot gnu.org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

             Bug #: 13942
           Summary: x86 acos inaccurate near 1
           Product: glibc
           Version: 2.15
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: jsm28@gcc.gnu.org
    Classification: Unclassified


http://sourceware.org/ml/libc-help/2008-09/msg00031.html (mentioned in a recent
bug) reports an issue with acos on x86 being inaccurate near 1 because of using
sqrt (1 - x*x).  This is still present with current sources, as illustrated by
the test:

#include <math.h>
#include <stdio.h>

volatile double d = 0x0.ffffffff8p0;

int
main (void)
{
  volatile double r = acos (d);
  printf ("%.14a\n", r);
  return 0;
}

This prints 0x1.000000002aaab0p-16 when an accurate result computed by GNU MPFR
is 0x1.000000000aaab0p-16.

It looks like this issue will affect x86 acos and acosl and x86_64 acosl.  For
acosf the squaring (in x87 64-bit precision) will always be exact so there is
no problem.  For asin functions and for acos functions near -1 I don't think
there will be significant loss of accuracy from this issue (the result will be
near pi/2 or pi and the problem is for results near 0) but probably all the
functions (other than asinf/acosf) should be fixed and tests added to the
testsuite for all these cases.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[bugdal at aerifal dot cx]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

Rich Felker <bugdal at aerifal dot cx> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |bugdal at aerifal dot cx

--- Comment #1 from Rich Felker <bugdal at aerifal dot cx> 2012-04-04 01:16:15 UTC ---
We experienced the same issue in musl's implementation, and fixed it by using
the following formula for acos in the asm:

acos(x) = 2*atan2(sqrt(1-x), sqrt(1+x))

Assuming the right hand side is evaluated in extended precision, it seems to
give <1ulp error for double precision inputs. We were doing the division
(rather than relying on the FPU's 2-argument atan) so we could avoid one sqrt,
but it was causing incorrect results with downward rounding direction on some
corner cases, so we switched to using two sqrts.

asin does not seem to suffer from precision issues with the usual formula
involving x².

Here's our asm:

http://git.etalabs.net/cgi-bin/gitweb.cgi?p=musl;a=blob;f=src/math/i386/acos.s;h=bfff0c5c9a89d3332ea7df6a8e5d3cd0a68a01f3;hb=HEAD

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #2 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-04-04 09:20:34 UTC ---
I expect the approach suggested in the libc-help posting of using 
(1-x)*(1+x) would be fine as a fix.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[bugdal at aerifal dot cx]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #3 from Rich Felker <bugdal at aerifal dot cx> 2012-04-04 14:13:00 UTC ---
Hi Joseph,

I did some preliminary tests, and that formula does not seem to fix it.

Given the following, which has been written in C to use extended precision for
intermediate results:

  double x = 0x0.ffffffffffffp0;
  printf("%.17g\n", acos(x));
  printf("%.17g\n", (double)atanl(sqrtl(1-(long double)x*x)));
  printf("%.17g\n", (double)atanl(sqrtl((1-(long double)x)*(1+(long
double)x))));

I'm getting:

8.4293697021788083e-08
8.4293697021787858e-08
8.4293697021787791e-08

The first result (our current implementation in musl) agrees with Wolfram
Alpha's value for acos(0.999999999999996447286321199499070644378662109375).
Both of the latter results have relatively huge error. So unless I screwed up
something in doing the intermediate results at extended precision, I don't
think your solution works.

Note that our algorithm avoids performing any multiplication or division that
would lose precision. The only loss of precision takes place at the sqrt
operations (which will be correctly rounded for extended precision) and the
atan operation.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #4 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-04-04 14:24:05 UTC ---
On Wed, 4 Apr 2012, bugdal at aerifal dot cx wrote:

>   printf("%.17g\n", (double)atanl(sqrtl(1-(long double)x*x)));
>   printf("%.17g\n", (double)atanl(sqrtl((1-(long double)x)*(1+(long
> double)x))));

That's not the formula.  It's atan2l (sqrtl (...), x) that the code is 
effectively using (and the point is to make what's inside the sqrtl more 
accurate).

> Note that our algorithm avoids performing any multiplication or division that
> would lose precision. The only loss of precision takes place at the sqrt
> operations (which will be correctly rounded for extended precision) and the
> atan operation.

Depending on the input, 1+x or 1-x may lose precision as well for the long 
double functions such as acosl, but I don't think that loss is 
significant for the final result.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[bugdal at aerifal dot cx]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #5 from Rich Felker <bugdal at aerifal dot cx> 2012-04-04 23:49:46 UTC ---
Indeed, I screwed up that test. It seems to work with the error corrected,
though I have not done any rigorous testing. One member of our development team
commented however that your version is expected to give wrong results for x=+-1
in downward rounding mode.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #6 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-04-05 00:12:54 UTC ---
fabs instructions could be used to avoid issues with sqrt(-0) = -0 (if 
that's the round-downward issue you mean) and may be cheaper than taking 
two square roots.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[bugdal at aerifal dot cx]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #7 from Rich Felker <bugdal at aerifal dot cx> 2012-04-05 00:30:44 UTC ---
Already considered that. Using fabs will make the function accept inputs
outside the domain without returning NAN.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #8 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-04-05 09:33:32 UTC ---
On Thu, 5 Apr 2012, bugdal at aerifal dot cx wrote:

> Already considered that. Using fabs will make the function accept inputs
> outside the domain without returning NAN.

I'm thinking of fabs after taking the square root, not before.  (And 
normally - except for the _LIB_VERSION == _IEEE_ case - the wrappers such 
as math/w_acos.c will have handled arguments outside the domain before 
this code is called at all.)

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[bugdal at aerifal dot cx]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

--- Comment #9 from Rich Felker <bugdal at aerifal dot cx> 2012-04-05 14:13:34 UTC ---
We were reportedly getting some >1ulp errors with your version too, but can't
reproduce them at the moment. I'll keep you posted if we find anything, but it
looks like your version (with the abs after the sqrt) is likely the best.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[jsm28 at gcc dot gnu.org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13942

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |FIXED

--- Comment #10 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-04-30 18:58:21 UTC ---
Fixed by:

commit adfbc8ac9e192b6e3007f7a47852df937afa2145
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Mon Apr 30 18:56:39 2012 +0000

    Fix x86 acos near 1 (bug 13942).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13942] x86 acos inaccurate near 1

[fweimer at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=13942

Florian Weimer <fweimer at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
              Flags|                            |security-

-- 
You are receiving this mail because:
You are on the CC list for the bug.