449,422 Members | 1,263 Online
Need help? Post your question and get tips & solutions from a community of 449,422 IT Pros & Developers. It's quick & easy.

# Comparing doubles

 P: n/a Hi everyone, To determine equality of two doubles a and b the following is often done: bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD ); } But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. Since an double consists of mantissa, exponent and sign I assume that the "best" way in this case would be to just "ignore" the last few (e.g. 4-8 bits) of the mantissa to check for equality (at least for the x86 architecture, which follow IEEE 754). Have someone here ever tried to implement a similar approach? If yes, which experiences have been made? However my first approach to implement this failed. Can anyone tell what I might have done wrong? union IEEErepresentation { long long man : 52; int exp : 11; int sign : 1; }; union IEEEdouble { double d; IEEErepresentation c; }; const int number_of_bits_to_ignore = 8; bool isEqual (double a, double b) { IEEEdouble aa, bb; aa.d = a; bb.d = b; return ((aa.c.man << number_of_bits_to_ignore) == (bb.c.man << number_of_bits_to_ignore)); } int main() { string a = isEqual (1324.5678901234, 1324.5678901256 ) ? "true" : "false"; string b = isEqual (134.5678901234, 134.5678901256 ) ? "true" : "false"; string c = isEqual (.5678901234, .5678901256 ) ? "true" : "false"; string d = isEqual (1324.5678901234, 124.5678901256 ) ? "true" : "false"; string e = isEqual (134.5678901234, 134.5178901256 ) ? "true" : "false"; string f = isEqual (.5678901234, .3678901256 ) ? "true" : "false"; cout << "a = " << a << endl; cout << "b = " << b << endl; cout << "c = " << c << endl; cout << "d = " << d << endl; cout << "e = " << e << endl; cout << "f = " << f << endl; return 0; } Thanks in advance, Thomas Kowalski Jul 9 '07 #1
27 Replies

 P: n/a Thomas Kowalski a écrit : Hi everyone, To determine equality of two doubles a and b the following is often done: bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD ); } But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. std::numeric_limits::epsilon() is a good starting point. Since an double consists of mantissa, exponent and sign I assume that the "best" way in this case would be to just "ignore" the last few (e.g. 4-8 bits) of the mantissa to check for equality (at least for the x86 architecture, which follow IEEE 754). Have someone here ever tried to implement a similar approach? If yes, which experiences have been made? I did try to make an switchsign() fuction by reseting the SIGN bit. It was less efficient than a=-a; And you realize of course that you may succeed on your compiler and architecture but it will likely fail on others'. > However my first approach to implement this failed. Can anyone tell what I might have done wrong? union IEEErepresentation { long long man : 52; int exp : 11; int sign : 1; }; Here, you must mean struct. > union IEEEdouble { double d; IEEErepresentation c; }; const int number_of_bits_to_ignore = 8; bool isEqual (double a, double b) { IEEEdouble aa, bb; aa.d = a; bb.d = b; return ((aa.c.man << number_of_bits_to_ignore) == (bb.c.man << number_of_bits_to_ignore)); } Should be >>. Little endian inverts bytes, not bits Michael Jul 9 '07 #2

 P: n/a On Jul 9, 12:07 pm, Thomas Kowalski

 P: n/a "Thomas Kowalski" wrote: Hi everyone, To determine equality of two doubles a and b the following is often done: bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD ); } But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. Since an double consists of mantissa, exponent and sign I assume that the "best" way in this case would be to just "ignore" the last few (e.g. 4-8 bits) of the mantissa to check for equality (at least for the x86 architecture, which follow IEEE 754). Have someone here ever tried to implement a similar approach? If yes, which experiences have been made? I think bit fiddling should be a last resort, it forces anyone encountering the code to go through the same torturous process you went through when you coded it. Have you looked to see what you can do with the notion of "relative error"? Jul 9 '07 #4

 P: n/a Why is this a union?> Right, this should of course be a struct. Regards, Tk Jul 9 '07 #5

 P: n/a On Jul 9, 2:43 pm, "osmium" ::epsilon() ); } Jul 9 '07 #6

 P: n/a But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. std::numeric_limits::epsilon() is a good starting point. Does it work for this test set? Maybe I did something wrong, but for me its not working well. ASSERT ( isEqual (1234567890123456, 1234567890123457 ) ); ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) ); ASSERT ( isEqual (0.000001234567890123456, 0.000001234567890123457 ) ); I did try to make an switchsign() fuction by reseting the SIGN bit. It was less efficient than a=-a; And you realize of course that you may succeed on your compiler and architecture but it will likely fail on others'. I do. It's a pity that there is no better way to do this (or at least I don't know it). Here, you must mean struct. Yes, true. Should be >>. Little endian inverts bytes, not bits With ">>" its not really working either. I assume there is a problem with shifting a 52 bit number? Regards, Thomas Jul 9 '07 #7

 P: n/a Thomas Kowalski a écrit : On Jul 9, 2:43 pm, "osmium" Have you looked to see what you can do with the notion of "relative error"? Something like the following don't work really well, either for small numbers (e.g. 0.0000031234434). bool isEqual ( double a, double b ) { return ( fabs ((a-b)/b < std::numeric_limits::epsilon() ); } If a = b + 10 and a!=b and b very big, the above assumption will be true but a!=b. The equality can be: abs(x - y) < std::numeric_limits::epsilon() And inequality less or equal can be of the same format without abs: x - y < std::numeric_limits::epsilon(); Michael Jul 9 '07 #8

 P: n/a Thomas Kowalski a écrit : >>But this a approach usually fails if comparing doubles of differentmagnitude since it's hard or not possible to find a suitable thresholdvalue. std::numeric_limits::epsilon() is a good starting point. Does it work for this test set? Maybe I did something wrong, but for me its not working well. ASSERT ( isEqual (1234567890123456, 1234567890123457 ) ); ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) ); ASSERT ( isEqual (0.000001234567890123456, 0.000001234567890123457 ) ); On my system, epsilon is 2.22045e-16. It is normal the first two assert trigger a fault, the third not. > >I did try to make an switchsign() fuction by reseting the SIGN bit. Itwas less efficient than a=-a;And you realize of course that you may succeed on your compiler andarchitecture but it will likely fail on others'. I do. It's a pity that there is no better way to do this (or at least I don't know it). >Here, you must mean struct. Yes, true. >Should be >>. Little endian inverts bytes, not bits With ">>" its not really working either. I assume there is a problem with shifting a 52 bit number? No reason for that. Perhaps you can try simply to mask the double: union dl { double d; unsigned long l; }; dl d1=0.000001234567890123456; dl d2=0.000001234567890123457; unsigned long mask = (~0)<

 P: n/a Michael DOUBEZ wrote: Thomas Kowalski a écrit : >>>But this a approach usually fails if comparing doubles of differentmagnitude since it's hard or not possible to find a suitablethreshold value.std::numeric_limits::epsilon() is a good starting point. Does it work for this test set? Maybe I did something wrong, but forme its not working well.ASSERT ( isEqual (1234567890123456, 1234567890123457 ) );ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) );ASSERT ( isEqual (0.000001234567890123456,0.000001234567890123457 ) ); On my system, epsilon is 2.22045e-16. It is normal the first two assert trigger a fault, the third not. >> >>I did try to make an switchsign() fuction by reseting the SIGN bit.It was less efficient than a=-a;And you realize of course that you may succeed on your compiler andarchitecture but it will likely fail on others'. I do. It's a pity that there is no better way to do this (or at leastI don't know it). >>Here, you must mean struct. Yes, true. >>Should be >>. Little endian inverts bytes, not bits With ">>" its not really working either. I assume there is a problemwith shifting a 52 bit number? No reason for that. Perhaps you can try simply to mask the double: union dl { double d; unsigned long l; }; dl d1=0.000001234567890123456; dl d2=0.000001234567890123457; unsigned long mask = (~0)<. V -- Please remove capital 'A's when replying by e-mail I do not respond to top-posted replies, please don't ask Jul 9 '07 #10

 P: n/a Victor Bazarov wrote: Michael DOUBEZ wrote: >Thomas Kowalski a écrit : >>>>But this a approach usually fails if comparing doubles ofdifferent magnitude since it's hard or not possible to find asuitable threshold value.std::numeric_limits::epsilon() is a good starting point.Does it work for this test set? Maybe I did something wrong, but forme its not working well.ASSERT ( isEqual (1234567890123456, 1234567890123457 ) );ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) );ASSERT ( isEqual (0.000001234567890123456,0.000001234567890123457 ) ); On my system, epsilon is 2.22045e-16. It is normal the first twoassert trigger a fault, the third not. >>>I did try to make an switchsign() fuction by reseting the SIGN bit.It was less efficient than a=-a;And you realize of course that you may succeed on your compiler andarchitecture but it will likely fail on others'.I do. It's a pity that there is no better way to do this (or atleast I don't know it). Here, you must mean struct.Yes, true.Should be >>. Little endian inverts bytes, not bitsWith ">>" its not really working either. I assume there is a problemwith shifting a 52 bit number? No reason for that. Perhaps you can try simply to mask the double:union dl{ double d; unsigned long l;};dl d1=0.000001234567890123456;dl d2=0.000001234567890123457;unsigned long mask = (~0)<. Of course, I've forgotten to mention that exponents need to be the same, IOW the mantissas obtained from 'frexp' need to be corrected for the minimal exponent before comparison... V -- Please remove capital 'A's when replying by e-mail I do not respond to top-posted replies, please don't ask Jul 9 '07 #11

 P: n/a Victor Bazarov wrote: Victor Bazarov wrote: >Michael DOUBEZ wrote: >>Thomas Kowalski a écrit :>But this a approach usually fails if comparing doubles of>different magnitude since it's hard or not possible to find a>suitable threshold value.std::numeric_limits::epsilon() is a good starting point.Does it work for this test set? Maybe I did something wrong, butfor me its not working well.ASSERT ( isEqual (1234567890123456, 1234567890123457 ) );ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) );ASSERT ( isEqual (0.000001234567890123456,0.000001234567890123457 ) );On my system, epsilon is 2.22045e-16. It is normal the first twoassert trigger a fault, the third not. I did try to make an switchsign() fuction by reseting the SIGNbit. It was less efficient than a=-a;>And you realize of course that you may succeed on your compilerand architecture but it will likely fail on others'.I do. It's a pity that there is no better way to do this (or atleast I don't know it). Here, you must mean struct.Yes, true.Should be >>. Little endian inverts bytes, not bitsWith ">>" its not really working either. I assume there is aproblem with shifting a 52 bit number?No reason for that. Perhaps you can try simply to mask the double:union dl{ double d; unsigned long l;};dl d1=0.000001234567890123456;dl d2=0.000001234567890123457;unsigned long mask = (~0)<. Of course, I've forgotten to mention that exponents need to be the same, IOW the mantissas obtained from 'frexp' need to be corrected for the minimal exponent before comparison... Damn... I haven't paid attention before posting, have I? There is no need to correct the mantissas. Just compare the exponents and if they are not the same, you have your answer, if they are the same, compare the mantissas using the chosen epsilon. Sorry about the fuss... V -- Please remove capital 'A's when replying by e-mail I do not respond to top-posted replies, please don't ask Jul 9 '07 #12

 P: n/a Thomas Kowalski wrote: To determine equality of two doubles a and b the following is often done: bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD ); } But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. It is indeed true that the classic comparison fabs(a-b)

 P: n/a Juha Nieminen wrote: Thomas Kowalski wrote: >To determine equality of two doubles a and b the following is oftendone:bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD );}But this a approach usually fails if comparing doubles of differentmagnitude since it's hard or not possible to find a suitablethreshold value. It is indeed true that the classic comparison fabs(a-b) The magnitude of a value (ie. how many digits it has at the left of the decimal point) can be calculated with floor(log10(fabs(value))). Right. Of *a* value. Now if you divide the value by pow(10, magnitude) you get a value which is scaled to the range 0-1. Actually, it's scaled to the range 0.1-1. In other words: svalue = value/pow(10, floor(log10(fabs(value)))); Now svalue is value scaled to 0-1, 0.1-1 which makes it easy to compare the n significant digits of the value. Yes. Compare to what? You can divide the other value with that same divisor and then do the fabs(svalue1-svalue2)

 P: n/a On 9 Jul, 12:12, Michael DOUBEZ ::epsilon() is a good starting point. IMO its useful to be able to supply the epsilon . Being given a hard coded one seems arbitrary. You can also get more elaborate: #include //compare against some acceptable error range // epsilon can be computed // and may therefore be negative // hence use of abs inline bool eq_( double lhs, double rhs, double epsilon) { return std::abs(lhs-rhs) <= std::abs(epsilon); } //overload for function to supply epsilon inline bool eq_( double lhs, double rhs, double (*f)() ) { return std::abs(lhs-rhs) <= std::abs(f()); } //overload for function to supply epsilon based on inputs inline bool eq_( double lhs, double rhs, double (*f)(double, double) ) { return std::abs(lhs-rhs) <= std::abs(f(lhs,rhs)); } // get some percent based on x,y template double exp10range( double x, double y) { double avg =(std::abs(x) + std::abs(y))/2; return avg * std::pow(10.,N); } #include int main() { double x = -7.; double y = -7.007; std::cout << eq_(x,y, exp10range<-1) << '\n'; std::cout << eq_(x,y, exp10range<-2) << '\n'; std::cout << eq_(x,y, exp10range<-3) << '\n'; std::cout << eq_(x,y, exp10range<-4) << '\n'; } regards Andy Little Jul 9 '07 #15

 P: n/a kwikius a écrit : On 9 Jul, 12:12, Michael DOUBEZ Thomas Kowalski a écrit : >>Hi everyone,To determine equality of two doubles a and b the following is oftendone:bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD );}But this a approach usually fails if comparing doubles of differentmagnitude since it's hard or not possible to find a suitable thresholdvalue. std::numeric_limits::epsilon() is a good starting point. IMO its useful to be able to supply the epsilon . Being given a hard coded one seems arbitrary. You can also get more elaborate: [snip] I see what you mean but in the case of a customizable threshold, I wouldn't call it an equality but a proximity. All boils down to "what do the OP mean by isEqual(double,double) ?" as Victor Bazarov pointed out. Michael Jul 10 '07 #16

 P: n/a It is indeed true that the classic comparison fabs(a-b)

 P: n/a It is indeed true that the classic comparison fabs(a-b)

 P: n/a Hello everyone, thanks to the input I got in this thread I figured out some possible solution for the problem might be something like this: const double THRESHOLD = 4096.0*std::numeric_limits::epsilon(); // or any other multiple of epsilon. bool isEqual(double a, double b) { return fabs(1.0-(b/a)) < THRESHOLD; } or in case that a might be zero: bool isEqual(double a, double b) { if (a = 0.0) { return fabs(b) < THRESHOLD; } else { return fabs(1.0-(b/a)) < THRESHOLD; } } The question on hand is now. Does anyone have an idea whether and how this (or something similiar) could be speeded up (using C, not resorting to assembler) ? Regards, Thomas Kowalski Jul 10 '07 #19

 P: n/a On 10 Jul, 07:48, Michael DOUBEZ

 P: n/a Thomas Kowalski wrote: Hi everyone, To determine equality of two doubles a and b the following is often done: bool isEqual ( double a, double b ) { return ( fabs (a-b) < THRESHOLD ); } But this a approach usually fails if comparing doubles of different magnitude since it's hard or not possible to find a suitable threshold value. Since an double consists of mantissa, exponent and sign I assume that the "best" way in this case would be to just "ignore" the last few (e.g. 4-8 bits) of the mantissa to check for equality (at least for the x86 architecture, which follow IEEE 754). Have someone here ever tried to implement a similar approach? If yes, which experiences have been made? However my first approach to implement this failed. Can anyone tell what I might have done wrong? union IEEErepresentation { long long man : 52; int exp : 11; int sign : 1; }; union IEEEdouble { double d; IEEErepresentation c; }; const int number_of_bits_to_ignore = 8; bool isEqual (double a, double b) { IEEEdouble aa, bb; aa.d = a; bb.d = b; return ((aa.c.man << number_of_bits_to_ignore) == (bb.c.man << number_of_bits_to_ignore)); } int main() { string a = isEqual (1324.5678901234, 1324.5678901256 ) ? "true" : "false"; string b = isEqual (134.5678901234, 134.5678901256 ) ? "true" : "false"; string c = isEqual (.5678901234, .5678901256 ) ? "true" : "false"; string d = isEqual (1324.5678901234, 124.5678901256 ) ? "true" : "false"; string e = isEqual (134.5678901234, 134.5178901256 ) ? "true" : "false"; string f = isEqual (.5678901234, .3678901256 ) ? "true" : "false"; cout << "a = " << a << endl; cout << "b = " << b << endl; cout << "c = " << c << endl; cout << "d = " << d << endl; cout << "e = " << e << endl; cout << "f = " << f << endl; return 0; } Thanks in advance, Thomas Kowalski The lcc-win32 compiler system provides fcmp: --------------------------------------------------------cut here /* fcmp Copyright (c) 1998-2000 Theodore C. Belding University of Michigan Center for the Study of Complex Systems This file is part of the fcmp distribution. fcmp is free software; you can redistribute and modify it under the terms of the GNU Library General Public License (LGPL), version 2 or later. This software comes with absolutely no warranty. See the file COPYING for details and terms of copying. File: fcmp.c Description: see fcmp.h and README files. Knuth's floating point comparison operators, from: Knuth, D. E. (1998). The Art of Computer Programming. Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley. Section 4.2.2, p. 233. ISBN 0-201-89684-2. Input parameters: x1, x2: numbers to be compared epsilon: determines tolerance epsilon should be carefully chosen based on the machine's precision, the observed magnitude of error, the desired precision, and the magnitude of the numbers to be compared. See the fcmp README file for more information. This routine may be used for both single-precision (float) and double-precision (double) floating-point numbers. */ #ifdef HAVE_CONFIG_H #include #endif #include "fcmp.h" #include int EXPORT fcmp(double x1, double x2, double epsilon) { int exponent; double delta; double difference; /* Get exponent(max(fabs(x1), fabs(x2))) and store it in exponent. */ /* If neither x1 nor x2 is 0, */ /* this is equivalent to max(exponent(x1), exponent(x2)). */ /* If either x1 or x2 is 0, its exponent returned by frexp would be 0, */ /* which is much larger than the exponents of numbers close to 0 in */ /* magnitude. But the exponent of 0 should be less than any number */ /* whose magnitude is greater than 0. */ /* So we only want to set exponent to 0 if both x1 and */ /* x2 are 0. Hence, the following works for all x1 and x2. */ frexp(fabs(x1) fabs(x2) ? x1 : x2, &exponent); /* Do the comparison. */ /* delta = epsilon * pow(2, exponent) */ /* Form a neighborhood around x2 of size delta in either direction. */ /* If x1 is within this delta neighborhood of x2, x1 == x2. */ /* Otherwise x1 x2 or x1 < x2, depending on which side of */ /* the neighborhood x1 is on. */ delta = ldexp(epsilon, exponent); difference = x1 - x2; if (difference delta) return 1; /* x1 x2 */ else if (difference < -delta) return -1; /* x1 < x2 */ else /* -delta <= difference <= delta */ return 0; /* x1 == x2 */ } #ifdef TEST #include #include #define ASSERT assert #include int main() { int result; ASSERT ( fcmp (1234567890123456, 1234567890123457 ,DBL_EPSILON) ); ASSERT ( fcmp (1.234567890123456, 1.234567890123457 ,DBL_EPSILON) ); ASSERT ( fcmp (0.000001234567890123456, 0.000001234567890123457,DBL_EPSILON ) ); return 0; } #endif ---------------------------------------------------------end of fcmp.c Jul 10 '07 #21

 P: n/a On Tue, 10 Jul 2007 01:47:00 -0700, kwikius wrote: On 10 Jul, 07:48, Michael DOUBEZ All boils down to "what do the OP mean by isEqual(double,double) ?" asVictor Bazarov pointed out. Suggested function names: in_the_right_field(x,y); //courtesy of an old Boatbuilder. It was his favourite saying. in_the_ballpark(x,y); so_far_off_the_line_you_cant_even_see_the_line(x,y ) // Joey in "Friends" on_another_planet(x,y); // usually used by others to describe myself ;-) void nail_him_Tony(x,y); // Favourite saying of a roofer I knew when a batten (or whatever) was in the right spot, according to him anyway;-). Actually was "Nail 'im Tone' ! ". Tony was his long suffering Boss :-) I have a series of functions nad_eq, nad_gt, nad_lt, ... as in "near-as- dammit". Also has (possibly quite appropriate) echoes of "nads" as in "a kick in the 'nads" (sorry for that). -- Lionel B Jul 10 '07 #22

 P: n/a On 10 Jul, 09:27, Thomas Kowalski ::epsilon(); // or any other multiple of epsilon. could do : template double threshold() { enum { mux = (1 << N) }; return mux * std::numeric_limits::epsilon(); } regards Andy Little Jul 10 '07 #23

 P: n/a On Jul 10, 9:36 am, Thomas Kowalski

 P: n/a On 2007-07-10 03:36:40 -0400, Thomas Kowalski To be more a bit more precise. I have two numbers of similar magnitude (nearly the same exponent, different mantissa). I perform some calculations like e.g. additions and multiplications. If I want to check whether the number is in some range (let's say 0.0 to 1.0 ... or any other arbitrary numbers) it's clear that thanks to rounding errors I might get the wrong result from time to time (with a few billion numbers it's nearly always we case). Therefore I need a way to just ignore the last few digits of the mantissa to eliminate rounding errors while I want to preserve the accuracy of the _isEqual_ operation ( the operation isEqual is expected to tell me whether the numbers I am comparing would be the same in case the calculation (additions etc.) would be precise). Now all that's left is to figure out which of the "last few digits" reflect rounding errors and which ones are actually part of the "precise" result. If you don't have a rigorous procedure for doing that, the results you get from jamming zeros (or ones, for that matter) into the "last few digits" are no better than the original values. Doing that requires analyzing your algorithm to see where roundoff errors are occurring. Having done that, you no longer need this sort of sledge hammer. -- -- Pete Roundhouse Consulting, Ltd. (www.versatilecoding.com) Author of "The Standard C++ Library Extensions: a Tutorial and Reference (www.petebecker.com/tr1book) Jul 10 '07 #25

 P: n/a On 2007-07-10 03:23:52 -0400, Thomas Kowalski That's exactly what I mean. In my case I can't tell the magnitude of a number. (Lets just define the magnitude as the exponent of the double, since clearly 2 doubles are not equal if the exponent is not equal). But there's no such requirement in the C++ standard, and in general, it's not true. 1.0x10^3 is equal to 0.1x10^4, despite their different exponents. Granted, most floating-point representations these days use normalized values, so representable values with different magnitudes have different values. But that's hardware specific, not guaranteed. -- -- Pete Roundhouse Consulting, Ltd. (www.versatilecoding.com) Author of "The Standard C++ Library Extensions: a Tutorial and Reference (www.petebecker.com/tr1book) Jul 10 '07 #26

 P: n/a On Jul 10, 9:23 am, Thomas Kowalski

 P: n/a On 9 Jul, 19:42, Juha Nieminen

### This discussion thread is closed

Replies have been disabled for this discussion.

### Similar topics

Browse more C / C++ Questions on Bytes