3

Short question: is the behavior of comparison operators defined at max precision? Meaning, if I have two numbers (x and y) that are identical up to precision, should I always expect the same answer when I do x < y?

I realize this might seem like a dumb question, but let me elaborate. I am working with double and I have an array of numbers like: 

0:  62536.5477752959
1:  62536.4840613718
2:  62536.4576412381
3:  62522.8487197062
4:  62536.5473896233
5:  62536.5467941254
6:  62527.3508907998
7:  62536.5477752959
8:  62517.5900098039
9:  62536.5477752959

Notice that entries 0, 7 and 9 have the same value.

When I do (something like this) : 

int low = 0, high = 0;
for(int i = 0; i < N; ++i) {
  if(x[i] < x[low])
    low = i;
  if(x[i] > x[high]) 
    high = i;
 }
 cout << "low: " << low << " high: " << high << endl;

I sometimes get: low: 8 high: 0, sometimes: low: 8 high: 7 I would have expected always the lowest index value.

Any ideas?

[edit missing braces.]

Randi
  • 320
  • 2
  • 11
  • How did you load the values into the array? – wally Sep 18 '17 at 15:51
  • 1
    Maybe replate/dupe: https://stackoverflow.com/questions/588004/is-floating-point-math-broken – NathanOliver Sep 18 '17 at 15:52
  • 2
    There's no such thing like a _specified precision_ for `double` values? Are you confusing that with their text representeation? – user0042 Sep 18 '17 at 15:53
  • There might be some braces missing in the code provided. – wally Sep 18 '17 at 15:58
  • `i` is never initialized, so that could be undefined behavior and give inconsistent results. – wally Sep 18 '17 at 16:01
  • The longer story is that this is happening in one of my optimization classes. The class has a a function pointer to an user provided cost-function that evaluates for example y = x^2. The class provides x and gets back y, which is the values of the arrays I am showing. @rex, thanks for the braces and initializaton i, this is not actually the code I am using, just for getting the point across – Randi Sep 18 '17 at 16:13
  • @Randi: Hum. My answer still stands. Out of interest, it's better to use x * x to evaluate the square of x rather than one of the std::pow overloads. – Bathsheba Sep 18 '17 at 16:23
  • @Bathsheba, Thank you for the answer, but I am still confused. I am outputting the values to a file using [max_digits](http://en.cppreference.com/w/cpp/types/numeric_limits/max_digits10). I don't understand why they would be different. x^2 is just an example, the actual cost functions are more complicated, involving pixel differences of images. When I look at the different returned values form the cost function they are all the same (in the printed manner mentioned above) – Randi Sep 18 '17 at 16:31
  • How are you looking at those values? – Bathsheba Sep 18 '17 at 16:42
  • `std::cout << std::setprecision(std::numeric_limits::max_digits10) << y[i] << std::endl;` T = double – Randi Sep 18 '17 at 17:01
  • @Randi, forgive the arrogance (?!), but my answer still stands. The problem, I conject, is elsewhere. But don't blame floating point relational operators. They are fine, trust me.I suggest you re-ask - the question in it's current form is still useful - with more specific criteria; an answer will emerge quite quickly. (Just two doubles, the recipe for constucting them, and the faulty comparison will be sufficient.) Note that your numbers are hardly *subnormal* so there's nothing funky there to consider. – Bathsheba Sep 18 '17 at 18:46
  • You were right @Bathsheba. What I didn't mention in the post above is that the actual cost function is partially computed in CUDA, which was the root of the [numerical issues](http://docs.nvidia.com/cuda/floating-point/index.html). Once I took care of those, the comparison operators did what they are supposed to do. Thank you. – Randi Sep 19 '17 at 13:16

2 Answers2

4

Yes it is, assuming IEEE754 for your floating point types. Any two double values x and y say are such that exactly one of x < y, x == y, or x > y holds, with the exception of some edge cases like +Inf, -Inf, or NaN.

Confusion starts to arise when you use decimal notation to represent floating point values; e.g. there is no such double as 62536.5477752959 (or any other one on your list for that matter).

The numbers that you present have been truncated by your debugger / standard outputter, they are not the ones being used in the actual algorithm that you present. Be assured that the same decimal number always produces the same double, there is no arbitrary choice being made here: IEEE754 mandates that the closest double is picked.

For further reading, see Is floating point math broken?

Finally, replace int i with int i = 0. Currently the behaviour of your program is undefined.

Bathsheba
  • 231,907
  • 34
  • 361
  • 483
  • This doesn't really explain why the output of the program varies. We should assume that parsing the exact same string representation of a decimal number will produce the exact same `double`, hence from your three alternatives `x == y` should hold between numbers `0` and `7`. What is making the program think that `7` is smaller _sometimes_? – jdehesa Sep 18 '17 at 15:58
  • 1
    @jdehesa: Oops, yes I missed that. Have amended the answer. – Bathsheba Sep 18 '17 at 16:01
  • @Bathsheba, the values are being printed with the max # of digits possible. e.g. `std::cout << std::setprecision(std::numeric_limits::max_digits10) << y_[i] << std::endl;` where `T = double`. I assume that the same number of digits and the same values for each digit would lead as you said to the same output. But that is not what I am seeing. – Randi Sep 18 '17 at 18:47
  • I've addressed this in the question comments. – Bathsheba Sep 18 '17 at 18:48
2

You can use libraries to sidestep the behaviour/limitations of built-in types:

Live On Coliru

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <algorithm>
#include <iostream>

int main() {
    using Float = boost::multiprecision::cpp_dec_float_100;
    std::vector<Float> values = {
        Float{ "62536.5477752959" }, Float{ "62536.4840613718" }, Float{ "62536.4576412381" }, Float{ "62522.8487197062" },
        Float{ "62536.5473896233" }, Float{ "62536.5467941254" }, Float{ "62527.3508907998" }, Float{ "62536.5477752959" },
        Float{ "62517.5900098039" }, Float{ "62536.5477752959" },
    };

    auto hilo = std::minmax_element(values.begin(),values.end());

    std::cout << "low: " << *hilo.first << " high: " << *hilo.second << std::endl;

}

Prints

low: 62517.6 high: 62536.5

To print indexes:

// indexes:
std::cout << "low: " << (hilo.first-values.begin()) << " high: " << (hilo.second-values.begin()) << std::endl;
sehe
  • 374,641
  • 47
  • 450
  • 633
  • The goal seems to be to find the index of max and min values, not the values themselves. Also, even though using a library may help, it doesn't really answers why the original program does not work as expected in some cases. – jdehesa Sep 18 '17 at 16:00
  • 1
    @jdehesa see edit - I'm aware that I didn't answer that, I was helpful in addition to the standard answers and comments. It's a beat-to-death subject but solutions are sometimes required – sehe Sep 18 '17 at 16:03
  • @sehe: I have you down as being in the top 5 smartest contributors on this site. But your answer seems a little incomplete to me. Do feel free to grab the starting text from mine - I'm sure you've learnt nothing from it! – Bathsheba Sep 18 '17 at 16:09
  • 1
    On the contrary, @Bathsheba. I consider my answer a mere appendix to yours :) – sehe Sep 18 '17 at 16:10
  • 1
    @sehe. Have an upvote all the same. To me though the humble `double` gets a bad press. It's a wonderful construct, capable of representing a Planck length, the diameter of the universe, and the distance from here to Pluto to the nearest centimetre. – Bathsheba Sep 18 '17 at 16:11
  • 1
    @Bathsheba Cheers! It's not the range that trips people up :) – sehe Sep 18 '17 at 16:12
  • @sehe thank you for the suggestion. I would like to avoid moving way from the build-in types if possible. – Randi Sep 19 '17 at 13:12