On Sat, Feb 4, 2012 at 11:20 AM, James Courtier-Dutton
<james.dut...@gmail.com> wrote:
> On 4 February 2012 00:06, Vincent Lefevre <vincent+...@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.
