Unwanted division operator behavior, what should I do?

Question

Problem description

During my fluid simulation, the physical time is marching as 0, 0.001, 0.002, ..., 4.598, 4.599, 4.6, 4.601, 4.602, .... Now I want to choose time = 0.1, 0.2, ..., 4.5, 4.6, ... from this time series and then do the further analysis. So I wrote the following code to judge if the fractpart hits zero.

But I am so surprised that I found the following two division methods are getting two different results, what should I do?

double param, fractpart, intpart;
double org = 4.6;
double ddd = 0.1;

// This is the correct one I need. I got intpart=46 and fractpart=0
// param = org*(1/ddd);

// This is not what I want. I got intpart=45 and fractpart=1
param = org/ddd;

fractpart = modf(param , &intpart);
Info<< "\n\nfractpart\t=\t"
    << fractpart
    << "\nAnd intpart\t=\t"
    << intpart
    << endl;

Why does it happen in this way?
And if you guys tolerate me a little bit, can I shout loudly: "Could C++ committee do something about this? Because this is confusing." :)

What is the best way to get a correct remainder to avoid the cut-off error effect? Is fmod a better solution? Thanks

Respond to the answer of

David Schwartz

double aTmp = 1;
double bTmp = 2;
double cTmp = 3;
double AAA = bTmp/cTmp;
double BBB = bTmp*(aTmp/cTmp);
Info<< "\n2/3\t=\t"
    << AAA
    << "\n2*(1/3)\t=\t"
    << BBB
    << endl;

And I got both ,

2/3     =       0.666667
2*(1/3) =       0.666667

Answer 1

Floating point values cannot exactly represent every possible number, so your numbers are being approximated. This results in different results when used in calculations.

If you need to compare floating point numbers, you should always use a small epsilon value rather than testing for equality. In your case I would round to the nearest integer (not round down), subtract that from the original value, and compare the abs() of the result against an epsilon.

If the question is, why does the sum differ, the simple answer is that they are different sums. For a longer explanation, here are the actual representations of the numbers involved:

             org:  4.5999999999999996 = 0x12666666666666 * 2^-50
             ddd: 0.10000000000000001 = 0x1999999999999a * 2^-56
           1/ddd:                  10 = 0x14000000000000 * 2^-49
   org * (1/ddd):                  46 = 0x17000000000000 * 2^-47
       org / ddd:  45.999999999999993 = 0x16ffffffffffff * 2^-47

You will see that neither input value is exactly represented in a double, each having been rounded up or down to the nearest value. org has been rounded down, because the next bit in the sequence would be 0. ddd has been rounded up, because the next bit in that sequence would be a 1.

Because of this, when mathematical operations are performed the rounding can either cancel, or accumulate, depending on the operation and how the original numbers have been rounded.

In this case, 1/0.1 happens to round neatly back to exactly 10.

Multiplying org by 10 happens to round up.

Dividing org by ddd happens to round down (I say 'happens to', but you're dividing a rounded-down number by a rounded-up number, so it's natural that the result is less).

Different inputs will round differently.

It's only a single bit of error, which can be easily ignored with even a tiny epsilon.

Answer 2

If I understand your question correctly, it's this: Why, with limited-precision arithmetic, is X/Y not the same is X * (1/Y)?

And the reason is simple: Consider, for example, using six digits of decimal precision. While this is not what doubles actually do, the concept is precisely the same.

With six decimal digits, 1/3 is .333333. But 2/3 is .666667. So:

2 / 3 = .666667  

2 * (1/3) = 2 * .333333 = .6666666  

That's just the nature of fixed-precision math. If you can't tolerate this behavior, don't use limited-precision types.