Re: Infinite loop with GCC 4.3

[Joe Neeman]

On Mon, 2008-05-26 at 10:09 +0200, Mats Bengtsson wrote:
> As is taught in any course on numerical methods, you should always be 
> very careful with
> comparing numerical values. For example you should use tests of the form
> if (abs(x-y) < epsilon)
> instead of
> if (x==y)
> where epsilon is a fixed value somewhat larger than machine precisions.
> 
> Following the same line of thought, you should use
> if (x >= y - epsilon)
> instead of
> if (x >= y )
> (though this is probably not mentioned as often in the course text books).
> and analogously it follows that you numerically cannot distinguish between
> if (x > y)
> and
> if (x >= y)
> 
> Here you show an example where you cannot rely on
> if (x > y)
> but have to do something more or less equivalent to
> if ( x > y - epsilon)
> to guarantee numerical stability.

The problem in this case (and in skyline.cc) is not numerical stability
but the consitency of comparisons. In fact, the algorithm in skyline.cc
was carefully written so that its correctness does not depend on the
rounding of floating point numbers (this was after Han-Wen complained --
quite rightly IMO -- that my previous version of skyline.cc was
difficult to understand because of the proliferation of fuzzy
comparisons).

Instead of fuzzy comparisons, the skyline merging algorithm relies on
the fact that a certain double will strictly increase with every loop
iteration. This is relatively easy to guarantee. The problem is that it
doesn't work, since GCC on x86 allows
if (x > y)
  assert (x > y)
to fail.

I would really like to avoid going back to fuzzy comparisons because it
makes the logic very difficult (for example, a previous version had to
deal with the possibility that building 'a' is above building 'b' at
points 'x' and 'x + 1' but the intersection of their roofs occurs
between 'x' and 'x + 1').

Slightly off-topic, but there are certainly situations in which equality
of floating point numbers is important. For example, 0.0 / 0.0 is NaN
but epsilon / 0.0 is INFINITY. Similarly, if max always used fuzzy
comparisons, you could never guarantee that
pow(max(x, 0.0), 0.5)
would return something meaningful.

Joe