fmod is inconsistent with division

Question

usually we got same result(not exactly same, But nearly same) for following code (suggest X > 0 , Y > 0):

code 1 : fmod(X,Y)
code 2 : X - Y * trunc(X/Y)

But for some value, they will be very different.

test code:

int main()
{
    cout << std::boolalpha;
    cout << std::setprecision(20);

    double X = 30.508474576271183309;
    double Y = 6.1016949152542370172;

    cout << "X: " << X << " Y:" << Y << '\n';
    cout << "X/Y: " << X / Y << '\n';
    cout << "X - trunc(X/Y)*Y: " << X - trunc(X / Y)*Y << '\n';
    cout << "fmod(X,Y): " << fmod(X, Y) << '\n';


    system("pause");
    return 0;
}

result:

X: 30.508474576271183309 Y:6.1016949152542370172
X/Y: 5
X - trunc(X/Y)*Y: 0
fmod(X,Y): 6.1016949152542352408

env: VS2017 Win10

UPDATE 1 : this problem is not same as c++ fmod returns 0.2 for fmod(0.6,0.2)

if you set X to 0.6 and Y to 0.2, you will get same result for code 1 and code 2. But in this situation they will be different.

UPDATE 2 : I will use another way to describe the problem.

double X, X > 0, known
double Y, Y > 0, known
int quotient, quotient > 0, unknown
double remainder, remainder > 0, unknown

formula 1: X == Y * quotient + remainder

for known X,Y , if I want calculate quotient and remainder (and make sure they satisfied with the formula 1), in c++, usually we have two ways to do that:

Method 1: quotient = trunc(X/Y), remainder = fmod(X,Y);
Method 2: quotient = trunc(X/Y), remainder = X - Y * trunc(X/Y);

for Method 1, it will output quotient and remainder but they will not satisfied with the formula 1 (for specific X and Y).

UPDATE 3 : I meet this problem when doing ADS-B CPR local decoding. There is a step called "Calculate the latitude index j" as followed

if we strictly follow it and just replace mod with fmod in C++, then we will get wrong result for specific Lat_ref and dLat.

https://stackoverflow.com/questions/38447548/c-fmod-returns-0-2-for-fmod0-6-0-2?rq=1 — James Picone, Aug 30 '17 at 07:50
The "make a string to `double`" is highly likely to give UB, since the alignment will be defined by the `char` type, not by alignment of `double`. Not really related to the results, just that it's not a particularly great construction. — Mats Petersson, Aug 30 '17 at 07:56
[Just write the number!](http://en.cppreference.com/w/cpp/language/floating_literal) `double X = 30.5084745762712; double Y = 6.10169491525424;` (with as many digits as you care to be precise to) — Caleth, Aug 30 '17 at 11:03
If I set X to 30.5084745762712, Y to 6.10169491525424. Then Code1 will output 3.5527136788005e-15, Code2 will output 2.66453525910038e-15. This is just a normal result. the difference between Code1 and Code2 is very small in this situation. @Caleth — Barrypp.zzx, Aug 30 '17 at 14:16
Is there any good way to represent a known double var without loss of precision ? — Barrypp.zzx, Aug 30 '17 at 14:25
Use more decimal digits. ~20 should suffice ~60 is way over. See also: [`std::numeric_limits::digits10`](http://en.cppreference.com/w/cpp/types/numeric_limits/digits10) — Caleth, Aug 30 '17 at 14:26
The thing about checking `double Z = fmod(X, Y)` is that you have to test **both** `abs(Z) < eps` **and** `abs(Z - Y) < eps` to see if the result is ~= 0 — Caleth, Aug 30 '17 at 14:34
I mean "without loss of precision". not losing any bit of precision. what is "eps" ? Why should I do that check ? @Caleth — Barrypp.zzx, Aug 30 '17 at 15:33
`eps` is some value below which you regard your doubles as equivalent. Because using bitwise equality,`double` arithmetic **doesn't** follow the mathematical laws of real arithmetic, e.g. addition isn't commutative. `a + (b + c)` can differ from `(a + b) + c` — Caleth, Aug 30 '17 at 16:24
The maximum that those sums can differ by is `max(a, b, c) * std::numeric_limits::epsilon` — Caleth, Aug 30 '17 at 16:25
You can use the "p" format, which uses hex-digits to precisely set every bit in a floating point value. — Mats Petersson, Aug 31 '17 at 06:56
If you do: `double val = X; for(int i = 0; i < 5; i++) { val -= Y; printf("val = %a\n", val); }` You will see that you are left with a "tiny value". Since the implementation of `fmod` a bunch of bit manipulation, not floating point calculations, there is no rounding of the small delta, and the result is a bigger error than you'd expect. — Mats Petersson, Aug 31 '17 at 07:48
@MatsPetersson So, the result of fmod is more accurate ? Could you tell me a better way to get quotient and remainder ? — Barrypp.zzx, Aug 31 '17 at 12:08
Neither is "more accurate", just "different". The key point here is that floating point math is "fiddly" and some values will have surprising results like this. If you use division `(X - trunc(X / Y)*Y)`, it will do some rounding of the result. The way (typically) that `fmod` is implemented is by bitfiddling. When I compile the above code with -m32 (so it uses the x87 FPU's fprem instruction), then the result is different. I have not studied the exact bit fiddling algorithm, but it clearly "sees" the small difference as bigger than it "should" - in this particular case. — Mats Petersson, Sep 01 '17 at 06:10
If you generally get the correct result with `(X - trunc(X / Y)*Y)`, then I would use that, unless performance difference with `fmod` is more important than the right result. Of course, if you "don't know what is the right result", I'm not sure if you can determine that one is better than the other... — Mats Petersson, Sep 01 '17 at 06:18

fmod is inconsistent with division

0 Answers0