CODEWITHSUNDEEP
c++floating-point
jdm
jdm

Reputation: 10150

Sorting/comparison of doubles in C(++) not stable?

I've ran into a pretty wierd problem with doubles. I have a list of floating point numbers (double) that are sorted in decreasing order. Later in my program I find however that they are not exactly sorted anymore. For example:

0.65801139819
0.6545651031    <-- a
0.65456513001   <-- b
0.64422968678

The two numbers in the middle are flipped. One might think that this problem lies in the representations of the numbers, and they are just printed out wrong. But I compare each number with the previous one using the same operator I use to sort them - there is no conversion to base 10 or similar going on:

double last_pt = 0;
for (int i = 0; i < npoints; i++) {
  if (last_pt && last_pt < track[i]->Pt()) {
    cout << "ERROR: new value " << track[i]->Pt()
         << " is higher than previous value " << last_pt << endl;
  }
  last_pt = track[i]->Pt();
}

The values are compared during sorting by

bool moreThan(const Track& a, const Track& b) {
  return a.Pt() > b.Pt();
}

and I made sure that they are always double, and not converted to float. Pt() returns a double. There are no NaNs in the list, and I don't touch the list after sorting.

Why is this, what's wrong with these numbers, and (how) can I sort the numbers so that they stay sorted?

Upvotes: 1

Views: 740

Answers (2)

James Kanze
James Kanze

Reputation: 154047

You're comparing the return value of a function. Floating point return values are returned in a floating point register, which has higher precision than a double. When comparing two such values (e.g. a.Pt() > b.Pt()), the compiler will call one of the functions, store the return value in an unnamed temporary of type double (thus rounding the results to double), then call the other function, and compare its results (still in the floating point register, and not rounded to double) with the stored value. This means that you can end up with cases where a.Pt() > b.Pt() and b.Pt() > a.Pt(), or a.Pt() > a.Pt(). Which will cause sort to get more than a little confused. (Formally, if we're talking about std::sort here, this results in undefined behavior, and I've heard of cases where it did cause a core dump.)

On the other hand, you say that Pt() "just returns a double field". If Pt() does no calculation what so ever; if it's only:

double Pt() const { return someDouble; }

, then this shouldn't be an issue (provided someDouble has type double). The extended precision can represent all possible double values exactly.

Upvotes: 1

Archie
Archie

Reputation: 6854

Are you sure you're not converting double to float at some time? Let us take a look at binary representation of these two numbers:

0 01111111110 0100111100100011001010000011110111010101101100010101
0 01111111110 0100111100100011001010010010010011111101011010001001

In double we've got 1 bit of sign, 11 bits of exponent and 53 bits of mantissa, while in float there's 1 bit of sign, 8 bit of exponent and 23 bits of mantissa. Notice that mantissa in both numbers are identical at their first 23 bits.

Depending on rounding method, there would be different behaviour. In case when bits >23 are just trimmed, these two numbers as float are identical:

0 011111110 01001111001000110010100 (trim: 00011110111010101101100010101)
0 011111110 01001111001000110010100 (trim: 10010010011111101011010001001)

Upvotes: 7

Related Questions