www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - floating point comparison basics

reply "ed" <sillymongrel gmail.com> writes:
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++):









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
next sibling parent "ed" <sillymongrel gmail.com> writes:
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
prev sibling parent reply "H. S. Teoh" <hsteoh quickfur.ath.cx> writes:
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
next sibling parent reply "ed" <sillymongrel gmail.com> writes:
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
parent "Gary Willoughby" <dev nomad.so> writes:
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.htm


This 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
prev sibling parent reply "bearophile" <bearophileHUGS lycos.com> writes:
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
parent "ed" <sillymongrel gmail.com> writes:
On Wednesday, 4 December 2013 at 01:52:09 UTC, bearophile wrote:

 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


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,
Ed
Dec 03 2013