Help, best way to compare doubles
Author 
Message 
Thomas Koen #1 / 31

Help, best way to compare doubles
[ FollowupTo: comp.lang.c.moderated ]
Quote:
>: Having a hard time to figure out a function on how to compare doubles. >: Any ideas? >No function necessary.. C has this really cool boolean equality operator >that does this sort of thing for you. Have a look:
[...] It is unwise to compare floating point numbers for equality without thinking hard first. Roundoff error means that mathematically equivalent expressions may yield different results for different orderings; for example, a+b+c may very well compare unequal to a+c+b (try this out with a=1.0, b=1.0 and c=1.E20, for example). Roundoff error can also mean that algorithms which should converge don't, so beware. (I even know one system where a/2*2 != a for one quarter of floating point numbers  no points for guessing it's a /370 :) The correct way to compare two floating point numbers for "equality" to use fabs(ab) < eps. The problem arises in determining what to use for eps, of course :). What's reasonable depends very much on your problem and on the floating point characteristics of your system (for example, do your values span many orders of magnitude? Are they close to unity? Could they have opposing signs?), so it's really impossible to tell without a fairly thorough analysis. 
The joy of engineering is to find a straight line on a double logarithmic diagram.

Wed, 01 Jul 1998 03:00:00 GMT 


Carlie Coa #2 / 31

Help, best way to compare doubles
Quote:
>[ FollowupTo: comp.lang.c.moderated ]
>>: Having a hard time to figure out a function on how to compare doubles. >>: Any ideas? ... >The correct way to compare two floating point numbers for "equality" to >use fabs(ab) < eps..
Or to use _normalized_ differences. First version (has a division and a square root): fabs( a  b )/sqrt( a*a + b*b ) < eps better (has no division): fabs( a  b ) < eps * sqrt( a*a + b*b ) better yet (neither division nor square root): ( t = a  b )*t < eps * ( a*a + b*b )

Thu, 02 Jul 1998 03:00:00 GMT 


Pierre Assel #3 / 31

Help, best way to compare doubles
Quote:
>>: Having a hard time to figure out a function on how to compare doubles. >>: Any ideas? >The correct way to compare two floating point numbers for "equality" to >use fabs(ab) < eps. The problem arises in determining what to use for >eps, of course :).
Sometimes I do this: if(tmp1= a+b, tmp2= (ab)/64, tmp2+= tmp1, tmp1==tmp2) { /* about equal */; } else { /* different */; } The idea is to test whether (ab) falls off the mantissa when scaled to the same exponent as (a+b). The 64 is for slop. Hopefully this code won't be placed in the middle of a tight loop... I believe the C compiler is forbidden to optimize the code into a mathematically equivalent but numerically incorrect test. Can it be done without the temporaries ? I don't read the standard anymore :)  Pierre Asselin, Westminster, Colorado

Thu, 02 Jul 1998 03:00:00 GMT 


Gene Wirchen #4 / 31

Help, best way to compare doubles
Quote:
>>[ FollowupTo: comp.lang.c.moderated ]
>>>: Having a hard time to figure out a function on how to compare doubles. >>>: Any ideas? >... >>The correct way to compare two floating point numbers for "equality" to >>use fabs(ab) < eps.. >Or to use _normalized_ differences. >First version (has a division and a square root): > fabs( a  b )/sqrt( a*a + b*b ) < eps >better (has no division): > fabs( a  b ) < eps * sqrt( a*a + b*b ) >better yet (neither division nor square root): > ( t = a  b )*t < eps * ( a*a + b*b )
Nice. Is there a general method you are using here? Details would be nice. If you couch it in terms of C, the moderator may whistle and look the other way <G>. Free lessons on whistling and averting ones eyes available here. Sincerely, Gene Wirchenko C Pronunciation Guide: y=x++; "wye equals ex plus plus semicolon" x=x++; "ex equals ex doublecross semicolon"

Fri, 03 Jul 1998 03:00:00 GMT 


k.. #5 / 31

Help, best way to compare doubles
Quote:
> Or to use _normalized_ differences. > First version (has a division and a square root): > fabs( a  b )/sqrt( a*a + b*b ) < eps > better (has no division): > fabs( a  b ) < eps * sqrt( a*a + b*b ) > better yet (neither division nor square root): > ( t = a  b )*t < eps * ( a*a + b*b )
^ Does any rule in C++ tell you that the latter t will be evaluated later than the expression "t = a  b"? [C++? Oh, well. I don't think there is one in C, anyway. mod]  /v\ Isaac To, karkeung  'y.a1j  MPhil, Computer Science  p:b>w>G u>G:S$h=R5{ 
URL: http://www.cs.hku.hk/~kkto  :t6!Ghttp://www.cs.hku.hk/~kkto  \^/

Fri, 03 Jul 1998 03:00:00 GMT 


Thomas Koen #6 / 31

Help, best way to compare doubles
Quote:
>>The correct way to compare two floating point numbers for "equality" to >>use fabs(ab) < eps.. >Or to use _normalized_ differences. >First version (has a division and a square root): > fabs( a  b )/sqrt( a*a + b*b ) < eps
This only works correctly if a*a + b*b does not overflow your machine number, and also that a*a+b*b does not underflow. Still some thought required... A cleaner overall solution would be to use interval arithmetic instead of conventional floating point variables. But by then, we'd probably be outside of what is normally considered the scope of C. 
The joy of engineering is to find a straight line on a double logarithmic diagram.

Fri, 03 Jul 1998 03:00:00 GMT 


Michael Smit #7 / 31

Help, best way to compare doubles
Quote:
> The correct way to compare two floating point numbers for "equality" to > use fabs(ab) < eps. The problem arises in determining what to use for > eps, of course :). What's reasonable depends very much on your problem > and on the floating point characteristics of your system (for example, > do your values span many orders of magnitude? Are they close to unity?
I often use macros like: #define TOLERENCE 1.0e6 #define EQUAL( a, b ) (fabs(ab)<=(fabs(a)/TOLERENCE)) By relating epsilon to the values used, the macro becomes more general. Note that the "<=" is important for the case a==b==0.  #####################################################################
Emmenjay Consulting Sofware Development System Administration/Management Applications Support PC Installation/Maintenance Training Technical Writing Windows/Unix General Consulting/Troubleshooting PO Box 909 Ph 018 240 704 Kensington 2033 AUSTRALIA #####################################################################

Fri, 03 Jul 1998 03:00:00 GMT 


Ed Ho #8 / 31

Help, best way to compare doubles
> > > Or to use _normalized_ differences. > > > First version (has a division and a square root): > > fabs( a  b )/sqrt( a*a + b*b ) < eps > > > better (has no division): > > fabs( a  b ) < eps * sqrt( a*a + b*b ) > > > better yet (neither division nor square root): > > ( t = a  b )*t < eps * ( a*a + b*b ) > ^ > > Does any rule in C++ tell you that the latter t will be evaluated later > than the expression "t = a  b"? > [C++? Oh, well. I don't think there is one in C, anyway. mod] > While we're at it, note that the final version (even after *that* bug is fixed) is _not_ equivalent to the ones which precede it. It *should* be something like: t = a  b; if ( t*t < eps*eps*(a*a + b*b) ) { /* 'a' and 'b' are "morally equal" */ } else { /* they're not */ } For small values of 'eps' (which are surely the most interesting), the original (flawed) approach is far too tolerant ... >  > /v\ > Isaac To, karkeung  'y.a1j  > MPhil, Computer Science  p:b>w>G u>G:S$h=R5{ 
> URL: http://www.cs.hku.hk/~kkto  :t6!Ghttp://www.cs.hku.hk/~kkto  > \^/  Ed Hook  Coppula eam, se non posit Computer Sciences Corporation  acceptera jocularum. NASA Langley Research Center  Me? Speak for my employer?...<*snort*>

Sun, 05 Jul 1998 03:00:00 GMT 


Adam Worra #9 / 31

Help, best way to compare doubles
>>> The correct way to compare two floating point numbers for "equality" to >>> use fabs(ab) < eps.. >> Or to use _normalized_ differences. >> First version (has a division and a square root): fabs( a  b )/sqrt( >> a*a + b*b ) < eps Thomas> This only works correctly if a*a + b*b does not overflow your Thomas> machine number, and also that a*a+b*b does not underflow. Still Thomas> some thought required... It also falls over if 'a' is less then 'eps', but b is exactly zero; you end up with a/a < eps , ==> 1.0 < eps. This is false, even though 'a' was within 'eps' of 'b'. I ended up with a big clunky macro for robustly dealing with all this, which uses two little macros for determining if a number is zero, (i.e. absolute error, (a<eps), and one for relative_error (as explained above). #define real_equal(a,b) \ ( ((a)==0.0) ? /* if (a==0.0) { */ \ ( ((b)==0.0) ? /* if (b==0.0) { */ \ TRUE /* return TRUE; */ \ : /* } else { */ \ is_zero (b) ) /* return is_zero(b); */ \ : /* } else { */ \ ( ((b)==0.0) ? /* if (b==0.0) { */ \ is_zero(a) /* return is_zero(a); */ \ : /* } else { */ \ eq_rel_error(a,b,EPSILON) /* return eq_rel_error; */ \ ) /* } */ \ ) /* } */ Does anybody have a better solution ? It seems a bit horrific !  Adam  =~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=
Computer Science Department, The University, Bristol, UK

Sun, 05 Jul 1998 03:00:00 GMT 


Paul Tombli #10 / 31

Help, best way to compare doubles
Quote:
>> The correct way to compare two floating point numbers for "equality" to >> use fabs(ab) < eps. The problem arises in determining what to use for >#define TOLERENCE 1.0e6 >#define EQUAL( a, b ) (fabs(ab)<=(fabs(a)/TOLERENCE))
I use a slightly more "expensive", and slightly less standard way. The advantage is that it works well for pairs of numbers no matter how large or small they are. I don't have the code I use here at home, but it's something like this. One difference between the real function I use and this one is that the real function has code to handle comparisons between a real zero and a very small positive or negative double, and where the exponents are off by one because one number is slightly less than a power of two and the other is slightly more. The function "frexp()" is a POSIX fuction that breaks a floating point number into a fraction and exponent. Boolean dequal(double x, double y) { int x_exp, y_exp; double x_mant, y_mant; x_mant = frexp(x, &x_exp); y_mant = frexp(y, &y_exp); if (x_exp != y_exp) return False; if (fabs(x_mant  y_mant) < 1.0e6) return True; return False; Quote: }

<a href="http://www.servtech.com/public/ptomblin/">My home page</a> "While you're waiting, read the free novel we sent you. It's a Spanish story about a guy named `Manual'"  Dilbert

Sun, 05 Jul 1998 03:00:00 GMT 


Michael J Ze #11 / 31

Help, best way to compare doubles
Quote:
>[ FollowupTo: comp.lang.c.moderated ]
>>: Having a hard time to figure out a function on how to compare doubles. >>: Any ideas? >>No function necessary.. C has this really cool boolean equality operator >>that does this sort of thing for you. Have a look: >It is unwise to compare floating point numbers for equality without >thinking hard first.
...as evidenced by the fact that many of the proposed "solutions" have potential overflow or underflow problems. Of course, it's often necessary to first make sure you've defined the algorithm for equality that you want, then make sure it works for all possible input. The algorithm you want is highly dependent on the application (whose context was either never posted or I missed it). For example, if you're doing a homework assignment and are waiting for a numerical calculation to stabilize, and you know more or less the range of input values for the comparison, then define an epsilon, use a simple test with fabs, and be done with it. If you're writing a database server though, you're likely to need a tunable comparison (one user wants 10 digits of precision, another wants 14), that's well documented, and works correctly for all input values. I can think of two broad categories of equality algorithms. I'll call them interval and difference, for want of knowing if they have official names. An interval method would partition the real numbers into contiguous subsets or intervals within which all real numbers would be defined as "equal" for the purpose of the aplication. For example, the interval called "1.65" might include all real numbers > 1.645 and <= 1.655. This has an advantage in doing some kinds of applications because it corresponds to how we often think of currency. In the US, we round monetary figures to the nearest .01 dollars, for example. The disadvantage of this method is that at all places in the application where conversion is done between these intervals and real numbers, there's a potential for unexpected behavior. 1.646 and 1.653 are considered "equal" but 1.646 and 1.644 are unequal. More specifically, a calculation that we expect to result in 1.645 could be labelled 1.64 or 1.65 depending on the order of the operands, etc. This method is often most appropriate for systems that convert input values once into a decimal format and performa calculations with them in that format. A difference method relies on the relative difference between two numbers. Instead of translating between real numbers and these intervals, the results are kept in the natural machine format. (Note however that the natural machine format already is an interval of the real numbers, but hopefully at a finer grain than what the application cares about.) When two of these numbers are compared, they are equal if their difference is less than some other number. This other number could be a fixed amount (like .01) or an amount relative to the numbers being compared (say, .01 times the number with the lowest magnitude). One advantage of this method is that there are no discontinuities like there is with the interval method. A disadvantage is that if A=B and B=C, A and C are not necessarily equal. One can even combine fixed and relative comparisons of the difference method for difference equality. It might be desireable to have a relative difference of .00001 times the number with the lowest magnitude for equality (essentially comparing to 5 digits), but for values close to 0 use a fixed difference like .00001. This ensures that results that you expect to be 0 compare as equal, even if one is twice the other. (For example, one result might be the lowest machine representable value, and anotehr result might be twice that value, yet you still probably want those to compare equal in many applications.) Once you've decided on the algorithm that's most appropriate for your application, you still have to implement it so that there's no chance of over or underflow. This can be nontrivial, especially if you also want to do it quickly in terms of runtime performance. Also remember that if you have a special function for equality of real numbers, you need to be careful if you're using comparison operators: if (a < b) might instead need to be written as: if ((!real_equal(a,b)) && (a<b)) to prevent unexpected results! As Thomas Koenig said, one should think carefully about what you need before answering this question. There is no single correct method of doing the comparison. Michael J. Zehr Sr. Software Engineer Kenan Systems Corporation

Sun, 05 Jul 1998 03:00:00 GMT 


n.. #12 / 31

Help, best way to compare doubles
Quote:
>>: Having a hard time to figure out a function on how to compare doubles. >>: Any ideas? >The correct way to compare two floating point numbers for "equality" to >use fabs(ab) < eps. The problem arises in determining what to use for >eps, of course :).
Back when I still did numerical programming (in fortran) I used the equilvalent of eps1 >= (fabs(ab) / fmax(fabs(a), fabs(b), eps2) If you are converging on a large number, then you want a relative equality test (agree to X decimal places), however, if the sequence happens to converge on zero than a relative test won't work, so below "eps2" you switch over to an absolute "fuzz factor"

Sun, 05 Jul 1998 03:00:00 GMT 


usene #13 / 31

Help, best way to compare doubles
Quote:
>>: Having a hard time to figure out a function on how to compare doubles. >>: Any ideas?
..>The correct way to compare two floating point numbers for "equality" to Quote: >use fabs(ab) < eps..
You'll be in trouble if you apply the transitive law to this: a approximately equal to b, b approximately equal to c will not imply a approximately equal to c. Any epsilon you choose in the absence of further information is likely to be arbitrary. If it's important that two floats are or are not equal, there's usually some other logical evidence (which you may have foolishly thrown away) elsewhere in the program telling you so. If there isn't any evidence, then more often than not it's not terribly important to tell the difference between the case of the floats being equal, and the case of the floats being closer together than epsilon, so use (fabs(ab) < eps) You should only mean what you say and say what you mean. If you can't rely on the statement (a == b) then don't say something else like (fabs(ab) < eps) instead, because you won't mean it, and you should only say what you mean, you see? Julian Todd.

Sun, 05 Jul 1998 03:00:00 GMT 


Peter Roozema #14 / 31

Help, best way to compare doubles
Quote:
> >>: Having a hard time to figure out a function on how to compare doubles. > >>: Any ideas? > ... > >The correct way to compare two floating point numbers for "equality" to > >use fabs(ab) < eps.. > Or to use _normalized_ differences. Better > First version (has a division and a square root): > fabs( a  b )/sqrt( a*a + b*b ) < eps
A usually faster variant (assuming absolute values are taken faster then multiplications: fabs( a  b )/( fabs(a) + fabs(b) ) < eps Avoiding the square also significantly reduces the chances of overflow or underflow during evaluation of the expressions. Try to compare the doubles 1e200 and 2e200 or 1e200 and 2e200 with both routines on a machine that uses IEEE double precision arithmetics. Squaring them will give overflow or underflow. In the file <float.h> some usefull constants are defined: DBL_MAX, DBL_MIN and especially DBL_EPSILON, the minimum positive floating point number x such that 1.0 + x != 1.0 In most cases you can chose eps a few times DBL_EPSILON when using normalised differences. Quote: > better (has no division): > fabs( a  b ) < eps * sqrt( a*a + b*b )
for completeness: fabs( a  b ) < eps * ( fabs(a) + fabs(b) ) Quote: > better yet (neither division nor square root): > ( t = a  b )*t < eps * ( a*a + b*b )
^ YUCK the order of evaluation is UNDEFINED does this t mean the NEW or the ORIGINAL value of t? (squaring argument remains also valid over here)  ===============================================================================
There's no computer like an ARMed computer! =============================================================================== ... Sorry... my mind has a few bad sectors.

Sun, 05 Jul 1998 03:00:00 GMT 


M. J. Saltzm #15 / 31

Help, best way to compare doubles
Slightly edited, Quote:
>>>The correct way to compare two floating point numbers for "equality" to >>>use fabs(ab) < eps.. >>Or to use _normalized_ differences. >> fabs( a  b )/sqrt( a*a + b*b ) < eps >> fabs( a  b ) < eps * sqrt( a*a + b*b ) >> ( t = a  b )*t < eps * ( a*a + b*b )
^^^^^^^^^^^^^^^ Others have pointed out that the assignment and reference to t results in undefined behavior, but I will anyway, in the hope that I can con the moderator into passing the rest of the article. Quote: > Nice. Is there a general method you are using here? Details >would be nice. If you couch it in terms of C, the moderator may >whistle and look the other way <G>.
To _normalize_ a vector (or a quantity computed from the components of a vector, as in this case) is to divide its components by the norm (or length) of the vector. If we treat the values (a,b) as a vector, its length by the usual "Euclidian" definition is sqrt(a^2 + b^2) (where ^ is used here to mean "raised to the power", a common notation in groups that are not comp.lang.c.*). The point is this: If you are interested in computing results to (say) six digits of accuracy, then you would want a = 100000 and b = 100000.5 to compare equal, but if you compare using fabs(a  b) < 1.0e6, they will not, whereas a = 0.1 and b = 0.1000005 will. The process of dividing by the norm results in the comparisons being made in a way that takes the number of digits of agreement into account while factoring out the influence of the magnitudes of the values. Some comments: (1) The Euclidian norm is not the only one that can be used. You can get the same effect (though with slightly different results for some comparisons) by using max(fabs(a),fabs(b)) or (fabs(a)+fabs(b)) as the norm. This saves a sqrt() call (at the cost of fabs() calls), and the "max norm" does not risk overflow. (2) You may not want to normalize comparisons of very small numbers. For example, if six digits is your goal, then you may want to accept 0.000001 and 0.0000015 as equal, but a normalized comparison will not do so. This situation is often dealt with by using max(1.0, fabs(a), fabs(b)) or (1.0 + fabs(a) + fabs(b)) as the "norm". I'm not sure if anything other than max(1.0, sqrt(a^2+b^2)) is needed to deal with the Euclidian version.  Matthew Saltzman Clemson University Math Sciences

Sun, 05 Jul 1998 03:00:00 GMT 


Page 1 of 3

[ 31 post ] 

Go to page:
[1]
[2] [3] 
