Bug 43490

Created attachment 20168 [details]
Trivial test program that outputs sin(x) and sinf(x) for various vals VS expected results

I haven't added tests for cos or tan etc, but I expect they will be broken as well. *My guess* is that the range reduction that gets X down to the equivalent in the range  [0, Pi/2] or [0, Pi/4] has broken as the errors get worse with increasing x.

duplicate of PR43405.

(In reply to comment #2)
> duplicate of PR43405.
> 
It's doesn't seem to be the same. I just tried the test source from  43405 and on the old system (gcc 3.4.6) (which I assumed was working) got:

         double precision: sin(1e22) = -0.8522008497671888
         quad precison: sin(1e22)=0.4626130407646018

while on the systems I was worried about, got
         double precision: sin(1e22) = 0.4626130407646017
         quad precison: sin(1e22)=0.4626130407646018

Using maple and computing the result to 30 decimal places, I get
                -.852200849767188801772705893753
so it looks like there is an additional problem.

> BUT sin(x) becomes progressively more inaccurate with increasing magnitude of
> x, as with the above version (on x86). At a guess, it would seem like something
> has broken the "range reduction" maths.  

Since pi is irrational, "range reduction" is inherently "broken" for large x. Try to compute sin(x) for a large value of x and its nearest values from below and above.

(In reply to comment #4)
> > BUT sin(x) becomes progressively more inaccurate with increasing magnitude of
> > x, as with the above version (on x86). At a guess, it would seem like something
> > has broken the "range reduction" maths.  
> 
> Since pi is irrational, "range reduction" is inherently "broken" for large x.
> Try to compute sin(x) for a large value of x and its nearest values from below
> and above.

Yes, I am aware of that, but the *standard* seems to be to assume that the source value is accurate, compute the range reduced result to sufficient precision, and then compute sine/cosine to the required ULP accuracy. (e.g. see Payne-Hanek reduction algorithm or any of the newer methods such as ftp://ftp.inria.fr/INRIA/publication/dienst/RR-4267.pdf).

It was working on earlier incarnations so is frustrating that it's now broken.

(In reply to comment #3)
> (In reply to comment #2)
> > duplicate of PR43405.

> Using maple and computing the result to 30 decimal places, I get
>                 -.852200849767188801772705893753
> so it looks like there is an additional problem.

you can test two different implementations:
1). '-O2 -m32 -fno-builtin' - force libm.so calls and test libc implementation.
2). '-O2 -m32' to test gcc compile-time evaluation.

This is a glibc problem.  Confirmed by LD_PRELOADing a different libm which
makes it work.

Please report to glibc bugzilla instead.

(In reply to comment #6)

> you can test two different implementations:
> 1). '-O2 -m32 -fno-builtin' - force libm.so calls and test libc implementation.
> 2). '-O2 -m32' to test gcc compile-time evaluation.
> 

FWIW I've just tried
gcc -o simplesintest -Wall -W -pedantic -std=c99  simplesintest.c -O2 -m32 -fno-builtin -lm
  and
gcc -o simplesintest -Wall -W -pedantic -std=c99  simplesintest.c -O2 -m32 -lm

and they've made no difference :-|

(In reply to comment #7)
> This is a glibc problem.  Confirmed by LD_PRELOADing a different libm which
> makes it work.
> 
> Please report to glibc bugzilla instead.
> 

The glibc team are saying that nothing has changed in the maths library (http://sourceware.org/bugzilla/show_bug.cgi?id=11422#c5) which brings me back here as it would appear that the only difference between the working and broken libraries are the versions of gcc used to compile them, e.g. 

working library: 
   GNU C Library stable release version 2.3.4, by Roland McGrath et al.
   Compiled by GNU CC version 3.4.6 20060404 (Red Hat 3.4.6-11).

example broken library:
    GNU C Library (EGLIBC) stable release version 2.10.1, 
    gcc version 4.4.1 (Ubuntu 4.4.1-4ubuntu9)

Both sin and sinf are assembly routines so can't probably be affected by
the gcc used to compile them.  See sysdeps/i386/fpu/s_sinf.S and
sysdeps/i386/fpu/s_sin.S.

I have not seen a working library for your testcase (2.9, built with
GCC 4.3.2 is broken).

This isn't a GCC bug but a glibc bug.  Clearly glibc uses fsin out of specs.

(In reply to comment #10)
> Both sin and sinf are assembly routines so can't probably be affected by
> the gcc used to compile them.  See sysdeps/i386/fpu/s_sinf.S and
> sysdeps/i386/fpu/s_sin.S.
> 
> I have not seen a working library for your testcase (2.9, built with
> GCC 4.3.2 is broken).

Hmm. The last lot of code I saw was all generic c. Interesting

> 
> This isn't a GCC bug but a glibc bug.  Clearly glibc uses fsin out of specs.

I can't see that *is* this issue unless fsin is also fundamentally broken. The Intel "Instruction Set Reference" states fsin works in the range -2^63 to +2^63 and I have just tested a value ~2^7 and got the following

sin(1240093312.000000) (i.e. sin(0x93d4a5 * 2^7))
        is computed as  -4.581222e-08 (-0x1.898622f829ffep-25)
        It SHOULD be ~  -4.581062e-08 (-0x1.8982a036335f2p-25)
        Relative error is 3.485e-03%

That is dreadfully inaccurate and corresponds to about 15 bits accuracy instead of an expected result of about 52~53 bits for a double.

Is there ANY way to get the gcc and glibc developers to talk directly to each other rather than have me ping-ponging back and forth between the two bugzilla databases?

(In reply to comment #11)
> (In reply to comment #10)

> and I have just tested a value ~2^7 and got the following
Sorry I meant to say ~2^30

Same error with glibc 2.5.

Note that fsin is not part of IEEE and thus does not guarantee 1ulp precision
(in fact it is known for its inaccuracies).

Please do not keep re-opening this bug.  It is not a GCC bug but a GLIBC
bug.  I don't think it is a GLIBC regression either, unless the x87 assembly
wasn't there in 2.3.4.

i've checked the Test0 from example a little on my intel Q9300 cpu.

first case: setup 32-bit code and break on libm.sin implementation.

Breakpoint 1, main (argc=1, argv=0xffffe3c4) at t.c:74
(gdb) b sin
Breakpoint 2 at 0xf7fa7110: file ../sysdeps/i386/fpu/s_sin.S, line 14.

Test 0:
(gdb) p testvalue
$1 = 43998769152

Breakpoint 2, sin () at ../sysdeps/i386/fpu/s_sin.S:14
(gdb) bt
#0  sin () at ../sysdeps/i386/fpu/s_sin.S:14
#1  0x08048682 in main (argc=1, argv=0xffffe3c4) at t.c:81

(gdb) info float
=>R7: Valid   0x4022a3e87f0000000000 +43998769152

11x RCSID("$NetBSD: s_sin.S,v 1.5 1995/05/09 00:25:54 jtc Exp $")
12x
13x ENTRY(__sin)
14x         fldl    4(%esp)
15x         fxam
16x         fstsw   %ax
17x         movb    $0x45, %dh
18x         andb    %ah, %dh
19x         cmpb    $0x05, %dh
20x         je      3f
21x 4:      fsin                  <==== here we calculate Test0.
22t>        fnstsw  %ax
23x         testl   $0x400,%eax
24x         jnz     1f
25x         ret                   <==== and return...
26x         .align ALIGNARG(4)
27x 1:      fldpi
28x         fadd    %st(0)
29x         fxch    %st(1)
30x 2:      fprem1
31x         fnstsw  %ax
32x         testl   $0x400,%eax
33x         jnz     2b
34x         fstp    %st(1)
35x         fsin
36x         ret
37x 3:

and the hardware result is:

(gdb) info float
=>R7: Valid   0xbfe38c41185b67ffffe4 -4.081936725100227664e-09

sin(43998769152.000000) (i.e. sin(0xa3e87f * 2^12))
        is computed as  -4.081937e-09 (-0x1.188230b6dp-28)
        It SHOULD be ~  -4.025292e-09 (-0x1.149dafd6b8987p-28)
        Relative error is 1.407% !!!!

as you can see for 32-bit implementation there's *NO* smart
argument range reduction.

second case: setup 64-bit code.

(gdb) break sin
Breakpoint 2 at 0x7ffff7b79e00: file ../sysdeps/ieee754/dbl-64/s_sin.c, line 90.
(gdb) c

Test 0:
Continuing.
Breakpoint 2, __sin (x=43998769152) at ../sysdeps/ieee754/dbl-64/s_sin.c:90

and here's some kind of range reduction for -105414350 <|x|< 281474976710656.

(gdb) fin
Run till exit from #0  __sin (x=43998769152) at ../sysdeps/ieee754/dbl-64/s_sin.c:235
0x0000000000400882 in main (argc=1, argv=0x7fffffffe298) at t.c:81
Value returned is $1 = -4.0252920638371055e-09

sin(43998769152.000000) (i.e. sin(0xa3e87f * 2^12))
        is computed as  -4.025292e-09 (-0x1.149dafd6b8987p-28)
        It SHOULD be ~  -4.025292e-09 (-0x1.149dafd6b8987p-28)
        Relative error is 0.000e+00%


finally, as you can see this NOT a gcc bug. it's a QOI bug in glibc ;)

Thanks. I'll go back to the glibc bugzilla.