Path: utzoo!attcan!uunet!mcsun!ukc!pyrltd!tetrauk!rick
From: rick@.tetrauk.UUCP (Rick Jones)
Newsgroups: comp.lang.c
Subject: Floating point non-exactness
Message-ID: <622@.tetrauk.UUCP>
Date: 30 Jul 90 10:52:12 GMT
Reply-To: rick@tetrauk.UUCP (Rick Jones)
Organization: Tetra Ltd., Maidenhead, UK
Lines: 59


(This is nothing to do with unsigned integers :-)

I don't know if this has been discussed in the dim past, but I haven't seen
anything since I have been reading News.

The problem is well-known - not every rational decimal number has an exact
representation in binary floating-point form.  This leads to anomalies such
as:
	a = .291 / 2.0 ;
	b = .1455 ;

	if (a == b)		/* false! */

Both expressions as expressed in decimal yield exactly .1455, but the computed
results may not in fact be equal.  This fails on my system when the 3 values
are generated from strings using atof().  What does and does not fail will
presumably depend on the FP representation used, and possibly on the actual
hardware and compiler.

The reason is that .1455 does not have an exact binary representation.  In this
case the nearest two values are .14549999999999999 and .14550000000000002.  The
two routes to the result get different nearest approximations.  Thus the two
results are not actually equal.

Is there a general solution to this problem without continuously rounding
everything?  I believe it's possible by applying a "maximum reliable precision"
concept to FP numbers.  For IEEE 64-bit numbers this seems to be 1 in 10^16, or
15 significant digits.  Using this, a function "fpcmp" can be written which
takes 2 arguments:

	fpcmp (double a, double b)

		returns 0 if a equals b within the reliable precision,
		else returns 1 if a > b, or -1 if a < b

Real problem:  how do you write fpcmp() ?

I've got one or two possible solutions, but they're not very nice, and
computationally inefficient.  Has anyone got a really good answer to this?

Associated problems:

Can you ensure that "fpcmp (a, b)" and "fpcmp (a-b, 0.0)" yield the same result
when a and b are very close?

What about a rounding function which will consistently round on "5"?  E.g.
using (s)printf to convert .1455 to 3 dec. places will give .145 and .146 for
the two different close values above.


While it's a good discussion subject, I suggest that anything extensive, or
with large chunks of code, you email to me.  I'll summarise such contributions
later.

-- 
Rick Jones					You gotta stand for something
Tetra Ltd.  Maidenhead, Berks			Or you'll fall for anything
rick@tetrauk.uucp (...!ukc!tetrauk.uucp!rick)	     - John Cougar Mellencamp