Bug 48140 - fmod() not accurate to double precision?
Summary: fmod() not accurate to double precision?
Status: RESOLVED INVALID
Alias: None
Product: gcc
Classification: Unclassified
Component: c (show other bugs)
Version: 4.5.2
: P3 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2011-03-15 20:35 UTC by Simon DeDeo
Modified: 2017-12-09 02:19 UTC (History)
0 users

See Also:
Host:
Target:
Build:
Known to work:
Known to fail:
Last reconfirmed:


Attachments
c code to produce the output above (260 bytes, application/octet-stream)
2011-03-15 20:35 UTC, Simon DeDeo
Details
.i file (Intel Core i7, gcc 4.5.2) (3.67 KB, application/octet-stream)
2011-03-15 20:46 UTC, Simon DeDeo
Details
.i file (Different system, gcc 4.1.2 (Gentoo 4.1.2 p1.1)) (4.70 KB, application/octet-stream)
2011-03-15 20:48 UTC, Simon DeDeo
Details

Note You need to log in before you can comment on or make changes to this bug.
Description Simon DeDeo 2011-03-15 20:35:23 UTC
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");
}
Comment 1 Simon DeDeo 2011-03-15 20:46:55 UTC
Created attachment 23671 [details]
.i file (Intel Core i7, gcc 4.5.2)
Comment 2 Simon DeDeo 2011-03-15 20:48:28 UTC
Created attachment 23672 [details]
.i file (Different system, gcc 4.1.2 (Gentoo 4.1.2 p1.1))
Comment 3 Richard Biener 2011-03-16 10:53:42 UTC
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.
Comment 4 Albert Chan 2017-12-09 02:19:28 UTC
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