Re: floating point comparison

[Martin@nezumi.demon.co.uk (Martin Tom Brown)]

In article <mheaney-ya023680003007972240570001@news.ni.net>
           mheaney@ni.net "Matthew Heaney" writes:

> In article <dewar.870304869@merv>, dewar@merv.cs.nyu.edu (Robert Dewar) wrote:
> 
> 
> >As for when it makes sense to compare floating-point values for equality,
> >that is a matter for the programmer to understand. Floating-point arithmetic
> >on modern IEEE machines is not some kind of approximate hocus-pocus, it
> >is a well defined arithmetic systenm, with well defined, well behaved
> >results, in which equality has a perfectly reasonable meaning.
> 
> OK, so if  do this
> 
> declare
>    m1, m2 : slope;
> begin
>    if m1 = m2 then
> 
> and I'm using IEEE floating point types (I'm on a SUN box), can I assume
> that the lines are parallel?

There are two obvious situations where you can get into trouble -

#1 m1 equal to m2 	but the lines are not parallel, and numerical 
			cancellation gives fewer significant bits.

You can try and control this loss of significance.

#2 m1 not equal to m2 	but the lines are actually parallel and the
			small difference is due to rounding errors.

On some systems this may be unavoidable hence you have to accept
numbers as equal if they differ by less than some limit. 

If "slope" is a single floating point number then you are in trouble
when a line becomes parallel to one of the two axes (division by 0).

As Robert has said IEEE is a well defined floating-point number system
and so for any given calculation it is possible to determine bounds.
How much thought you give to a problem depends on how much it matters - 
for example there is probably no need to worry if a plotter pen is 
a few microns out of position after drawing a 1m line.

As a concrete and I think simpler example of the difference between
algebraic identities and numerical IEEE FP approximations try 
evaluating the following in the explicit order shown :

	f1(x) =	1 + x + x^2 + ... + x^(N-1) + x^N

	f2(x) = x^N + x^(N-1) + x^(N-2) + ... + x + 1

	f3(x) = (1 - x^(N+1))/(1-x)

Barring typos I hope everyone agrees these are algebraic identities
(at least for finite N)

If you now try to evaluate each of them numerically for say N=10 
and x = 10, 3, 2, 1, 1/2, 1/3, 1/10 you should see some differences.
(There is a gotcha in this list)

Common programmer errors often arise from exact decimal fractions 
which do not have a finite representation in the binary system.

Regards,
-- 
Martin Brown  <martin@nezumi.demon.co.uk>     __                CIS: 71651,470
Scientific Software Consultancy             /^,,)__/