CODEWITHSUNDEEP
c++
Daniel
Daniel

Reputation: 2726

Unwanted division operator behavior, what should I do?

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

Upvotes: 6

Views: 254

Answers (4)

JasonD
JasonD

Reputation: 16612

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.

Upvotes: 7

Dims
Dims

Reputation: 51229

You can't represent 4.6 precisely: http://www.binaryconvert.com/result_double.html?decimal=052046054

Use rounding before separating integer and fraction parts.

UPDATE

You may wish to use rational class from Boost library: http://www.boost.org/doc/libs/1_52_0/libs/rational/rational.html

CONCERNING YOUR TASK

To find required double take precision into account, for example, to find 4.6 calculate "closeness" to it:

double time;

...

double epsilon = 0.001;

if( abs(time-4.6) <= epsilon ) {
    // found!
}

Upvotes: 0

Zane
Zane

Reputation: 926

Hm not really sure what you want to achieve, but if you want get a value and then want to do some refine in the range of 1/1000, why not use integers instead of floats/doubles?

You would have a divisor, which is 1000, and have values that you iterate over that you need to multiply by your divisor.

So you would get something like

double org = ... // comes from somewhere
int divisor = 1000;
int referenceValue = org * div;
for (size_t step = referenceValue - 10; step < referenceValue + 10; ++step) {
   // use   (double) step / divisor to feed to your algorithm
}

Upvotes: 0

David Schwartz
David Schwartz

Reputation: 182837

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.

Upvotes: 1

Related Questions