Floating point non-exactness

[Michael C. Grant]

In article <11117 at alice.UUCP>, ark at alice.UUCP (Andrew Koenig) writes:
> In article <5467 at quanta.eng.ohio-state.edu>, rob at baloo.eng.ohio-state.edu (Rob Carriere) writes:
>> In other words, usually what you want is something like
>> int fcmp( double a, double b, double eps )
>> {
>>    if (  fabs( a-b ) < eps  ) return 0;
>>    if (  a > b  ) return 1;
>>    return -1;
>> }
> Probably not.
> 
> The trouble is that if a and b are big enough, then a-b will either be
> exactly 0 or will be larger than eps.  (Here, `big enough' is approximately
> eps * 2^n, where n is the number of significant bits in a floating-point
> freaction.)  There is also the possibility that a-b might overflow,
> in which case the comparison would surely have been valid if done directly.

Perhaps a solution to this problem would be an implementation-
dependent library function fcmp, in which 'eps' is not some absolute
number, but a fraction of the two nearly-equal numbers themselves:

	if (fabs(a-b)<(eps*a)) return 0; 
	(assuming a is almost equal to b)

The reason it would need to be a library function is that implementing
it efficiently would involve extracting the mantissa and exponent and
examining them separately.  For example, assuming a standard which
requires normalization, this algorithm might work (note this is 
just a near-equality test, not a complete comparison):

	if (exponent_a!=exponent_b) return false
	if abs(mantissa_a-mantissa_b)<error_limit return true
	   else return false

Of course, if you don't care about portable code just write
this function yourself!

Michael C. Grant
no .sig, no .identity