Bug 13658 - sincos() is incorrect for large inputs on x86_64
Summary: sincos() is incorrect for large inputs on x86_64
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.13
: P2 normal
Target Milestone: ---
Assignee: Carlos O'Donell
URL:
Keywords: glibc_2.15
: 13837 (view as bug list)
Depends on:
Blocks: 13851 13852 13854 15563
  Show dependency treegraph
 
Reported: 2012-02-03 14:41 UTC by Vincent Lefèvre
Modified: 2014-06-27 09:57 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 Vincent Lefèvre 2012-02-03 14:41:59 UTC
sincos() is inaccurate for large inputs on x86_64: with glibc 2.13,

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

int main (void)
{
  volatile double x = 1.0e22;
  double s1, s2, c1;

  sincos (x, &s1, &c1);
  s2 = sin (x);
  printf ("s1 = %.17g\n", s1);
  printf ("s2 = %.17g\n", s2);
  return 0;
}

outputs:

s1 = 0.46261304076460175
s2 = -0.85220084976718879

(s2 is the correct value). I suppose that contrary to the other trig functions, glibc uses the hardware sincos instruction, which has never been meant to be used directly by a C library (the hardware elementary functions of the x86 processors were designed for small inputs, and they must not be used by code where inputs can be large, like here). The sincos() function can simply be implemented by a call to sin() and a call to cos() on this target.

Ditto for sincosf() and sincosl().

Note: x86 (32 bits) has the same problem, but it has been claimed that users don't care about correctness on this target.
Comment 1 Andreas Jaeger 2012-02-03 15:00:02 UTC
Btw. x86-64 uses the same fsincos instruction like x86.
Comment 2 Vincent Lefèvre 2012-02-03 15:09:07 UTC
For the reference about the hardware trig instructions, you can see "Intel® 64 and IA-32 Architectures Software Developer's Manual Volume 1: Basic Architecture" on:

  http://www.intel.com/content/www/us/en/processors/architectures-software-developer-manuals.html

in particular Section 8.1.3.2, which says:

"The FPTAN, FSIN, FCOS, and FSINCOS instructions set the C2 flag to 1 to indicate that the source operand is beyond the allowable range of ±2^63 and clear the C2 flag if the source operand is within the allowable range."

So, outside the interval [-2^63,+2^63] ("allowable range"), these instructions must not be used (or they can be used, but with a fallback if the C2 flag is set to 1). But note that the glibc implementation is more accurate, even with (very probably) correct rounding, so that it is better to use it anyway.
Comment 3 Carlos O'Donell 2012-02-13 04:53:41 UTC
Andreas posted a possible patch for this issue here:
http://cygwin.com/ml/libc-alpha/2012-02/msg00194.html

GCC discussion is here:
http://gcc.gnu.org/ml/gcc/2012-02/msg00072.html
Comment 4 Andreas Jaeger 2012-03-15 13:24:59 UTC
This is fixed now for sincos, sin and cos on i386 and x86-64. I'm opening separate bugs for float and long double.
Comment 5 Andreas Jaeger 2012-03-15 13:31:20 UTC
*** Bug 13837 has been marked as a duplicate of this bug. ***
Comment 6 Jackie Rosen 2014-02-16 16:56:16 UTC Comment hidden (spam)