Re: weird optimization in sin+cos, x86 backend

[Richard Guenther <richard dot guenther at gmail dot com>]

On Sat, Feb 4, 2012 at 11:20 AM, James Courtier-Dutton
<james.dutton@gmail.com> wrote:
> On 4 February 2012 00:06, Vincent Lefevre <vincent+gcc@vinc17.org> wrote:
>> On 2012-02-03 17:40:05 +0100, Dominique Dhumieres wrote:
>>> While I fail to see how the "correct value" of
>>> cos(4.47460300787e+182)+sin(4.47460300787e+182)
>>> can be defined in the 'double' world, cos^2(x)+sin^2(x)=1 and
>>> sin(2*x)=2*sin(x)*cos(x) seems to be verified (at least for this value)
>>> even if the actual values of sin and cos depends on the optimisation level.
>>
>> Actually this depends on the context. It is even worse: the value
>> of sin(some_value) depends on the context. Consider the following
>> program:
>>
>> #include <stdio.h>
>> #include <math.h>
>>
>> int main (void)
>> {
>> ?double x, c, s;
>> ?volatile double v;
>>
>> ?x = 1.0e22;
>> ?s = sin (x);
>> ?printf ("sin(%.17g) = %.17g\n", x, s);
>>
>> ?v = x;
>> ?x = v;
>> ?c = cos (x);
>> ?s = sin (x);
>> ?printf ("sin(%.17g) = %.17g\n", x, s);
>>
>> ?return c == 0;
>> }
>>
>> With "gcc -O" on x86_64 with glibc 2.13, one gets:
>>
>> sin(1e+22) = -0.85220084976718879
>> sin(1e+22) = 0.46261304076460175
>>
>> In the program, you can replace 1.0e22 by 4.47460300787e+182 or
>> whatever large constant you want. You'll get the same problem.
>>
>
> Ok, so the value -0.85... seems to be the correct value.
>
> If you convert all the doubles into long doubles.
> You get:
> gcc -O0 -c -o sincos.o sincos.c
> gcc sincos.o -lm
> sinl(10000000000000000000000.00000000000000000) = 0.46261304076460176
> sinl(10000000000000000000000.00000000000000000) = 0.46261304076460176
> gcc sincos.o -lm
> gcc -O2 -c -o sincos.o sincos.c
> sin(10000000000000000000000.00000000000000000) = -0.85220084976718880
> sin(10000000000000000000000.00000000000000000) = 0.46261304076460176
>
> So, I agree with Vincent Lefevre, my earlier statements are wrong.
> There is definitely something going wrong either in gcc or libm.

Note that you are comparing a constant folded sin() result against sincos()
(or sin() and cos()).  Use

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

int main (void)
{
  double x, c, s;
  volatile double v;

  x = 1.0e22;
  v = x;
  x = v;
  s = sin (x);
  printf ("sin(%.17g) = %.17g\n", x, s);

  v = x;
  x = v;
  c = cos (x);
  s = sin (x);
  printf ("sin(%.17g) = %.17g\n", x, s);

  return c == 0;
}

to compare libm sin against sincos.  Works for me, btw.  Note that I remember
that at some point sincos() was just fsincos even on x86_64 (ugh).

Richard.