[Bug c/48140] New: fmod() not accurate to double precision?

[sdedeo at post dot harvard.edu]

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=48140

           Summary: fmod() not accurate to double precision?
           Product: gcc
           Version: 4.5.2
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: c
        AssignedTo: unassig...@gcc.gnu.org
        ReportedBy: sde...@post.harvard.edu


Created attachment 23670
  --> http://gcc.gnu.org/bugzilla/attachment.cgi?id=23670
c code to produce the output above

I am not sure if this is a "good" bug to submit to gcc, since it involves the
math library. I have found it occurs on my own machine (Intel Core i7) with
gcc4.5 and gcc4.6, as well as on other machines I have access to that use gcc. 

Again, I apologize if this is either not a real error, or if it is related to
other floating point errors that people report. I have tried to look through
the database. The error below appears even without optimization flags.

In any case, I am finding that fmod, given a double precision set of inputs,
produces answers that are only good to (approximately) single precision. In
other words, the output of fmod is a double, but the answer is not accurate. It
seems to show up for large values of x in fmod(x,y).

Below I compare the outputs to Mathematica (which can handle arbitrary
precision.)

10^9 mod 2pi:      0.5773954624831035
Mathematica value: 0.5773954235013852
[this is accurate only to 10^-8]

sin(10^9 mod 2pi): 0.5458434821109812
sin(10^9):         0.5458434494486994
Mathematica value: 0.5458434494486996
[the sin function handles the large number gracefully to give standard 10^-16
precision]

#include <stdio.h>
#include <math.h>

#define TWOPI 6.2831853071795864769252867665590057683943387987502

int main() {
    double big, twopi;

    big=1000000000.0;

    printf("10^9 mod 2pi:      %#17.16g\n",  fmod(big, TWOPI));
    printf("Mathematica value: 0.5773954235013852\n");
    printf("\n");
    printf("sin(10^9 mod 2pi): %#17.16g\nsin(10^9):         %#17.16g\n",
sin(fmod(big, TWOPI)), sin(big));
    printf("Mathematica value: 0.5458434494486996\n");
}