Re: 1.5.24: sin() bug

["Lev Bishop" <lev dot bishop at gmail dot com>]

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/