Reputation: 10150
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
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
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