sinl (and probably others) are not computed correctly. At least for large inputs. Please consider the following simple testcase: $ cat sintest.c #include <stdio.h> #include <math.h> int main (void) { double arg = 1e22; long double larg = 1e22L; printf("double precision: sin(1e22) = %.16lf\n", sin(arg)); printf("quad precison: sin(1e22)=%.16Lf\n", sinl(larg)); return 0; } $gcc sintest.c $./a.out double precision: sin(1e22) = -0.8522008497671888 quad precison: sin(1e22)=0.4626130407646018 This behavior is demonstrated by gcc-4.4.3 and gcc-4.5 (current snapshot) I am using ubuntu on a 64 bit linux machine (intel). However, the same results were obtained on different machines.
Created attachment 20131 [details] testcase as a standalone file
this is a bug in glibc-2.11.1/sysdeps/x86_64/fpu/s_sinl.S
Glibc is a separate project, see http://sources.redhat.com/bugzilla/
more details... intel (24319101.pdf) manual describe requirements for fsin opcode: "Calculates the sine of the source operand in register ST(0) and stores the result in ST(0). The source operand must be given in radians and must be within the range −2^63 to +2^63." the 1e+22 is greater than allowed 2^63 ~ 9.22e+18. so this bug is rejected by hardware ;)
The very same code compiled by the Intel C compiler runs as expected. Moreover, the prototype of sinl is as follows long double sinl(long double x); and 1e22 definitely withing the bounds of long double. (In reply to comment #4) > more details... > > intel (24319101.pdf) manual describe requirements for fsin opcode: > > "Calculates the sine of the source operand in register ST(0) and stores > the result in ST(0). The source operand must be given in radians and must > be within the range −2^63 to +2^63." > > the 1e+22 is greater than allowed 2^63 ~ 9.22e+18. > > so this bug is rejected by hardware ;) > (In reply to comment #4) > more details... > > intel (24319101.pdf) manual describe requirements for fsin opcode: > > "Calculates the sine of the source operand in register ST(0) and stores > the result in ST(0). The source operand must be given in radians and must > be within the range −2^63 to +2^63." > > the 1e+22 is greater than allowed 2^63 ~ 9.22e+18. > > so this bug is rejected by hardware ;) > (In reply to comment #4) > more details... > > intel (24319101.pdf) manual describe requirements for fsin opcode: > > "Calculates the sine of the source operand in register ST(0) and stores > the result in ST(0). The source operand must be given in radians and must > be within the range −2^63 to +2^63." > > the 1e+22 is greater than allowed 2^63 ~ 9.22e+18. > > so this bug is rejected by hardware ;) >
Also: 1e22 is not exactly representable as a floating point number. By consequence, 1e22 is different numbers when stored as a double or a long double, and we should expect different results when applying the sine to it. W.
Subject: Re: sinl is not computed correctly On Thu, 18 Mar 2010, bangerth at gmail dot com wrote: > Also: 1e22 is not exactly representable as a floating point number. By 5**22 is smaller than 2**53, so 1e22 (= 5**22 * 2**22) *is* exactly representable in double or formats with more bits than double.
(In reply to comment #6) > Also: 1e22 is not exactly representable as a floating point number. By > consequence, 1e22 is different numbers when stored as a double or a > long double, and we should expect different results when applying the > sine to it. > > W. > This is *not* a problem of binary representation. I tried to use exactly the same argument (one time double, another time *promoted* to long double). The function is not implemented correctly, or, at least for some numbers in legal range.
(In reply to comment #7) > Subject: Re: sinl is not computed correctly > > On Thu, 18 Mar 2010, bangerth at gmail dot com wrote: > > > Also: 1e22 is not exactly representable as a floating point number. By > > 5**22 is smaller than 2**53, so 1e22 (= 5**22 * 2**22) *is* exactly > representable in double or formats with more bits than double. > 1e22 is not representable exactly in the binary format. Not because of its magnitude, just because most numbers cannot be represented exactly in computers.