Floating point non-exactness

[carlos at urz.unibas.ch]

In article <622 at .tetrauk.UUCP>, rick at .tetrauk.UUCP (Rick Jones) writes:
> 
> (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 at tetrauk.uucp (...!ukc!tetrauk.uucp!rick)	     - John Cougar Mellencamp