# Choosing the right epsilon for comparing doubles

Discussion in 'C++' started by walkietalkie74, Feb 2, 2014.

1. ### walkietalkie74Guest

Hi all,

Sorry if this question has been already asked in the past (I bet so). I would appreciate if you could point me in the right direction.

I do not generally compare doubles because I am well aware of the round-offerrors that make such a comparison meaningless.

However, I am trying to write my program in a test-driven way, so I'd like to write a tester that checks that the result of a specific calculation is actually zero.

My code is the one below. It's calculating a straight line equation based on a point and a vector (in space, a straight line is the intersection of two planes, hence the getPlane1() and getPlane2()), and is checking whether the obtained equation is correct by calculating arbitrary points that are supposed to be on the same straight line, and checking whether their coordinates solve the plane equations:

bool GeometryTester::runTests() {

Point p(2, -2, 1);
Vector v(2, 4, 2);
StraightLine straightLine = Geometry::getStraightLine(p, v);

bool result = true;

for (double scalar = -10.0; scalar <= 10.0; scalar += 0.11131719) {
Point calcPoint = Geometry::getPoint(p, v, scalar);

double s1 = straightLine.getPlane1().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

double s2 = straightLine.getPlane2().solve(calcPoint.getX(),
calcPoint.getY(),
calcPoint.getZ());

result = result &&
MathComparer::equalsZero(s1);

result = result &&
MathComparer::equalsZero(s2);
}

return result;
}

For now, in MathComparer::equalsZero(), I am using a (wrong - I know) direct == comparison between the argument passed in and 0.0, but I'd like toreplace it with the appropriate comparison using an epsilon.

Am I correct in thinking that the epsilon should actually be chosen depending on the current calculation and the accuracy I am expecting from it?

I can see, for example, that the solve() method most of the times returns 0..0 but at times something like 6e-12 (which is absolutely sensible, given the imprecise nature of floating-pint types).

Would it be sensible to pass in an arbitrary epsilon to the equalsZero() function?

I have seen examples of usage of std::numerics<double>::epsilon but I am not sure how it should be used in my case...

Can you help me understand this stuff better?

Thanks!

walkietalkie74, Feb 2, 2014

2. ### JaxGuest

Walkietalkie as I'm sure you know there's more than one pair of planes
which can describe a given line..... so are there some special constraints
on the 2 planes you obtain which define the straight line? Just wondering.

Anyways you may want to consider always choosing the error (epsilon) to be
a few orders of magnitude greater than the specified inaccuracy of your
floating-point hardware. Alternatively, and I prefer this, epsilon could
be based on the magnitude of the numbers you are comparing rather than an
absolute constant.

Jax, Feb 2, 2014

3. ### walkietalkie74Guest

Yes, I know that there can be infinite pairs of planes which can describe agiven line. This is the actual algorithm used in the method that builds the equations of the planes, given a point and a vector:

StraightLine Geometry::getStraightLine(const Point& p, const Vector& v) {
StraightLine straightLine;

double vx = v.getvx();
double vy = v.getvy();
double vz = v.getvz();

double xp = p.getX();
double yp = p.getY();
double zp = p.getZ();

straightLine.setPlane1(Plane(vy, -vx, 0, vx*yp - vy*xp));
straightLine.setPlane2(Plane(0, vz, -vy, vy*zp - vz*yp));

return straightLine;
}

Do you see any errors?
I used the GeometryTester class in my previous post to test that the equations are correct. Almost all the points that I calculate are solutions for the equations, with the exceptions of a few for which the solve() functions returns something like 0.000000000006....(which should be an approximation of zero - though I agree that it's quite a big error).
Thanks. I have already discarded the idea of an absolute constant. I am aware of the existence of std::numerics<double>::epsilon and I've found a few examples around showing how to relate it to the magnitude of the numbers involved in the calculations.
Someone suggested to use it the following way:

a and b are "equal" if:
std::abs(a - b) <= std::abs(a) * std::numeric_limits<double>::epsilon

It seems using the same algorithm shown in section 4.2.2 of The Art of Computer Programming (D. Knuth), but it doesn't seem to work.

walkietalkie74, Feb 2, 2014
4. ### walkietalkie74Guest

Thanks for your answer. The problem with choosing an epsilon dictated by gut feeling is that it would make the tester class necessarily unreliable. I wouldn't be 100% sure that the epsilon would be correct in all circumstances, using any numbers. If a new developer adds a test case which fails, I have no way to figure out whether it's because of the wrong epsilon chosen orthe new test condition itself...
Thanks for the suggestions. Really precious. I am afraid you're right. Thatdifference seems quite a big error to me as well. I might well be doing what you say, because I've analysed the values at the intermediate steps of the algorithm and, for example, I've got something like the following:

38.398935600000016 -
3.3989356000000184 -
35

walkietalkie74, Feb 2, 2014
5. ### walkietalkie74Guest

Thanks. That's another option I'd thought about.
I will google it and see if I can find any libraries.

walkietalkie74, Feb 2, 2014
6. ### Stefan RamGuest

A floating point number, like say, 10000.1, has an IEEE
example. This actually represents a whole range of numbers
from 0x1.3880ccccccccc8p13 to 0x1.3880ccccccccd7p13 or
something like that.

This error of the plane parameters gives an error
(uncertainity) of the planes and of the line parameters. The
points of the line have to be points of the plane within the
limits of those errors. The details of error propagation
(propagation of distributions) are a standard subject of
numerical mathematics, I believe.

And strictly for every triple of a line and two planes,
we do not get a yes/no-answer but a probability that it
is correct.

But for any plane-like object in the physical world, we
cannot measure its parameters up to 15 digits, so the
measurement error usually exceeds this IEEE represenation
error. Possibly, this larger error also has to be taken
into consideration.

Stefan Ram, Feb 2, 2014
7. ### walkietalkie74Guest

Hi Stefan. Thanks. Yes, I know that the floating point values are inherently imprecise and that any calculation with them can propagate errors. I wouldn't normally compare two doubles for equality. This necessity arises from the fact that I need to write a tester class. My three-dimensional objects will eventually be projected and displayed, so I am not really interested in having a perfect accuracy. I will try to rearrange the algorithm.

walkietalkie74, Feb 2, 2014
8. ### JaxGuest

Walkie.... I'm a complete tyro at C++ so I'm not able to comment on your
examples.

Nice to see people like you still use Knuth/TAOCP although it's highly
theoretical and section 4.2.2 only contains a handful of theorems (four I
think) whose results are most probably incorporated into library subroutines.

Why are you using vector geometry for this? If your problem had dynamic
properties (such as aircraft dynamics in 3D) then I could understand. But
something as static as your straight line seems to be crying out for a
Cartesian approach. Just saying!

Jax, Feb 3, 2014
9. ### red floydGuest

red floyd, Feb 4, 2014