Bug 11422

Summary: sin(x) (actually probably all trig) is inaccurate for large x
Product: glibc Reporter: Simon Fenney <simon.fenney>
Component: libcAssignee: Ulrich Drepper <drepper.fsp>
Status: RESOLVED INVALID    
Severity: normal CC: glibc-bugs, vr5
Priority: P2 Flags: fweimer: security-
Version: 2.10   
Target Milestone: ---   
Host: i486-linux-gnu Target: i486-linux-gnu
Build: Last reconfirmed:
Attachments: simple test program to illustrate some error cases (gzipped)

Description Simon Fenney 2010-03-23 14:13:46 UTC
Was reported to (http://gcc.gnu.org/bugzilla/show_bug.cgi?id=43490) but was then
sent here...

sin(x) (and sinf(x)) *used* to work accurately on earlier gcc versions, e.g.
    gcc 3.4.6 20060404 (Red Hat 3.4.6-9) ON x86_64-redhat-linux
    gcc 4.2.1  on Macintosh
OR
    cross compiling for ARM Cortex-A8 with gcc version 4.3.2

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. 

Systems that were found to be broken:
    gcc 4.4.1 (Ubuntu 4.4.1-4ubuntu9)
    gcc 4.4.3 20100108 (prerelease) (Debian 4.4.2-9)
    GNU C 4.3.2 (i686-pc-linux-gnu) (Gentoo 4.3.2-r3 p1.6, pie-10.1.5


Build command: gcc -v -save-temps -o simplesintest  -g -DDEBUG=1 -DDEBUG_INFO=1
-Wall -W -pedantic -std=c99  simplesintest.c -lm

Command line: ./simplesintest

Typical output is:
=========================
Test 0:
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% !!!!

sinf(43998769152.000000) (i.e. sinf(0xa3e87f * 2^12))
        is computed as  -4.081937e-09 (-0x1.18823p-28)
        It SHOULD be ~  -4.025292e-09 (-0x1.149dbp-28)
        Relative error is 1.407% !!!!


Test 1:
sin(9903547467890318699652972544.000000) (i.e. sin(0x800017 * 2^70))
        is computed as  2.763719e-01 (0x1.1b013a2290cc8p-2)
        It SHOULD be ~  2.003433e-01 (0x1.9a4d9074cecaap-3)
        Relative error is -37.949% !!!!

sinf(9903547467890318699652972544.000000) (i.e. sinf(0x800017 * 2^70))
        is computed as  2.763719e-01 (0x1.1b013ap-2)
        It SHOULD be ~  2.003433e-01 (0x1.9a4d9p-3)
        Relative error is -37.949% !!!!


Test 2:
sin(8773115793420245409943394342404096.000000) (i.e. sin(0xd84625 * 2^89))
        is computed as  7.399661e-01 (0x1.7adcd596e1d56p-1)
        It SHOULD be ~  2.044974e-08 (0x1.5f52ea84f120cp-26)
        Relative error is -3618461696.000% !!!!

sinf(8773115793420245409943394342404096.000000) (i.e. sinf(0xd84625 * 2^89))
        is computed as  7.399661e-01 (0x1.7adcd6p-1)
        It SHOULD be ~  2.044974e-08 (0x1.5f52eap-26)
        Relative error is -3618461696.000% !!!!


Test 3:
sin(10633823966279326983230456482242756608.000000) (i.e. sin(0x800000 * 2^100))
        is computed as  -7.480374e-01 (-0x1.7efec078c6a44p-1)
        It SHOULD be ~  -4.205416e-02 (-0x1.5881f6eeb6a8dp-5)
        Relative error is 1678.748% !!!!

sinf(10633823966279326983230456482242756608.000000) (i.e. sinf(0x800000 *
2^100))
        is computed as  -7.480373e-01 (-0x1.7efecp-1)
        It SHOULD be ~  -4.205416e-02 (-0x1.5881f6p-5)
        Relative error is 1678.748% !!!!
=========================
Comment 1 Simon Fenney 2010-03-23 14:16:59 UTC
Created attachment 4673 [details]
simple test program to illustrate some error cases  (gzipped)
Comment 2 Simon Fenney 2010-03-23 14:20:07 UTC
(In reply to comment #1)
> Created an attachment (id=4673)
> simple test program to illustrate some error cases  (gzipped) 
> 

Forgot to mention that one of the older systems that works has libc version  2.3.4
Comment 3 Ulrich Drepper 2010-03-23 19:38:31 UTC
The x86 code forever used the FPU instructions.  If you think they changed then
complain to the CPU manufacturers.
Comment 4 Simon Fenney 2010-03-24 10:46:02 UTC
With all due respect, it can't just be changes in the CPU. I just ran code
compiled with  
  gcc version 3.4.6 20060404 (Red Hat 3.4.6-11)
with libc version 2.3.4
i.e.
   GNU C Library stable release version 2.3.4, by Roland McGrath et al.
   Copyright (C) 2005 Free Software Foundation, Inc.
   This is free software; see the source for copying conditions.
 
   Compiled by GNU CC version 3.4.6 20060404 (Red Hat 3.4.6-11).
   Compiled on a Linux 2.4.20 system on 2009-06-01.

On a 4 month old processor:
       vendor_id       : GenuineIntel
       cpu family      : 6
       model           : 26
       model name      : Intel(R) Xeon(R) CPU           L5520  @ 2.27GHz
       stepping        : 5
       cpu MHz         : 2268.000

and *THAT* works. (I have also tried this build with older and non-Intel CPUs)


The bug, OTOH, occurs with (amongst others) 
    gcc version 4.4.1 (Ubuntu 4.4.1-4ubuntu9)

    GNU C Library (EGLIBC) stable release version 2.10.1, 
    Compiled by GNU CC version 4.4.1.
    Compiled on a Linux >>2.6.24-23-server<< system on 2010-01-03.

with both new CPUs (e.g. Atom) and those that are several years old, eg
   "Intel(R) Pentium(R) 4 CPU 2.80GHz"

Something else *must* have changed, either in the source code for the maths
library OR in the way gcc compiles it. The gcc guys are claiming it's not their
responsibility (http://gcc.gnu.org/bugzilla/show_bug.cgi?id=43490#c7) but I
don't think it's fair to lay the blame on the CPU manufacturers.
Comment 5 Ulrich Drepper 2010-03-24 11:15:45 UTC
I wrote nothing changed in the code and that exclusively the FPU instructions
are used on x86.  Whatever other problems you have are your own doing.  This is
no place for you to get programming advise.  If you don't know about the way the
x86 FPU works with it's temporary 80 bit values I suggest you read up or you
stop relying on that level of precision on x86.
Comment 6 Simon Fenney 2010-03-24 12:13:13 UTC
Please don't get annoyed, but I am well aware of the joy of x86 80-bit float HW.

All I'm trying to establish is why the standard maths library (and by the look
of it is the range reduction (a la Payne-Haneck) which for all intents and
purposes are integer operations though often implemented using doubles)) is no
longer functioning correctly. 

Since you are saying that this source code has not changed, then I should be
able to go back to the gcc team and state that the problem is at their end and
stop this buck-passing.
Comment 7 Simon Fenney 2010-03-24 14:42:18 UTC
Sigh... the gcc guys have come back with 
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=43490#c13  etc





Comment 8 Ulrich Drepper 2010-03-24 22:54:09 UTC
Stop reopening the bug.  I wrote many times that nothing changed.  Whatever you
imagined to have worked in the past is the same today if the processor behaves
the same.

Of course the FPU instructions aren't accurate.  But this is nothing new and
will never change for x86.  Use x86-64 if you want more accuracy.
Comment 9 Simon Fenney 2010-03-25 11:48:31 UTC
Dear Ulrich, 
For my own sanity, I created a vintage SuSE 8 system and now see that sin was
indeed broken back then. I apologise for the assertion that it had been working
on older systems.

However, I don't believe x86-64 is an option as the code has to run on Atom
platforms. 

The thing that is effectively broken is the "range reduction" (i.e. mapping X to
the equivalent value between 0 and Pi/2 (or perhaps Pi/4) which is, in effect, a
(multi-word precision) integer operation. Can't the 32-bit builds use the same
code base that the x86-64 systems use (which presumably avoid using the dubious
fsin instruction on non trivial values)?

Regards
Simon
Comment 10 Ulrich Drepper 2010-03-26 14:27:17 UTC
(In reply to comment #9)
> However, I don't believe x86-64 is an option as the code has to run on Atom
> platforms. 

There are 64-bit Atoms.


> Can't the 32-bit builds use the same
> code base that the x86-64 systems use (which presumably avoid using the dubious
> fsin instruction on non trivial values)?

No.  There are people asking for the opposite: using the inaccurate x86 code on
x86-64.  The code is what it is.  The good neas is that it's trivial to use your
own libm replacement once you write it.