Bug 13381 - rounding problem wih sincosl
Summary: rounding problem wih sincosl
Status: RESOLVED INVALID
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.14
: P2 normal
Target Milestone: ---
Assignee: Andreas Jaeger
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2011-11-04 12:50 UTC by gastineau
Modified: 2014-06-27 11:39 UTC (History)
4 users (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description gastineau 2011-11-04 12:50:55 UTC
With some input values, the sincosl function and sinl returns a different value for the sin part. 

This problem occurs with glibc 2.14.

The code source is at the end of the message.

The output with these compilers is :
gcc -m64 bugsincos.c -O0 -o bugsincos -lm
./bugsincos
diffC=0.000000e+00 C1=0xf.e8f0776670a18eap-4 C=0xf.e8f0776670a18eap-4
diffS=-6.776264e-21 S1=0xd.903ba52345a81a4p-7 S=0xd.903ba52345a81a5p-7

diffS must be 0 if same algorithm is used inside the math library.

If I compile with "-O3", the problem doesn't occur because it doesn't use the math library.
I give some system information below and the source code.

Mickaël

xxxxx$ gcc  -m64 bugsincos.c -O0 -o bugsincos -lm
xxxxx$ ./bugsincos
diffC=0.000000e+00 C1=0xf.e8f0776670a18eap-4 C=0xf.e8f0776670a18eap-4
diffS=-6.776264e-21 S1=0xd.903ba52345a81a4p-7 S=0xd.903ba52345a81a5p-7
xxxxx$ nm bugsincos | grep sincosl
                 U sincosl@@GLIBC_2.2.5
xxxxx$ rpm -qa | grep glibc
glibc-static-2.14-5.x86_64
glibc-headers-2.14-5.x86_64
glibc-devel-2.14-5.x86_64
glibc-2.14-5.i686
glibc-common-2.14-5.x86_64
glibc-2.14-5.x86_64
xxxxx$ gcc --version
gcc (GCC) 4.6.1 20110908 (Red Hat 4.6.1-9)
Copyright (C) 2011 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


-------------- cut here----------------
#include <mathimf.h>
#include <stdio.h>

int main()
{
long double X=0xd.96c3941c1094ed5p-7L;
long double C=cosl(X);
long double S=sinl(X);
long double C1, S1;
sincosl(X, &S1, &C1);
printf("diffC=%Le C1=%La C=%La\n", C1-C, C1, C);
printf("diffS=%Le S1=%La S=%La\n", S1-S, S1, S);
return 0;

}
-------------- cut here----------------
Comment 1 gastineau 2011-11-07 16:17:24 UTC
The correct source code is :


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

int main()
{
long double X=0xd.96c3941c1094ed5p-7L;
long double C=cosl(X);
long double S=sinl(X);
long double C1, S1;
sincosl(X, &S1, &C1);
printf("diffC=%Le C1=%La C=%La\n", C1-C, C1, C);
printf("diffS=%Le S1=%La S=%La\n", S1-S, S1, S);
return 0;

}
Comment 2 Andreas Schwab 2011-11-08 07:47:15 UTC
You'll have to complain to Intel.
Comment 3 gastineau 2011-11-08 08:36:20 UTC
I don't understand. I use the gcc compiler and the glibc. 
You should check the second source code.  The first source code was an error.
 The second source code produce the error.
Comment 4 Andreas Schwab 2011-11-08 08:58:45 UTC
Everything is implemented in hardware.
Comment 5 Vincent Lefèvre 2011-11-09 00:48:59 UTC
The AMD Opteron doesn't have this "bug".

Now, there's a much more important problem than the rounding, on both Intel and AMD processors: the fact that their implementation is broken due to the lack of accurate range reduction. Try with X = 1.0e22 and compare the results obtained with -O0 (hardware implementation) and -O3 (obtained with MPFR, I suppose, thus correct).

In short, the hardware implementation must not be used. If the user wants a fast implementation but incorrect results on large values, he can still use compiler options like -ffast-math, with which the compiler would generate hardware instructions directly.
Comment 6 Ulrich Drepper 2011-12-22 16:33:15 UTC
All the more reason to complain to Intel.
Comment 7 gastineau 2011-12-22 16:49:07 UTC
if the user does not specify "-ffast-math", the user wants to get the same result with a call to sincos as two distincts calls to sin and cos, even if the hardware gives different results.

if you do not want correct this bug, you must specify in the documentation that the result of sincos may be not the same as two calls to sin and cos. This is a very important information for developers who take care about accuracy.
Comment 8 Ulrich Drepper 2011-12-22 17:45:01 UTC
If you don't understand how computers work then don't bother other people.  The function is implemented by a single instruction.  It's entire Intel's responsibility to provide the result.
Comment 9 gastineau 2011-12-22 19:31:39 UTC
The current glibc documentation says :

In many applications where sin and cos are used, the sine and cosine of the same angle are needed at the same time. It is more efficient to compute them simultaneously, so the library provides a 
function to do that.

— Function: void sincos (double x, double *sinx, double *cosx)
— Function: void sincosf (float x, float *sinx, float *cosx)
— Function: void sincosl (long double x, long double *sinx, long double *cosx)
These functions return the sine of x in *sinx and the cosine of x in *cos, where x is given in radians. Both values, *sinx and *cosx, are in the range of -1 to 1.

This function is a GNU extension. Portable programs should be prepared to cope with its absence.



So, if it's a GNU extension, it's not of the responsability of the user to check that it gives a different result, as the documentation says no thing about the implementation.

Mickaël




Le 22/12/11 18:45, drepper.fsp at gmail dot com a écrit :
> http://sourceware.org/bugzilla/show_bug.cgi?id=13381
>
> Ulrich Drepper<drepper.fsp at gmail dot com>  changed:
>
>             What    |Removed                     |Added
> ----------------------------------------------------------------------------
>               Status|REOPENED                    |RESOLVED
>           Resolution|                            |INVALID
>
> --- Comment #8 from Ulrich Drepper<drepper.fsp at gmail dot com>  2011-12-22 17:45:01 UTC ---
> If you don't understand how computers work then don't bother other people.  The
> function is implemented by a single instruction.  It's entire Intel's
> responsibility to provide the result.
>
Comment 10 Vincent Lefèvre 2011-12-23 11:01:30 UTC
I agree that this bug here can be seen as invalid, as long as the error is not greater than the documented error bound. It's even "normal" that sincosl doesn't return the same result as sinl and cosl, as the algorithm is different and the goal of such a hardware implementation is to provide a result quickly, though not highly accurate.

Now, there could be two valid related bugs:

1. As these hardware functions seem to be completely wrong for large arguments (on both Intel and AMD), the glibc should not use them in such cases. At least, it should provide an option so that it can be compiled to use a software implementation. Note: though relatively fast hardware range reduction is possible, I'm sure that was the intent of these hardware instructions when they were introduced in the 80's (?). I mean that blindly using these hardware instructions in a system library (i.e. without knowing the application) could be seen as a misuse of these instructions.

2. The IEEE 754-2008 standard recommends correct rounding. So, a RFE to have correct rounding would be valid, but perhaps wontfix for the extended precision (not standard, its primary goal being to ease accurate implementations of the double precision). For the double precision, unless this has been fixed, sincos is affected by this problem, which affects code using sin and cos compiled by GCC; see

  http://www.vinc17.net/research/slides/sieste2010.pdf

Slides 11-15 for this point.
Comment 11 gastineau 2011-12-23 11:20:24 UTC
The documentation does not mention that the algorithm is different and that sincosl could return a different result that cosl and sinl.
So I think that the documentation may be updated to specify that result could be differents.



Le 23/12/11 12:01, vincent-srcware at vinc17 dot net a écrit :
> http://sourceware.org/bugzilla/show_bug.cgi?id=13381
>
> --- Comment #10 from Vincent Lefèvre<vincent-srcware at vinc17 dot net>  2011-12-23 11:01:30 UTC ---
> I agree that this bug here can be seen as invalid, as long as the error is not
> greater than the documented error bound. It's even "normal" that sincosl
> doesn't return the same result as sinl and cosl, as the algorithm is different
> and the goal of such a hardware implementation is to provide a result quickly,
> though not highly accurate.
>
> Now, there could be two valid related bugs:
>
> 1. As these hardware functions seem to be completely wrong for large arguments
> (on both Intel and AMD), the glibc should not use them in such cases. At least,
> it should provide an option so that it can be compiled to use a software
> implementation. Note: though relatively fast hardware range reduction is
> possible, I'm sure that was the intent of these hardware instructions when they
> were introduced in the 80's (?). I mean that blindly using these hardware
> instructions in a system library (i.e. without knowing the application) could
> be seen as a misuse of these instructions.
>
> 2. The IEEE 754-2008 standard recommends correct rounding. So, a RFE to have
> correct rounding would be valid, but perhaps wontfix for the extended precision
> (not standard, its primary goal being to ease accurate implementations of the
> double precision). For the double precision, unless this has been fixed, sincos
> is affected by this problem, which affects code using sin and cos compiled by
> GCC; see
>
>    http://www.vinc17.net/research/slides/sieste2010.pdf
>
> Slides 11-15 for this point.
>
Comment 12 Vincent Lefèvre 2011-12-23 11:41:46 UTC
(In reply to comment #11)
> The documentation does not mention that the algorithm is different and that
> sincosl could return a different result that cosl and sinl.

It doesn't have to. The C standard doesn't even require that the following program always return true on finite values of c (i.e. x and y can be different).

int test (double c)
{
  double x, y;
  x = sin(c);
  y = sin(c);
  return x == y;
}

See my SIESTE slides to see how one can easily affect the result returned by sin().
Comment 13 martyn.j.corden 2012-01-09 20:02:37 UTC
Yes, this is likely an issue of slightly different algorithms in the math library. The main purpose of a sincos function is to give results more quickly/efficiently that two separate calls to sin and cos, and that might be better achieved with a different algorithm.
      The same issue arose in the math library provided with the Intel compiler, libimf. Following a customer report, a change was made for version 12.1 of the Intel compiler, (in Intel Composer XE 2011 SP1), so that sincosl would return the same results as calls to sinl and cosl. Presumably, similar changes could be made to the GNU libm if this is felt to be an important issue, otherwise the sincos documentation should be updated to warn of the potential difference.  If you have access to Intel’s libimf, it is also possible to link to it when compiling with gcc, to work around this issue.

Martyn Corden      (Intel Developer Support)
Comment 14 Vincent Lefèvre 2012-02-03 14:49:31 UTC
Note: for the correctness problem on large arguments (which is a much more important problem than the one reported here), I've reported bug 13658. A consequence might be that the problem reported here would no longer appear, but without any guarantee (it is not about correct rounding yet).