Created attachment 23670 [details] 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"); }
Created attachment 23671 [details] .i file (Intel Core i7, gcc 4.5.2)
Created attachment 23672 [details] .i file (Different system, gcc 4.1.2 (Gentoo 4.1.2 p1.1))
This is at most a bug in the C library. Note that CPUs do not have native support for fmod or remainder so the implementations need to weight speed against precision.
I think fmod is accurate Gcc cannot do exact fmod(1e9, 2*pi) it can only do fmod(1e9, 0x1.921f54442d18p+2) = 0.57739546248310347 which is exactly what gcc returns in fact, it is VERY accurate If we do remainder calculation in steps: quotient = floor(1e9 / (2 pi)) = 159154943 remainder = 1e9 - quotient * 0x1.921f54442d18p+2 = 0.57739543914794922 doing in steps under-counted 210,184,384 ULP