Erroneous results in trigonometric functions for > double-precision values

[daniel.is.fischer]

For large inputs, sin, cos and tan return the input value:

GHCi, version 6.8.1: http://www.haskell.org/ghc/  :? for help
Loading package base ... linking ... done.
Prelude> sin 1e20
1.0e20
Prelude> cos 1e20
1.0e20
Prelude> cos 1e19
1.0e19
Prelude> cos 1e18
0.11965025504785125
Prelude> tan 1e19
1.0e19

[ajd]

The problem is that the Double type cannot store that much precision and so GHC can't calculate trig functions at all accurately. (This is actually documented, see http://www.haskell.org/ghc/docs/latest/html/users_guide/wrong-compilee.html. Since this is the way the compiler is documented to work, I'm changing the status to feature request.)

From my understanding, in order to solve this, one would need to implement an arbitrary precision floating-point system in GHC (currently, we only have arbitrary precision integers through GMP). This could be done through something like MPFR, but I'm guessing that would be a huge job and potentially plagued with the same problems as GMP, as the latter is a dependency of the former.

[Frederik]

I thought that the value of 'cos' and 'sin' should be between -1 and 1. Is 1e20 really the closest we can afford to get to this interval? gcc does much better:

$ cat t.c
#include <math.h>
#include <stdio.h>

int main() {
  printf("%lf\n",cos(1e20));
  return 0;
}
$ gcc t.c -lm
$ ./a.out    
-0.664912

This answer may not be very meaningful, but at least it is in the proper range.

[daniel.is.fischer]

I reported it because the behaviour was different up to 6.6.1, didn't expect it to be an intentional change (not that I'd rely on even the sign of the result for such large inputs).

[Frederik]

I see that GHC 6.4.2 gives the same results as C

Prelude> cos(1e20)
-0.6649117899070088

Ajd, I hope it is OK to change this back to a bug, it certainly seems more than a feature request...

[igloo]

I can reproduce this on x86/Linux, but not amd64/Linux.

[igloo]

OK, the problem is that the fsin instruction (and, I assume, its friends) only take arguments with absolute value <= 2^63. The reason it works in C is that they look to see if the out-of-range flag is set when they call it, and if so take the sine of the remainder after dividing by 2 * pi (only it's slightly more complicated, as they really take the partial remainder, which means they might have to do it multiple times):

0xf7f810f0 <sin+0>:     fldl   0x4(%esp)
0xf7f810f4 <sin+4>:     fsin
0xf7f810f6 <sin+6>:     fnstsw %ax
0xf7f810f8 <sin+8>:     test   $0x400,%eax
0xf7f810fd <sin+13>:    jne    0xf7f81100 <sin+16>
0xf7f810ff <sin+15>:    ret    
0xf7f81100 <sin+16>:    fldpi  
0xf7f81102 <sin+18>:    fadd   %st(0),%st
0xf7f81104 <sin+20>:    fxch   %st(1)
0xf7f81106 <sin+22>:    fprem1 
0xf7f81108 <sin+24>:    fnstsw %ax
0xf7f8110a <sin+26>:    test   $0x400,%eax
0xf7f8110f <sin+31>:    jne    0xf7f81106 <sin+22>
0xf7f81111 <sin+33>:    fstp   %st(1)
0xf7f81113 <sin+35>:    fsin   
0xf7f81115 <sin+37>:    ret    

I guess we should do the same.

[igloo]

Fixed in HEAD and 6.8:

Fri May  2 17:08:41 PDT 2008  Ian Lynagh <igloo@earth.li>
  * Fix sin/cos/tan on x86; trac #2059
  If the value is > 2^63 then we need to work out its value mod 2pi,
  and apply the operation to that instead.