6

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
Daniel
  • 2,576
  • 7
  • 37
  • 51
  • 3
    What are you expecting? What are you observing? – NPE Dec 17 '12 at 16:56
  • What are you trying to do? – Mysticial Dec 17 '12 at 16:56
  • 1
    One of the problems that you are seeing is that 0.1 is hard to represent in floats or doubles http://stackoverflow.com/questions/3448777/how-to-represent-0-1-in-floating-point-arithmetic-and-decimal Your example may work better with other values of `ddd`. – DanS Dec 17 '12 at 17:21
  • @DanS, are you implying that `0.1` is the **culprit**? – Daniel Dec 17 '12 at 19:08
  • @Daniel: Like `1/3` is a repeating decimal and has no exact representation in decimal, `1/10` has no exact representation in binary -- it's `.000110011001100110011001...` – David Schwartz Dec 17 '12 at 19:37
  • @Daniel, it's not really that 0.1 is the culprit, but it exposes an issue. The problem is relying on an exact floating point value. – DanS Dec 18 '12 at 12:06

4 Answers4

7

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.

JasonD
  • 16,464
  • 2
  • 29
  • 44
1

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.

David Schwartz
  • 179,497
  • 17
  • 214
  • 278
  • I like this answer a lot actually! – Mooing Duck Dec 17 '12 at 18:10
  • Guess what, see my original post, because I cant format these results here in the comments. – Daniel Dec 17 '12 at 18:30
  • I'm not sure what I'm supposed to see there. I never said that every limited-precision representation must always give different answers for every case. Only that it is inherent in limited-precision representations that they sometimes give different answers. I showed one case where they do give different answers, you showed one where they don't. (They are different cases because my case involved decimal precision and yours involved binary precision.) – David Schwartz Dec 17 '12 at 18:41
  • (Do you agree that with six-digit decimal precision, `1/3` is represented as `.333333` and `2/3` is represented as `.6666667`? And do you not agree that if you first calculate `1/3` with six-digit decimal precision and then multiply by 2, you will get `.666666`?) – David Schwartz Dec 17 '12 at 18:43
  • I FULLY agree, but it seems the computer doesn't. – Daniel Dec 17 '12 at 19:05
  • Why do you say the computer doesn't? You never did `1/3` in six-digit decimal precision, so you haven't checked if it agrees or disagrees. You did `1/3` in `double`s. – David Schwartz Dec 17 '12 at 19:07
0

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!
}
Dims
  • 47,675
  • 117
  • 331
  • 600
  • It's not just 4.6, but from 0.1,0.2,0.3,...,4.4,4.5,4.6,4.7,.... So I have to use something divided by 0.1. – Daniel Dec 17 '12 at 18:17
  • And if 4.6 can not be represented precisely as you linked , that is 4.6 is always "4.59999999999999964472863211995E0". Then my first approach using `param = org*(1/ddd)` would be `4.59999999999999964472863211995E0 * (1/0.1) = 4.59999999999999964472863211995E0 * 10`, and then after modf, it will always be 45 and ~1, why my two approaches are getting two results? – Daniel Dec 17 '12 at 18:21
0

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
}
Zane
  • 926
  • 8
  • 21