On 8/29/07, Dmitry Karasik wrote: > I'd like to submit a bug in cygwin implementation of sin(). <snip> > the difference is in 7th digit, and is significant for double precision.
This is not a bug in newlib. The problem is in glibc and msvc and newlib is (more) correct in this case. Or, to put it another way, the problem is that the 8087 was designed a few years before the proper theory for doing argument reduction for trig functions was discovered, and most subsequent x86 compatible processors (except for the AMD K5) have decided to be bug-compatible with the original 8087. If you would do your experiment on another architecture (PowerPC, IA64, x86-64 (assuming glibc built using the default -msse -mfpmath=sse compiler options) etc) you would get the same result as from newlib. The C specifications don't require exact argument reduction, and so many implementations just return whatever the underlying hardware implements. This is in opposition to, say, Java, where the FP rules are more precise. See: http://blogs.sun.com/jag/entry/transcendental_meditation Analysis of your testcase: Closest IEEE-754 double to 3.1415926535897900074 is indeed 0x1.921fb54442d11p+1 (as correctly given by gcc on both cygwin and linux). Sin(0x1.921fb54442d110000000000000000) is according to Mathematica: 0x1.d1a62633145c06e0e690...p-49 the closest IEEE-754 double to this (round-to-even) is 0x1.d1a62633145c0p-49 as given by newlib. So, in this case, newlib is getting it right to less than 0.5ulp (better than 1ulp is all that is usually possible, due to the table-maker's dilemma). The info file for glibc ("errors in math functions" node) claims that on ix86 sincos() should be accurate to 1 ulp (it doesn't give an error for sin() alone, but I reworked your test program to use sincos() and it didn't change the answer). But in this case, glibc is off by 0x2633145c0 ulps (about 10^10 ulps). Lev -- Unsubscribe info: http://cygwin.com/ml/#unsubscribe-simple Problem reports: http://cygwin.com/problems.html Documentation: http://cygwin.com/docs.html FAQ: http://cygwin.com/faq/
