## 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.
"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++):

# 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
"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
"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

(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
"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
```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.

```
Dec 04 2013
"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
"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