## digitalmars.D.learn - floating point comparison basics

- ed (23/23) Dec 03 2013 Hi All,
- ed (26/26) Dec 03 2013 OK, I've had a look into the code in std.math finally understand
- H. S. Teoh (58/74) Dec 03 2013 The first rule of floating-point comparisons is that you never use ==.
- ed (15/15) Dec 03 2013 On Tuesday, 3 December 2013 at 23:17:29 UTC, H. S. Teoh wrote:
- Gary Willoughby (4/6) Dec 04 2013 This is the exact source i used to implement the assertApprox
- bearophile (9/11) Dec 03 2013 So I suggested to disallow the == among FP numbers in D (and
- ed (7/18) Dec 03 2013 Thanks, I will look into it.

Hi All, I'm learning programming and chose D because it is the best :D But, I've hit floating point numbers and I'm stuck on some of the basics. What is the proper way to do floating point comparisons, in particular I need to check if a value is zero? For example, given "real x = someCalculatingFunction();" how do I check if X is zero in a robust way. if(x == 0.0) {} // <-- Will this work as expected? I see some code from others in my class doing this (mostly C++): # if(fabs(x) < 1e-10) {} // <-- I've heard this is bad, is it? # if(fabs(x) < std::numeric_limits<double>::min()) {} # if(fabs(x) < std::numeric_limits<double>::epsilon()) {} # if(!(fabs(x) > 0.0)) {} # double eps = 1/(1/x); # if(!(fabs(eps) > 0.0)) {} So you can see we're all confused! Thanks, Ed PS: I have read this great article and the links it provides: http://dlang.org/d-floating-point.html Most of it makes sense but I'm struggling to tie it all together when it comes time to apply it.

Dec 03 2013

OK, I've had a look into the code in std.math finally understand it, I think. I was confused about maxRelativeError and maxAbsoluteError then I found this also: http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm (...If you want to count numbers near zero but of opposite sign as being equal then you need to add a maxAbsoluteError check also. ...) So my understanding is this: 1. If comparing to 0.0 then "fabs(lhs) < maxAbsoluteError" is OK 2. If comparing two floating point values X, Y: "fabs((lhs-rhs)/rhs) < maxRelativeError" is OK. Reading elsewhere on the web (bad idea I know) I get the impression the ordering of lhs and rhs matters. Some code, as in the link above, is written as: if(fabs(rhs) > fabs(lhs)) { relErr = fabs((lhs-rhs)/rhs); } else { fabs((lhs-rhs)/lhs); } if(relErr < maxRelErr) { return true;} else {return false;} Why does Phobos not check lhs < rhs? Does this imply that in some cases: approxEqual(X, Y) != approxEqual(Y, X) ?? Thanks, Ed

Dec 03 2013

On Tue, Dec 03, 2013 at 11:03:48PM +0100, ed wrote:Hi All, I'm learning programming and chose D because it is the best :D But, I've hit floating point numbers and I'm stuck on some of the basics. What is the proper way to do floating point comparisons, in particular I need to check if a value is zero?The first rule of floating-point comparisons is that you never use ==. Well, not *literally* never (there are some cases where it's useful), but you should never use == by default, and every time you do, you'd better have a good reason for it. As for why, see below.For example, given "real x = someCalculatingFunction();" how do I check if X is zero in a robust way. if(x == 0.0) {} // <-- Will this work as expected?Most likely, it will not. Unless you explicitly set x to 0.0 somewhere. If x is the result of some complex computations, most likely it will not be *exactly* 0.0. The correct way to compare floats is to write: if (abs(x - y) < Epsilon) { // x and y are approximately equal } for some small-enough value of Epsilon. Or, in your case: if (abs(x) < Epsilon) { // x is approximately zero } So the first thing to know about floating-point is that it's only an approximation, and because of that, (1) it is NOT the same as the real numbers in mathematics, and (2) operations on floating-point values do not always follow the same rules as mathematics. For example: float a = 1.0 / 5.0; assert(a == 0.2); // <--- this will FAIL This is because 1/5 in binary has a non-terminating digit expansion (much like 1/3 in decimal has a non-terminating digit expansion 0.3333...). Since we can only store a finite number of digits in a float, the digits have to be truncated past a certain point, and that introduces a slight round-off error. The round-off error introduced by the division operation in 1.0 / 5.0 is slightly different from the round-off error introduced by converting the literal 0.2 into binary, so 1.0 / 5.0 == 0.2 fails to hold. Another gotcha is that (x+y)+z is not always the same as x+(y+z), unlike in mathematics. If x is a very large number relative to y, (x+y) could be truncated to just x, so writing (x+y)+z becomes the same as writing x+z; but if z is of intermediate magnitude, then (y+z) could be a value different from z, and so x+(y+z) will produce a different answer than (x+y)+z. [...]PS: I have read this great article and the links it provides: http://dlang.org/d-floating-point.html Most of it makes sense but I'm struggling to tie it all together when it comes time to apply it.Two rules of thumb with floating-point values: (1) Never use == unless you have a good reason for it (and if you don't know what constitutes a good reason, you don't have one, so don't use it). Instead, compare abs(a-b) with a small constant value, conventionally named Epsilon, that represents "close enough" for your purposes (and this will differ depending on what you're trying to do in your program). (2) Don't assume that floating-point operations behave the same way as mathematical operators. For example, x*x - y*y and (x+y)*(x-y) are the same thing in math, but in floating-point arithmetic, the former is vulnerable to catastrophic cancellation (which may produce garbled results for certain inputs), whereas the latter will give a reasonably accurate answer for all inputs. When in doubt, consult well-researched resources on floating-point arithmetic, such as: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html T -- Truth, Sir, is a cow which will give [skeptics] no more milk, and so they are gone to milk the bull. -- Sam. Johnson

Dec 03 2013

On Tuesday, 3 December 2013 at 23:17:29 UTC, H. S. Teoh wrote: Thanks for the reply and the link, it gives me more confidence that I'm understanding things a bit better. I finally feel like I've opened the lid on the black box of float precision. I found a great set of articles which I'm working through now. http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm It seems that Phobos approxEqual is (one of the) agreed upon methods from all sources I've read. Another involves computing the ULPs as integers and comparing those for equality. Both approaches are very similar in terms of code though. Next year I have a whole semester on numerical computing, covering floating point numbers, machine imprecision and how to deal with etc....I'm determined to ace that subject :D Cheers, Ed

Dec 03 2013

On Tuesday, 3 December 2013 at 23:50:34 UTC, ed wrote:I found a great set of articles which I'm working through now. http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htmThis is the exact source i used to implement the assertApprox function in the DUnit toolkit. See the unit tests for examples. https://github.com/nomad-software/dunit/blob/master/source/dunit/toolkit.d#L42-L112

Dec 04 2013

H. S. Teoh:The first rule of floating-point comparisons is that you never use ==.So I suggested to disallow the == among FP numbers in D (and allow FP comparisons with "is" or with specialized functions, including one function that does what == does today). For the original poster, for general FP comparisons I suggest std.math.feqrel, it's not fast, it's not const-correct, but it's good. Bye, bearophile

Dec 03 2013

On Wednesday, 4 December 2013 at 01:52:09 UTC, bearophile wrote:H. S. Teoh:Thanks, I will look into it. I've implemented my comparison function and it gives the same results as approxEquals ... yay! ...only took me several hours to achieve about 10 LOC but I finally understand it Cheers, EdThe first rule of floating-point comparisons is that you never use ==.So I suggested to disallow the == among FP numbers in D (and allow FP comparisons with "is" or with specialized functions, including one function that does what == does today). For the original poster, for general FP comparisons I suggest std.math.feqrel, it's not fast, it's not const-correct, but it's good. Bye, bearophile

Dec 03 2013